001:       SUBROUTINE CPOTF2( UPLO, N, A, LDA, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          UPLO
009:       INTEGER            INFO, LDA, N
010: *     ..
011: *     .. Array Arguments ..
012:       COMPLEX            A( LDA, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  CPOTF2 computes the Cholesky factorization of a complex Hermitian
019: *  positive definite matrix A.
020: *
021: *  The factorization has the form
022: *     A = U' * U ,  if UPLO = 'U', or
023: *     A = L  * L',  if UPLO = 'L',
024: *  where U is an upper triangular matrix and L is lower triangular.
025: *
026: *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  UPLO    (input) CHARACTER*1
032: *          Specifies whether the upper or lower triangular part of the
033: *          Hermitian matrix A is stored.
034: *          = 'U':  Upper triangular
035: *          = 'L':  Lower triangular
036: *
037: *  N       (input) INTEGER
038: *          The order of the matrix A.  N >= 0.
039: *
040: *  A       (input/output) COMPLEX array, dimension (LDA,N)
041: *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
042: *          n by n upper triangular part of A contains the upper
043: *          triangular part of the matrix A, and the strictly lower
044: *          triangular part of A is not referenced.  If UPLO = 'L', the
045: *          leading n by n lower triangular part of A contains the lower
046: *          triangular part of the matrix A, and the strictly upper
047: *          triangular part of A is not referenced.
048: *
049: *          On exit, if INFO = 0, the factor U or L from the Cholesky
050: *          factorization A = U'*U  or A = L*L'.
051: *
052: *  LDA     (input) INTEGER
053: *          The leading dimension of the array A.  LDA >= max(1,N).
054: *
055: *  INFO    (output) INTEGER
056: *          = 0: successful exit
057: *          < 0: if INFO = -k, the k-th argument had an illegal value
058: *          > 0: if INFO = k, the leading minor of order k is not
059: *               positive definite, and the factorization could not be
060: *               completed.
061: *
062: *  =====================================================================
063: *
064: *     .. Parameters ..
065:       REAL               ONE, ZERO
066:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
067:       COMPLEX            CONE
068:       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
069: *     ..
070: *     .. Local Scalars ..
071:       LOGICAL            UPPER
072:       INTEGER            J
073:       REAL               AJJ
074: *     ..
075: *     .. External Functions ..
076:       LOGICAL            LSAME, SISNAN
077:       COMPLEX            CDOTC
078:       EXTERNAL           LSAME, CDOTC, SISNAN
079: *     ..
080: *     .. External Subroutines ..
081:       EXTERNAL           CGEMV, CLACGV, CSSCAL, XERBLA
082: *     ..
083: *     .. Intrinsic Functions ..
084:       INTRINSIC          MAX, REAL, SQRT
085: *     ..
086: *     .. Executable Statements ..
087: *
088: *     Test the input parameters.
089: *
090:       INFO = 0
091:       UPPER = LSAME( UPLO, 'U' )
092:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
093:          INFO = -1
094:       ELSE IF( N.LT.0 ) THEN
095:          INFO = -2
096:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
097:          INFO = -4
098:       END IF
099:       IF( INFO.NE.0 ) THEN
100:          CALL XERBLA( 'CPOTF2', -INFO )
101:          RETURN
102:       END IF
103: *
104: *     Quick return if possible
105: *
106:       IF( N.EQ.0 )
107:      $   RETURN
108: *
109:       IF( UPPER ) THEN
110: *
111: *        Compute the Cholesky factorization A = U'*U.
112: *
113:          DO 10 J = 1, N
114: *
115: *           Compute U(J,J) and test for non-positive-definiteness.
116: *
117:             AJJ = REAL( A( J, J ) ) - CDOTC( J-1, A( 1, J ), 1,
118:      $            A( 1, J ), 1 )
119:             IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
120:                A( J, J ) = AJJ
121:                GO TO 30
122:             END IF
123:             AJJ = SQRT( AJJ )
124:             A( J, J ) = AJJ
125: *
126: *           Compute elements J+1:N of row J.
127: *
128:             IF( J.LT.N ) THEN
129:                CALL CLACGV( J-1, A( 1, J ), 1 )
130:                CALL CGEMV( 'Transpose', J-1, N-J, -CONE, A( 1, J+1 ),
131:      $                     LDA, A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
132:                CALL CLACGV( J-1, A( 1, J ), 1 )
133:                CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
134:             END IF
135:    10    CONTINUE
136:       ELSE
137: *
138: *        Compute the Cholesky factorization A = L*L'.
139: *
140:          DO 20 J = 1, N
141: *
142: *           Compute L(J,J) and test for non-positive-definiteness.
143: *
144:             AJJ = REAL( A( J, J ) ) - CDOTC( J-1, A( J, 1 ), LDA,
145:      $            A( J, 1 ), LDA )
146:             IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
147:                A( J, J ) = AJJ
148:                GO TO 30
149:             END IF
150:             AJJ = SQRT( AJJ )
151:             A( J, J ) = AJJ
152: *
153: *           Compute elements J+1:N of column J.
154: *
155:             IF( J.LT.N ) THEN
156:                CALL CLACGV( J-1, A( J, 1 ), LDA )
157:                CALL CGEMV( 'No transpose', N-J, J-1, -CONE, A( J+1, 1 ),
158:      $                     LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
159:                CALL CLACGV( J-1, A( J, 1 ), LDA )
160:                CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
161:             END IF
162:    20    CONTINUE
163:       END IF
164:       GO TO 40
165: *
166:    30 CONTINUE
167:       INFO = J
168: *
169:    40 CONTINUE
170:       RETURN
171: *
172: *     End of CPOTF2
173: *
174:       END
175: