 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2.2) -- 00004 * 00005 * -- Contributed by Craig Lucas, University of Manchester / NAG Ltd. -- 00006 * -- June 2010 -- 00007 * 00008 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00009 * 00010 * .. Scalar Arguments .. 00011 DOUBLE PRECISION TOL 00012 INTEGER INFO, LDA, N, RANK 00013 CHARACTER UPLO 00014 * .. 00015 * .. Array Arguments .. 00016 COMPLEX*16 A( LDA, * ) 00017 DOUBLE PRECISION WORK( 2*N ) 00018 INTEGER PIV( N ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * ZPSTRF computes the Cholesky factorization with complete 00025 * pivoting of a complex Hermitian positive semidefinite matrix A. 00026 * 00027 * The factorization has the form 00028 * P**T * A * P = U**H * U , if UPLO = 'U', 00029 * P**T * A * P = L * L**H, if UPLO = 'L', 00030 * where U is an upper triangular matrix and L is lower triangular, and 00031 * P is stored as vector PIV. 00032 * 00033 * This algorithm does not attempt to check that A is positive 00034 * semidefinite. This version of the algorithm calls level 3 BLAS. 00035 * 00036 * Arguments 00037 * ========= 00038 * 00039 * UPLO (input) CHARACTER*1 00040 * Specifies whether the upper or lower triangular part of the 00041 * symmetric matrix A is stored. 00042 * = 'U': Upper triangular 00043 * = 'L': Lower triangular 00044 * 00045 * N (input) INTEGER 00046 * The order of the matrix A. N >= 0. 00047 * 00048 * A (input/output) COMPLEX*16 array, dimension (LDA,N) 00049 * On entry, the symmetric matrix A. If UPLO = 'U', the leading 00050 * n by n upper triangular part of A contains the upper 00051 * triangular part of the matrix A, and the strictly lower 00052 * triangular part of A is not referenced. If UPLO = 'L', the 00053 * leading n by n lower triangular part of A contains the lower 00054 * triangular part of the matrix A, and the strictly upper 00055 * triangular part of A is not referenced. 00056 * 00057 * On exit, if INFO = 0, the factor U or L from the Cholesky 00058 * factorization as above. 00059 * 00060 * LDA (input) INTEGER 00061 * The leading dimension of the array A. LDA >= max(1,N). 00062 * 00063 * PIV (output) INTEGER array, dimension (N) 00064 * PIV is such that the nonzero entries are P( PIV(K), K ) = 1. 00065 * 00066 * RANK (output) INTEGER 00067 * The rank of A given by the number of steps the algorithm 00068 * completed. 00069 * 00070 * TOL (input) DOUBLE PRECISION 00071 * User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) ) 00072 * will be used. The algorithm terminates at the (K-1)st step 00073 * if the pivot <= TOL. 00074 * 00075 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N) 00076 * Work space. 00077 * 00078 * INFO (output) INTEGER 00079 * < 0: If INFO = -K, the K-th argument had an illegal value, 00080 * = 0: algorithm completed successfully, and 00081 * > 0: the matrix A is either rank deficient with computed rank 00082 * as returned in RANK, or is indefinite. See Section 7 of 00083 * LAPACK Working Note #161 for further information. 00084 * 00085 * ===================================================================== 00086 * 00087 * .. Parameters .. 00088 DOUBLE PRECISION ONE, ZERO 00089 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00090 COMPLEX*16 CONE 00091 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 00092 * .. 00093 * .. Local Scalars .. 00094 COMPLEX*16 ZTEMP 00095 DOUBLE PRECISION AJJ, DSTOP, DTEMP 00096 INTEGER I, ITEMP, J, JB, K, NB, PVT 00097 LOGICAL UPPER 00098 * .. 00099 * .. External Functions .. 00100 DOUBLE PRECISION DLAMCH 00101 INTEGER ILAENV 00102 LOGICAL LSAME, DISNAN 00103 EXTERNAL DLAMCH, ILAENV, LSAME, DISNAN 00104 * .. 00105 * .. External Subroutines .. 00106 EXTERNAL ZDSCAL, ZGEMV, ZHERK, ZLACGV, ZPSTF2, ZSWAP, 00107 $ XERBLA 00108 * .. 00109 * .. Intrinsic Functions .. 00110 INTRINSIC DBLE, DCONJG, MAX, MIN, SQRT, MAXLOC 00111 * .. 00112 * .. Executable Statements .. 00113 * 00114 * Test the input parameters. 00115 * 00116 INFO = 0 00117 UPPER = LSAME( UPLO, 'U' ) 00118 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00119 INFO = -1 00120 ELSE IF( N.LT.0 ) THEN 00121 INFO = -2 00122 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00123 INFO = -4 00124 END IF 00125 IF( INFO.NE.0 ) THEN 00126 CALL XERBLA( 'ZPSTRF', -INFO ) 00127 RETURN 00128 END IF 00129 * 00130 * Quick return if possible 00131 * 00132 IF( N.EQ.0 ) 00133 $ RETURN 00134 * 00135 * Get block size 00136 * 00137 NB = ILAENV( 1, 'ZPOTRF', UPLO, N, -1, -1, -1 ) 00138 IF( NB.LE.1 .OR. NB.GE.N ) THEN 00139 * 00140 * Use unblocked code 00141 * 00142 CALL ZPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK, 00143 $ INFO ) 00144 GO TO 230 00145 * 00146 ELSE 00147 * 00148 * Initialize PIV 00149 * 00150 DO 100 I = 1, N 00151 PIV( I ) = I 00152 100 CONTINUE 00153 * 00154 * Compute stopping value 00155 * 00156 DO 110 I = 1, N 00157 WORK( I ) = DBLE( A( I, I ) ) 00158 110 CONTINUE 00159 PVT = MAXLOC( WORK( 1:N ), 1 ) 00160 AJJ = DBLE( A( PVT, PVT ) ) 00161 IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN 00162 RANK = 0 00163 INFO = 1 00164 GO TO 230 00165 END IF 00166 * 00167 * Compute stopping value if not supplied 00168 * 00169 IF( TOL.LT.ZERO ) THEN 00170 DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ 00171 ELSE 00172 DSTOP = TOL 00173 END IF 00174 * 00175 * 00176 IF( UPPER ) THEN 00177 * 00178 * Compute the Cholesky factorization P**T * A * P = U**H * U 00179 * 00180 DO 160 K = 1, N, NB 00181 * 00182 * Account for last block not being NB wide 00183 * 00184 JB = MIN( NB, N-K+1 ) 00185 * 00186 * Set relevant part of first half of WORK to zero, 00187 * holds dot products 00188 * 00189 DO 120 I = K, N 00190 WORK( I ) = 0 00191 120 CONTINUE 00192 * 00193 DO 150 J = K, K + JB - 1 00194 * 00195 * Find pivot, test for exit, else swap rows and columns 00196 * Update dot products, compute possible pivots which are 00197 * stored in the second half of WORK 00198 * 00199 DO 130 I = J, N 00200 * 00201 IF( J.GT.K ) THEN 00202 WORK( I ) = WORK( I ) + 00203 $ DBLE( DCONJG( A( J-1, I ) )* 00204 $ A( J-1, I ) ) 00205 END IF 00206 WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I ) 00207 * 00208 130 CONTINUE 00209 * 00210 IF( J.GT.1 ) THEN 00211 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 00212 PVT = ITEMP + J - 1 00213 AJJ = WORK( N+PVT ) 00214 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN 00215 A( J, J ) = AJJ 00216 GO TO 220 00217 END IF 00218 END IF 00219 * 00220 IF( J.NE.PVT ) THEN 00221 * 00222 * Pivot OK, so can now swap pivot rows and columns 00223 * 00224 A( PVT, PVT ) = A( J, J ) 00225 CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 ) 00226 IF( PVT.LT.N ) 00227 $ CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA, 00228 $ A( PVT, PVT+1 ), LDA ) 00229 DO 140 I = J + 1, PVT - 1 00230 ZTEMP = DCONJG( A( J, I ) ) 00231 A( J, I ) = DCONJG( A( I, PVT ) ) 00232 A( I, PVT ) = ZTEMP 00233 140 CONTINUE 00234 A( J, PVT ) = DCONJG( A( J, PVT ) ) 00235 * 00236 * Swap dot products and PIV 00237 * 00238 DTEMP = WORK( J ) 00239 WORK( J ) = WORK( PVT ) 00240 WORK( PVT ) = DTEMP 00241 ITEMP = PIV( PVT ) 00242 PIV( PVT ) = PIV( J ) 00243 PIV( J ) = ITEMP 00244 END IF 00245 * 00246 AJJ = SQRT( AJJ ) 00247 A( J, J ) = AJJ 00248 * 00249 * Compute elements J+1:N of row J. 00250 * 00251 IF( J.LT.N ) THEN 00252 CALL ZLACGV( J-1, A( 1, J ), 1 ) 00253 CALL ZGEMV( 'Trans', J-K, N-J, -CONE, A( K, J+1 ), 00254 $ LDA, A( K, J ), 1, CONE, A( J, J+1 ), 00255 $ LDA ) 00256 CALL ZLACGV( J-1, A( 1, J ), 1 ) 00257 CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 00258 END IF 00259 * 00260 150 CONTINUE 00261 * 00262 * Update trailing matrix, J already incremented 00263 * 00264 IF( K+JB.LE.N ) THEN 00265 CALL ZHERK( 'Upper', 'Conj Trans', N-J+1, JB, -ONE, 00266 $ A( K, J ), LDA, ONE, A( J, J ), LDA ) 00267 END IF 00268 * 00269 160 CONTINUE 00270 * 00271 ELSE 00272 * 00273 * Compute the Cholesky factorization P**T * A * P = L * L**H 00274 * 00275 DO 210 K = 1, N, NB 00276 * 00277 * Account for last block not being NB wide 00278 * 00279 JB = MIN( NB, N-K+1 ) 00280 * 00281 * Set relevant part of first half of WORK to zero, 00282 * holds dot products 00283 * 00284 DO 170 I = K, N 00285 WORK( I ) = 0 00286 170 CONTINUE 00287 * 00288 DO 200 J = K, K + JB - 1 00289 * 00290 * Find pivot, test for exit, else swap rows and columns 00291 * Update dot products, compute possible pivots which are 00292 * stored in the second half of WORK 00293 * 00294 DO 180 I = J, N 00295 * 00296 IF( J.GT.K ) THEN 00297 WORK( I ) = WORK( I ) + 00298 $ DBLE( DCONJG( A( I, J-1 ) )* 00299 $ A( I, J-1 ) ) 00300 END IF 00301 WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I ) 00302 * 00303 180 CONTINUE 00304 * 00305 IF( J.GT.1 ) THEN 00306 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 00307 PVT = ITEMP + J - 1 00308 AJJ = WORK( N+PVT ) 00309 IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN 00310 A( J, J ) = AJJ 00311 GO TO 220 00312 END IF 00313 END IF 00314 * 00315 IF( J.NE.PVT ) THEN 00316 * 00317 * Pivot OK, so can now swap pivot rows and columns 00318 * 00319 A( PVT, PVT ) = A( J, J ) 00320 CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA ) 00321 IF( PVT.LT.N ) 00322 $ CALL ZSWAP( N-PVT, A( PVT+1, J ), 1, 00323 $ A( PVT+1, PVT ), 1 ) 00324 DO 190 I = J + 1, PVT - 1 00325 ZTEMP = DCONJG( A( I, J ) ) 00326 A( I, J ) = DCONJG( A( PVT, I ) ) 00327 A( PVT, I ) = ZTEMP 00328 190 CONTINUE 00329 A( PVT, J ) = DCONJG( A( PVT, J ) ) 00330 * 00331 * 00332 * Swap dot products and PIV 00333 * 00334 DTEMP = WORK( J ) 00335 WORK( J ) = WORK( PVT ) 00336 WORK( PVT ) = DTEMP 00337 ITEMP = PIV( PVT ) 00338 PIV( PVT ) = PIV( J ) 00339 PIV( J ) = ITEMP 00340 END IF 00341 * 00342 AJJ = SQRT( AJJ ) 00343 A( J, J ) = AJJ 00344 * 00345 * Compute elements J+1:N of column J. 00346 * 00347 IF( J.LT.N ) THEN 00348 CALL ZLACGV( J-1, A( J, 1 ), LDA ) 00349 CALL ZGEMV( 'No Trans', N-J, J-K, -CONE, 00350 $ A( J+1, K ), LDA, A( J, K ), LDA, CONE, 00351 $ A( J+1, J ), 1 ) 00352 CALL ZLACGV( J-1, A( J, 1 ), LDA ) 00353 CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 00354 END IF 00355 * 00356 200 CONTINUE 00357 * 00358 * Update trailing matrix, J already incremented 00359 * 00360 IF( K+JB.LE.N ) THEN 00361 CALL ZHERK( 'Lower', 'No Trans', N-J+1, JB, -ONE, 00362 $ A( J, K ), LDA, ONE, A( J, J ), LDA ) 00363 END IF 00364 * 00365 210 CONTINUE 00366 * 00367 END IF 00368 END IF 00369 * 00370 * Ran to completion, A has full rank 00371 * 00372 RANK = N 00373 * 00374 GO TO 230 00375 220 CONTINUE 00376 * 00377 * Rank is the number of steps completed. Set INFO = 1 to signal 00378 * that the factorization cannot be used to solve a system. 00379 * 00380 RANK = J - 1 00381 INFO = 1 00382 * 00383 230 CONTINUE 00384 RETURN 00385 * 00386 * End of ZPSTRF 00387 * 00388 END
1.7.3 
