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