LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SPOTRF( UPLO, N, A, LDA, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * -- April 2011 -- 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 * SPOTRF computes the Cholesky factorization of a real symmetric 00020 * positive definite matrix A. 00021 * 00022 * The factorization has the form 00023 * A = U**T * U, if UPLO = 'U', or 00024 * A = L * L**T, if UPLO = 'L', 00025 * where U is an upper triangular matrix and L is lower triangular. 00026 * 00027 * This is the block version of the algorithm, calling Level 3 BLAS. 00028 * 00029 * Arguments 00030 * ========= 00031 * 00032 * UPLO (input) CHARACTER*1 00033 * = 'U': Upper triangle of A is stored; 00034 * = 'L': Lower triangle of A is stored. 00035 * 00036 * N (input) INTEGER 00037 * The order of the matrix A. N >= 0. 00038 * 00039 * A (input/output) REAL array, dimension (LDA,N) 00040 * On entry, the symmetric matrix A. If UPLO = 'U', the leading 00041 * N-by-N upper triangular part of A contains the upper 00042 * triangular part of the matrix A, and the strictly lower 00043 * triangular part of A is not referenced. If UPLO = 'L', the 00044 * leading N-by-N lower triangular part of A contains the lower 00045 * triangular part of the matrix A, and the strictly upper 00046 * triangular part of A is not referenced. 00047 * 00048 * On exit, if INFO = 0, the factor U or L from the Cholesky 00049 * factorization A = U**T*U or A = L*L**T. 00050 * 00051 * LDA (input) INTEGER 00052 * The leading dimension of the array A. LDA >= max(1,N). 00053 * 00054 * INFO (output) INTEGER 00055 * = 0: successful exit 00056 * < 0: if INFO = -i, the i-th argument had an illegal value 00057 * > 0: if INFO = i, the leading minor of order i is not 00058 * positive definite, and the factorization could not be 00059 * completed. 00060 * 00061 * ===================================================================== 00062 * 00063 * .. Parameters .. 00064 REAL ONE 00065 PARAMETER ( ONE = 1.0E+0 ) 00066 * .. 00067 * .. Local Scalars .. 00068 LOGICAL UPPER 00069 INTEGER J, JB, NB 00070 * .. 00071 * .. External Functions .. 00072 LOGICAL LSAME 00073 INTEGER ILAENV 00074 EXTERNAL LSAME, ILAENV 00075 * .. 00076 * .. External Subroutines .. 00077 EXTERNAL SGEMM, SPOTF2, SSYRK, STRSM, XERBLA 00078 * .. 00079 * .. Intrinsic Functions .. 00080 INTRINSIC MAX, MIN 00081 * .. 00082 * .. Executable Statements .. 00083 * 00084 * Test the input parameters. 00085 * 00086 INFO = 0 00087 UPPER = LSAME( UPLO, 'U' ) 00088 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00089 INFO = -1 00090 ELSE IF( N.LT.0 ) THEN 00091 INFO = -2 00092 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00093 INFO = -4 00094 END IF 00095 IF( INFO.NE.0 ) THEN 00096 CALL XERBLA( 'SPOTRF', -INFO ) 00097 RETURN 00098 END IF 00099 * 00100 * Quick return if possible 00101 * 00102 IF( N.EQ.0 ) 00103 $ RETURN 00104 * 00105 * Determine the block size for this environment. 00106 * 00107 NB = ILAENV( 1, 'SPOTRF', UPLO, N, -1, -1, -1 ) 00108 IF( NB.LE.1 .OR. NB.GE.N ) THEN 00109 * 00110 * Use unblocked code. 00111 * 00112 CALL SPOTF2( UPLO, N, A, LDA, INFO ) 00113 ELSE 00114 * 00115 * Use blocked code. 00116 * 00117 IF( UPPER ) THEN 00118 * 00119 * Compute the Cholesky factorization A = U**T*U. 00120 * 00121 DO 10 J = 1, N, NB 00122 * 00123 * Update and factorize the current diagonal block and test 00124 * for non-positive-definiteness. 00125 * 00126 JB = MIN( NB, N-J+1 ) 00127 CALL SSYRK( 'Upper', 'Transpose', JB, J-1, -ONE, 00128 $ A( 1, J ), LDA, ONE, A( J, J ), LDA ) 00129 CALL SPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) 00130 IF( INFO.NE.0 ) 00131 $ GO TO 30 00132 IF( J+JB.LE.N ) THEN 00133 * 00134 * Compute the current block row. 00135 * 00136 CALL SGEMM( 'Transpose', 'No transpose', JB, N-J-JB+1, 00137 $ J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ), 00138 $ LDA, ONE, A( J, J+JB ), LDA ) 00139 CALL STRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', 00140 $ JB, N-J-JB+1, ONE, A( J, J ), LDA, 00141 $ A( J, J+JB ), LDA ) 00142 END IF 00143 10 CONTINUE 00144 * 00145 ELSE 00146 * 00147 * Compute the Cholesky factorization A = L*L**T. 00148 * 00149 DO 20 J = 1, N, NB 00150 * 00151 * Update and factorize the current diagonal block and test 00152 * for non-positive-definiteness. 00153 * 00154 JB = MIN( NB, N-J+1 ) 00155 CALL SSYRK( 'Lower', 'No transpose', JB, J-1, -ONE, 00156 $ A( J, 1 ), LDA, ONE, A( J, J ), LDA ) 00157 CALL SPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) 00158 IF( INFO.NE.0 ) 00159 $ GO TO 30 00160 IF( J+JB.LE.N ) THEN 00161 * 00162 * Compute the current block column. 00163 * 00164 CALL SGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, 00165 $ J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ), 00166 $ LDA, ONE, A( J+JB, J ), LDA ) 00167 CALL STRSM( 'Right', 'Lower', 'Transpose', 'Non-unit', 00168 $ N-J-JB+1, JB, ONE, A( J, J ), LDA, 00169 $ A( J+JB, J ), LDA ) 00170 END IF 00171 20 CONTINUE 00172 END IF 00173 END IF 00174 GO TO 40 00175 * 00176 30 CONTINUE 00177 INFO = INFO + J - 1 00178 * 00179 40 CONTINUE 00180 RETURN 00181 * 00182 * End of SPOTRF 00183 * 00184 END