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:       COMPLEX            TMP
072:       LOGICAL            UPPER
073:       INTEGER            J
074:       REAL               AJJ
075: *     ..
076: *     .. External Functions ..
077:       LOGICAL            LSAME, SISNAN
078: *      COMPLEX            CDOTC
079: *      EXTERNAL           LSAME, CDOTC, SISNAN
080:       EXTERNAL           LSAME, SISNAN
081: *     ..
082: *     .. External Subroutines ..
083:       EXTERNAL           CGEMV, CLACGV, CSSCAL, XERBLA
084: *     ..
085: *     .. Intrinsic Functions ..
086:       INTRINSIC          MAX, REAL, SQRT
087: *     ..
088: *     .. Executable Statements ..
089: *
090: *     Test the input parameters.
091: *
092:       INFO = 0
093:       UPPER = LSAME( UPLO, 'U' )
094:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
095:          INFO = -1
096:       ELSE IF( N.LT.0 ) THEN
097:          INFO = -2
098:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
099:          INFO = -4
100:       END IF
101:       IF( INFO.NE.0 ) THEN
102:          CALL XERBLA( 'CPOTF2', -INFO )
103:          RETURN
104:       END IF
105: *
106: *     Quick return if possible
107: *
108:       IF( N.EQ.0 )
109:      $   RETURN
110: *
111:       IF( UPPER ) THEN
112: *
113: *        Compute the Cholesky factorization A = U'*U.
114: *
115:          DO 10 J = 1, N
116: *
117: *           Compute U(J,J) and test for non-positive-definiteness.
118: *
119:             CALL PLASMA_CDOTCSUB( J-1, A( 1, J ), 1, A( 1, J ), 1, 
120:      $                           TMP )
121:             AJJ = REAL( A( J, J ) ) - TMP
122:             IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
123:                A( J, J ) = AJJ
124:                GO TO 30
125:             END IF
126:             AJJ = SQRT( AJJ )
127:             A( J, J ) = AJJ
128: *
129: *           Compute elements J+1:N of row J.
130: *
131:             IF( J.LT.N ) THEN
132:                CALL CLACGV( J-1, A( 1, J ), 1 )
133:                CALL CGEMV( 'Transpose', J-1, N-J, -CONE, A( 1, J+1 ),
134:      $                     LDA, A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
135:                CALL CLACGV( J-1, A( 1, J ), 1 )
136:                CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
137:             END IF
138:    10    CONTINUE
139:       ELSE
140: *
141: *        Compute the Cholesky factorization A = L*L'.
142: *
143:          DO 20 J = 1, N
144: *
145: *           Compute L(J,J) and test for non-positive-definiteness.
146: *
147:             CALL PLASMA_CDOTCSUB( J-1, A( J, 1 ), LDA, A( J, 1 ), 
148:      $                           LDA, TMP )
149:             AJJ = REAL( A( J, J ) ) - TMP
150:             IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
151:                A( J, J ) = AJJ
152:                GO TO 30
153:             END IF
154:             AJJ = SQRT( AJJ )
155:             A( J, J ) = AJJ
156: *
157: *           Compute elements J+1:N of column J.
158: *
159:             IF( J.LT.N ) THEN
160:                CALL CLACGV( J-1, A( J, 1 ), LDA )
161:                CALL CGEMV( 'No transpose', N-J, J-1, -CONE, A( J+1, 1 ),
162:      $                     LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
163:                CALL CLACGV( J-1, A( J, 1 ), LDA )
164:                CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
165:             END IF
166:    20    CONTINUE
167:       END IF
168:       GO TO 40
169: *
170:    30 CONTINUE
171:       INFO = J
172: *
173:    40 CONTINUE
174:       RETURN
175: *
176: *     End of CPOTF2
177: *
178:       END
179: