00001 SUBROUTINE SPOTF2( UPLO, N, A, LDA, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER INFO, LDA, N 00011 * .. 00012 * .. Array Arguments .. 00013 REAL A( LDA, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * SPOTF2 computes the Cholesky factorization of a real symmetric 00020 * positive definite matrix A. 00021 * 00022 * The factorization has the form 00023 * A = U' * U , if UPLO = 'U', or 00024 * A = L * L', if UPLO = 'L', 00025 * where U is an upper triangular matrix and L is lower triangular. 00026 * 00027 * This is the unblocked version of the algorithm, calling Level 2 BLAS. 00028 * 00029 * Arguments 00030 * ========= 00031 * 00032 * UPLO (input) CHARACTER*1 00033 * Specifies whether the upper or lower triangular part of the 00034 * symmetric matrix A is stored. 00035 * = 'U': Upper triangular 00036 * = 'L': Lower triangular 00037 * 00038 * N (input) INTEGER 00039 * The order of the matrix A. N >= 0. 00040 * 00041 * A (input/output) REAL array, dimension (LDA,N) 00042 * On entry, the symmetric matrix A. If UPLO = 'U', the leading 00043 * n by n upper triangular part of A contains the upper 00044 * triangular part of the matrix A, and the strictly lower 00045 * triangular part of A is not referenced. If UPLO = 'L', the 00046 * leading n by n lower triangular part of A contains the lower 00047 * triangular part of the matrix A, and the strictly upper 00048 * triangular part of A is not referenced. 00049 * 00050 * On exit, if INFO = 0, the factor U or L from the Cholesky 00051 * factorization A = U'*U or A = L*L'. 00052 * 00053 * LDA (input) INTEGER 00054 * The leading dimension of the array A. LDA >= max(1,N). 00055 * 00056 * INFO (output) INTEGER 00057 * = 0: successful exit 00058 * < 0: if INFO = -k, the k-th argument had an illegal value 00059 * > 0: if INFO = k, the leading minor of order k is not 00060 * positive definite, and the factorization could not be 00061 * completed. 00062 * 00063 * ===================================================================== 00064 * 00065 * .. Parameters .. 00066 REAL ONE, ZERO 00067 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00068 * .. 00069 * .. Local Scalars .. 00070 LOGICAL UPPER 00071 INTEGER J 00072 REAL AJJ 00073 * .. 00074 * .. External Functions .. 00075 LOGICAL LSAME, SISNAN 00076 REAL SDOT 00077 EXTERNAL LSAME, SDOT, SISNAN 00078 * .. 00079 * .. External Subroutines .. 00080 EXTERNAL SGEMV, SSCAL, XERBLA 00081 * .. 00082 * .. Intrinsic Functions .. 00083 INTRINSIC MAX, SQRT 00084 * .. 00085 * .. Executable Statements .. 00086 * 00087 * Test the input parameters. 00088 * 00089 INFO = 0 00090 UPPER = LSAME( UPLO, 'U' ) 00091 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00092 INFO = -1 00093 ELSE IF( N.LT.0 ) THEN 00094 INFO = -2 00095 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00096 INFO = -4 00097 END IF 00098 IF( INFO.NE.0 ) THEN 00099 CALL XERBLA( 'SPOTF2', -INFO ) 00100 RETURN 00101 END IF 00102 * 00103 * Quick return if possible 00104 * 00105 IF( N.EQ.0 ) 00106 $ RETURN 00107 * 00108 IF( UPPER ) THEN 00109 * 00110 * Compute the Cholesky factorization A = U'*U. 00111 * 00112 DO 10 J = 1, N 00113 * 00114 * Compute U(J,J) and test for non-positive-definiteness. 00115 * 00116 AJJ = A( J, J ) - SDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 ) 00117 IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN 00118 A( J, J ) = AJJ 00119 GO TO 30 00120 END IF 00121 AJJ = SQRT( AJJ ) 00122 A( J, J ) = AJJ 00123 * 00124 * Compute elements J+1:N of row J. 00125 * 00126 IF( J.LT.N ) THEN 00127 CALL SGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ), 00128 $ LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA ) 00129 CALL SSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 00130 END IF 00131 10 CONTINUE 00132 ELSE 00133 * 00134 * Compute the Cholesky factorization A = L*L'. 00135 * 00136 DO 20 J = 1, N 00137 * 00138 * Compute L(J,J) and test for non-positive-definiteness. 00139 * 00140 AJJ = A( J, J ) - SDOT( J-1, A( J, 1 ), LDA, A( J, 1 ), 00141 $ LDA ) 00142 IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN 00143 A( J, J ) = AJJ 00144 GO TO 30 00145 END IF 00146 AJJ = SQRT( AJJ ) 00147 A( J, J ) = AJJ 00148 * 00149 * Compute elements J+1:N of column J. 00150 * 00151 IF( J.LT.N ) THEN 00152 CALL SGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ), 00153 $ LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 ) 00154 CALL SSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 00155 END IF 00156 20 CONTINUE 00157 END IF 00158 GO TO 40 00159 * 00160 30 CONTINUE 00161 INFO = J 00162 * 00163 40 CONTINUE 00164 RETURN 00165 * 00166 * End of SPOTF2 00167 * 00168 END
1.7.2