LAPACK 3.3.0
|
00001 SUBROUTINE CPOTF2( 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 COMPLEX A( LDA, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CPOTF2 computes the Cholesky factorization of a complex Hermitian 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 * Hermitian 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) COMPLEX array, dimension (LDA,N) 00042 * On entry, the Hermitian 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 COMPLEX CONE 00069 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00070 * .. 00071 * .. Local Scalars .. 00072 LOGICAL UPPER 00073 INTEGER J 00074 REAL AJJ 00075 * .. 00076 * .. External Functions .. 00077 LOGICAL LSAME, SISNAN 00078 COMPLEX CDOTC 00079 EXTERNAL LSAME, CDOTC, SISNAN 00080 * .. 00081 * .. External Subroutines .. 00082 EXTERNAL CGEMV, CLACGV, CSSCAL, XERBLA 00083 * .. 00084 * .. Intrinsic Functions .. 00085 INTRINSIC MAX, REAL, SQRT 00086 * .. 00087 * .. Executable Statements .. 00088 * 00089 * Test the input parameters. 00090 * 00091 INFO = 0 00092 UPPER = LSAME( UPLO, 'U' ) 00093 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00094 INFO = -1 00095 ELSE IF( N.LT.0 ) THEN 00096 INFO = -2 00097 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00098 INFO = -4 00099 END IF 00100 IF( INFO.NE.0 ) THEN 00101 CALL XERBLA( 'CPOTF2', -INFO ) 00102 RETURN 00103 END IF 00104 * 00105 * Quick return if possible 00106 * 00107 IF( N.EQ.0 ) 00108 $ RETURN 00109 * 00110 IF( UPPER ) THEN 00111 * 00112 * Compute the Cholesky factorization A = U'*U. 00113 * 00114 DO 10 J = 1, N 00115 * 00116 * Compute U(J,J) and test for non-positive-definiteness. 00117 * 00118 AJJ = REAL( A( J, J ) ) - CDOTC( J-1, A( 1, J ), 1, 00119 $ A( 1, J ), 1 ) 00120 IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN 00121 A( J, J ) = AJJ 00122 GO TO 30 00123 END IF 00124 AJJ = SQRT( AJJ ) 00125 A( J, J ) = AJJ 00126 * 00127 * Compute elements J+1:N of row J. 00128 * 00129 IF( J.LT.N ) THEN 00130 CALL CLACGV( J-1, A( 1, J ), 1 ) 00131 CALL CGEMV( 'Transpose', J-1, N-J, -CONE, A( 1, J+1 ), 00132 $ LDA, A( 1, J ), 1, CONE, A( J, J+1 ), LDA ) 00133 CALL CLACGV( J-1, A( 1, J ), 1 ) 00134 CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 00135 END IF 00136 10 CONTINUE 00137 ELSE 00138 * 00139 * Compute the Cholesky factorization A = L*L'. 00140 * 00141 DO 20 J = 1, N 00142 * 00143 * Compute L(J,J) and test for non-positive-definiteness. 00144 * 00145 AJJ = REAL( A( J, J ) ) - CDOTC( J-1, A( J, 1 ), LDA, 00146 $ A( J, 1 ), LDA ) 00147 IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN 00148 A( J, J ) = AJJ 00149 GO TO 30 00150 END IF 00151 AJJ = SQRT( AJJ ) 00152 A( J, J ) = AJJ 00153 * 00154 * Compute elements J+1:N of column J. 00155 * 00156 IF( J.LT.N ) THEN 00157 CALL CLACGV( J-1, A( J, 1 ), LDA ) 00158 CALL CGEMV( 'No transpose', N-J, J-1, -CONE, A( J+1, 1 ), 00159 $ LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 ) 00160 CALL CLACGV( J-1, A( J, 1 ), LDA ) 00161 CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 00162 END IF 00163 20 CONTINUE 00164 END IF 00165 GO TO 40 00166 * 00167 30 CONTINUE 00168 INFO = J 00169 * 00170 40 CONTINUE 00171 RETURN 00172 * 00173 * End of CPOTF2 00174 * 00175 END