1 SUBROUTINE CPOTRF( UPLO, N, A, LDA, INFO ) 2* 3* -- LAPACK routine (version 3.0) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* September 30, 1994 7* 8* .. Scalar Arguments .. 9 CHARACTER UPLO 10 INTEGER INFO, LDA, N 11* .. 12* .. Array Arguments .. 13 COMPLEX A( LDA, * ) 14* .. 15* 16* Purpose 17* ======= 18* 19* CPOTRF computes the Cholesky factorization of a complex Hermitian 20* positive definite matrix A. 21* 22* The factorization has the form 23* A = U**H * U, if UPLO = 'U', or 24* A = L * L**H, if UPLO = 'L', 25* where U is an upper triangular matrix and L is lower triangular. 26* 27* This is the block version of the algorithm, calling Level 3 BLAS. 28* 29* Arguments 30* ========= 31* 32* UPLO (input) CHARACTER*1 33* = 'U': Upper triangle of A is stored; 34* = 'L': Lower triangle of A is stored. 35* 36* N (input) INTEGER 37* The order of the matrix A. N >= 0. 38* 39* A (input/output) COMPLEX array, dimension (LDA,N) 40* On entry, the Hermitian matrix A. If UPLO = 'U', the leading 41* N-by-N upper triangular part of A contains the upper 42* triangular part of the matrix A, and the strictly lower 43* triangular part of A is not referenced. If UPLO = 'L', the 44* leading N-by-N lower triangular part of A contains the lower 45* triangular part of the matrix A, and the strictly upper 46* triangular part of A is not referenced. 47* 48* On exit, if INFO = 0, the factor U or L from the Cholesky 49* factorization A = U**H*U or A = L*L**H. 50* 51* LDA (input) INTEGER 52* The leading dimension of the array A. LDA >= max(1,N). 53* 54* INFO (output) INTEGER 55* = 0: successful exit 56* < 0: if INFO = -i, the i-th argument had an illegal value 57* > 0: if INFO = i, the leading minor of order i is not 58* positive definite, and the factorization could not be 59* completed. 60* 61* ===================================================================== 62* 63* .. Parameters .. 64 REAL ONE 65 COMPLEX CONE 66 PARAMETER ( ONE = 1.0E+0, CONE = ( 1.0E+0, 0.0E+0 ) ) 67* .. 68* .. Local Scalars .. 69 LOGICAL UPPER 70 INTEGER J, JB, NB 71* .. 72* .. External Functions .. 73 LOGICAL LSAME 74 INTEGER ILAENV 75 EXTERNAL LSAME, ILAENV 76* .. 77* .. External Subroutines .. 78 EXTERNAL CGEMM, CHERK, CPOTF2, CTRSM, XERBLA 79* .. 80* .. Intrinsic Functions .. 81 INTRINSIC MAX, MIN 82* .. 83* .. Executable Statements .. 84* 85* Test the input parameters. 86* 87 INFO = 0 88 UPPER = LSAME( UPLO, 'U' ) 89 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 90 INFO = -1 91 ELSE IF( N.LT.0 ) THEN 92 INFO = -2 93 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 94 INFO = -4 95 END IF 96 IF( INFO.NE.0 ) THEN 97 CALL XERBLA( 'CPOTRF', -INFO ) 98 RETURN 99 END IF 100* 101* Quick return if possible 102* 103 IF( N.EQ.0 ) 104 $ RETURN 105* 106* Determine the block size for this environment. 107* 108 NB = ILAENV( 1, 'CPOTRF', UPLO, N, -1, -1, -1 ) 109 IF( NB.LE.1 .OR. NB.GE.N ) THEN 110* 111* Use unblocked code. 112* 113 CALL CPOTF2( UPLO, N, A, LDA, INFO ) 114 ELSE 115* 116* Use blocked code. 117* 118 IF( UPPER ) THEN 119* 120* Compute the Cholesky factorization A = U'*U. 121* 122 DO 10 J = 1, N, NB 123* 124* Update and factorize the current diagonal block and test 125* for non-positive-definiteness. 126* 127 JB = MIN( NB, N-J+1 ) 128 CALL CHERK( 'Upper', 'Conjugate transpose', JB, J-1, 129 $ -ONE, A( 1, J ), LDA, ONE, A( J, J ), LDA ) 130 CALL CPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) 131 IF( INFO.NE.0 ) 132 $ GO TO 30 133 IF( J+JB.LE.N ) THEN 134* 135* Compute the current block row. 136* 137 CALL CGEMM( 'Conjugate transpose', 'No transpose', JB, 138 $ N-J-JB+1, J-1, -CONE, A( 1, J ), LDA, 139 $ A( 1, J+JB ), LDA, CONE, A( J, J+JB ), 140 $ LDA ) 141 CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose', 142 $ 'Non-unit', JB, N-J-JB+1, CONE, A( J, J ), 143 $ LDA, A( J, J+JB ), LDA ) 144 END IF 145 10 CONTINUE 146* 147 ELSE 148* 149* Compute the Cholesky factorization A = L*L'. 150* 151 DO 20 J = 1, N, NB 152* 153* Update and factorize the current diagonal block and test 154* for non-positive-definiteness. 155* 156 JB = MIN( NB, N-J+1 ) 157 CALL CHERK( 'Lower', 'No transpose', JB, J-1, -ONE, 158 $ A( J, 1 ), LDA, ONE, A( J, J ), LDA ) 159 CALL CPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) 160 IF( INFO.NE.0 ) 161 $ GO TO 30 162 IF( J+JB.LE.N ) THEN 163* 164* Compute the current block column. 165* 166 CALL CGEMM( 'No transpose', 'Conjugate transpose', 167 $ N-J-JB+1, JB, J-1, -CONE, A( J+JB, 1 ), 168 $ LDA, A( J, 1 ), LDA, CONE, A( J+JB, J ), 169 $ LDA ) 170 CALL CTRSM( 'Right', 'Lower', 'Conjugate transpose', 171 $ 'Non-unit', N-J-JB+1, JB, CONE, A( J, J ), 172 $ LDA, A( J+JB, J ), LDA ) 173 END IF 174 20 CONTINUE 175 END IF 176 END IF 177 GO TO 40 178* 179 30 CONTINUE 180 INFO = INFO + J - 1 181* 182 40 CONTINUE 183 RETURN 184* 185* End of CPOTRF 186* 187 END 188