00001 SUBROUTINE CPPTRF( UPLO, N, AP, 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, N 00011 * .. 00012 * .. Array Arguments .. 00013 COMPLEX AP( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CPPTRF computes the Cholesky factorization of a complex Hermitian 00020 * positive definite matrix A stored in packed format. 00021 * 00022 * The factorization has the form 00023 * A = U**H * U, if UPLO = 'U', or 00024 * A = L * L**H, if UPLO = 'L', 00025 * where U is an upper triangular matrix and L is lower triangular. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * UPLO (input) CHARACTER*1 00031 * = 'U': Upper triangle of A is stored; 00032 * = 'L': Lower triangle of A is stored. 00033 * 00034 * N (input) INTEGER 00035 * The order of the matrix A. N >= 0. 00036 * 00037 * AP (input/output) COMPLEX array, dimension (N*(N+1)/2) 00038 * On entry, the upper or lower triangle of the Hermitian matrix 00039 * A, packed columnwise in a linear array. The j-th column of A 00040 * is stored in the array AP as follows: 00041 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00042 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 00043 * See below for further details. 00044 * 00045 * On exit, if INFO = 0, the triangular factor U or L from the 00046 * Cholesky factorization A = U**H*U or A = L*L**H, in the same 00047 * storage format as A. 00048 * 00049 * INFO (output) INTEGER 00050 * = 0: successful exit 00051 * < 0: if INFO = -i, the i-th argument had an illegal value 00052 * > 0: if INFO = i, the leading minor of order i is not 00053 * positive definite, and the factorization could not be 00054 * completed. 00055 * 00056 * Further Details 00057 * =============== 00058 * 00059 * The packed storage scheme is illustrated by the following example 00060 * when N = 4, UPLO = 'U': 00061 * 00062 * Two-dimensional storage of the Hermitian matrix A: 00063 * 00064 * a11 a12 a13 a14 00065 * a22 a23 a24 00066 * a33 a34 (aij = conjg(aji)) 00067 * a44 00068 * 00069 * Packed storage of the upper triangle of A: 00070 * 00071 * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] 00072 * 00073 * ===================================================================== 00074 * 00075 * .. Parameters .. 00076 REAL ZERO, ONE 00077 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00078 * .. 00079 * .. Local Scalars .. 00080 LOGICAL UPPER 00081 INTEGER J, JC, JJ 00082 REAL AJJ 00083 * .. 00084 * .. External Functions .. 00085 LOGICAL LSAME 00086 COMPLEX CDOTC 00087 EXTERNAL LSAME, CDOTC 00088 * .. 00089 * .. External Subroutines .. 00090 EXTERNAL CHPR, CSSCAL, CTPSV, XERBLA 00091 * .. 00092 * .. Intrinsic Functions .. 00093 INTRINSIC REAL, SQRT 00094 * .. 00095 * .. Executable Statements .. 00096 * 00097 * Test the input parameters. 00098 * 00099 INFO = 0 00100 UPPER = LSAME( UPLO, 'U' ) 00101 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00102 INFO = -1 00103 ELSE IF( N.LT.0 ) THEN 00104 INFO = -2 00105 END IF 00106 IF( INFO.NE.0 ) THEN 00107 CALL XERBLA( 'CPPTRF', -INFO ) 00108 RETURN 00109 END IF 00110 * 00111 * Quick return if possible 00112 * 00113 IF( N.EQ.0 ) 00114 $ RETURN 00115 * 00116 IF( UPPER ) THEN 00117 * 00118 * Compute the Cholesky factorization A = U'*U. 00119 * 00120 JJ = 0 00121 DO 10 J = 1, N 00122 JC = JJ + 1 00123 JJ = JJ + J 00124 * 00125 * Compute elements 1:J-1 of column J. 00126 * 00127 IF( J.GT.1 ) 00128 $ CALL CTPSV( 'Upper', 'Conjugate transpose', 'Non-unit', 00129 $ J-1, AP, AP( JC ), 1 ) 00130 * 00131 * Compute U(J,J) and test for non-positive-definiteness. 00132 * 00133 AJJ = REAL( AP( JJ ) ) - CDOTC( J-1, AP( JC ), 1, AP( JC ), 00134 $ 1 ) 00135 IF( AJJ.LE.ZERO ) THEN 00136 AP( JJ ) = AJJ 00137 GO TO 30 00138 END IF 00139 AP( JJ ) = SQRT( AJJ ) 00140 10 CONTINUE 00141 ELSE 00142 * 00143 * Compute the Cholesky factorization A = L*L'. 00144 * 00145 JJ = 1 00146 DO 20 J = 1, N 00147 * 00148 * Compute L(J,J) and test for non-positive-definiteness. 00149 * 00150 AJJ = REAL( AP( JJ ) ) 00151 IF( AJJ.LE.ZERO ) THEN 00152 AP( JJ ) = AJJ 00153 GO TO 30 00154 END IF 00155 AJJ = SQRT( AJJ ) 00156 AP( JJ ) = AJJ 00157 * 00158 * Compute elements J+1:N of column J and update the trailing 00159 * submatrix. 00160 * 00161 IF( J.LT.N ) THEN 00162 CALL CSSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 ) 00163 CALL CHPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1, 00164 $ AP( JJ+N-J+1 ) ) 00165 JJ = JJ + N - J + 1 00166 END IF 00167 20 CONTINUE 00168 END IF 00169 GO TO 40 00170 * 00171 30 CONTINUE 00172 INFO = J 00173 * 00174 40 CONTINUE 00175 RETURN 00176 * 00177 * End of CPPTRF 00178 * 00179 END
1.7.2