Annotation of rpl/lapack/lapack/zpstf2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
! 2: *
! 3: * -- LAPACK PROTOTYPE routine (version 3.2.2) --
! 4: * Craig Lucas, University of Manchester / NAG Ltd.
! 5: * October, 2008
! 6: *
! 7: * .. Scalar Arguments ..
! 8: DOUBLE PRECISION TOL
! 9: INTEGER INFO, LDA, N, RANK
! 10: CHARACTER UPLO
! 11: * ..
! 12: * .. Array Arguments ..
! 13: COMPLEX*16 A( LDA, * )
! 14: DOUBLE PRECISION WORK( 2*N )
! 15: INTEGER PIV( N )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * ZPSTF2 computes the Cholesky factorization with complete
! 22: * pivoting of a complex Hermitian positive semidefinite matrix A.
! 23: *
! 24: * The factorization has the form
! 25: * P' * A * P = U' * U , if UPLO = 'U',
! 26: * P' * A * P = L * L', if UPLO = 'L',
! 27: * where U is an upper triangular matrix and L is lower triangular, and
! 28: * P is stored as vector PIV.
! 29: *
! 30: * This algorithm does not attempt to check that A is positive
! 31: * semidefinite. This version of the algorithm calls level 2 BLAS.
! 32: *
! 33: * Arguments
! 34: * =========
! 35: *
! 36: * UPLO (input) CHARACTER*1
! 37: * Specifies whether the upper or lower triangular part of the
! 38: * symmetric matrix A is stored.
! 39: * = 'U': Upper triangular
! 40: * = 'L': Lower triangular
! 41: *
! 42: * N (input) INTEGER
! 43: * The order of the matrix A. N >= 0.
! 44: *
! 45: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
! 46: * On entry, the symmetric matrix A. If UPLO = 'U', the leading
! 47: * n by n upper triangular part of A contains the upper
! 48: * triangular part of the matrix A, and the strictly lower
! 49: * triangular part of A is not referenced. If UPLO = 'L', the
! 50: * leading n by n lower triangular part of A contains the lower
! 51: * triangular part of the matrix A, and the strictly upper
! 52: * triangular part of A is not referenced.
! 53: *
! 54: * On exit, if INFO = 0, the factor U or L from the Cholesky
! 55: * factorization as above.
! 56: *
! 57: * PIV (output) INTEGER array, dimension (N)
! 58: * PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
! 59: *
! 60: * RANK (output) INTEGER
! 61: * The rank of A given by the number of steps the algorithm
! 62: * completed.
! 63: *
! 64: * TOL (input) DOUBLE PRECISION
! 65: * User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
! 66: * will be used. The algorithm terminates at the (K-1)st step
! 67: * if the pivot <= TOL.
! 68: *
! 69: * LDA (input) INTEGER
! 70: * The leading dimension of the array A. LDA >= max(1,N).
! 71: *
! 72: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
! 73: * Work space.
! 74: *
! 75: * INFO (output) INTEGER
! 76: * < 0: If INFO = -K, the K-th argument had an illegal value,
! 77: * = 0: algorithm completed successfully, and
! 78: * > 0: the matrix A is either rank deficient with computed rank
! 79: * as returned in RANK, or is indefinite. See Section 7 of
! 80: * LAPACK Working Note #161 for further information.
! 81: *
! 82: * =====================================================================
! 83: *
! 84: * .. Parameters ..
! 85: DOUBLE PRECISION ONE, ZERO
! 86: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 87: COMPLEX*16 CONE
! 88: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
! 89: * ..
! 90: * .. Local Scalars ..
! 91: COMPLEX*16 ZTEMP
! 92: DOUBLE PRECISION AJJ, DSTOP, DTEMP
! 93: INTEGER I, ITEMP, J, PVT
! 94: LOGICAL UPPER
! 95: * ..
! 96: * .. External Functions ..
! 97: DOUBLE PRECISION DLAMCH
! 98: LOGICAL LSAME, DISNAN
! 99: EXTERNAL DLAMCH, LSAME, DISNAN
! 100: * ..
! 101: * .. External Subroutines ..
! 102: EXTERNAL ZDSCAL, ZGEMV, ZLACGV, ZSWAP, XERBLA
! 103: * ..
! 104: * .. Intrinsic Functions ..
! 105: INTRINSIC DBLE, DCONJG, MAX, SQRT
! 106: * ..
! 107: * .. Executable Statements ..
! 108: *
! 109: * Test the input parameters
! 110: *
! 111: INFO = 0
! 112: UPPER = LSAME( UPLO, 'U' )
! 113: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 114: INFO = -1
! 115: ELSE IF( N.LT.0 ) THEN
! 116: INFO = -2
! 117: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 118: INFO = -4
! 119: END IF
! 120: IF( INFO.NE.0 ) THEN
! 121: CALL XERBLA( 'ZPSTF2', -INFO )
! 122: RETURN
! 123: END IF
! 124: *
! 125: * Quick return if possible
! 126: *
! 127: IF( N.EQ.0 )
! 128: $ RETURN
! 129: *
! 130: * Initialize PIV
! 131: *
! 132: DO 100 I = 1, N
! 133: PIV( I ) = I
! 134: 100 CONTINUE
! 135: *
! 136: * Compute stopping value
! 137: *
! 138: DO 110 I = 1, N
! 139: WORK( I ) = DBLE( A( I, I ) )
! 140: 110 CONTINUE
! 141: PVT = MAXLOC( WORK( 1:N ), 1 )
! 142: AJJ = DBLE( A( PVT, PVT ) )
! 143: IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
! 144: RANK = 0
! 145: INFO = 1
! 146: GO TO 200
! 147: END IF
! 148: *
! 149: * Compute stopping value if not supplied
! 150: *
! 151: IF( TOL.LT.ZERO ) THEN
! 152: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
! 153: ELSE
! 154: DSTOP = TOL
! 155: END IF
! 156: *
! 157: * Set first half of WORK to zero, holds dot products
! 158: *
! 159: DO 120 I = 1, N
! 160: WORK( I ) = 0
! 161: 120 CONTINUE
! 162: *
! 163: IF( UPPER ) THEN
! 164: *
! 165: * Compute the Cholesky factorization P' * A * P = U' * U
! 166: *
! 167: DO 150 J = 1, N
! 168: *
! 169: * Find pivot, test for exit, else swap rows and columns
! 170: * Update dot products, compute possible pivots which are
! 171: * stored in the second half of WORK
! 172: *
! 173: DO 130 I = J, N
! 174: *
! 175: IF( J.GT.1 ) THEN
! 176: WORK( I ) = WORK( I ) +
! 177: $ DBLE( DCONJG( A( J-1, I ) )*
! 178: $ A( J-1, I ) )
! 179: END IF
! 180: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
! 181: *
! 182: 130 CONTINUE
! 183: *
! 184: IF( J.GT.1 ) THEN
! 185: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
! 186: PVT = ITEMP + J - 1
! 187: AJJ = WORK( N+PVT )
! 188: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
! 189: A( J, J ) = AJJ
! 190: GO TO 190
! 191: END IF
! 192: END IF
! 193: *
! 194: IF( J.NE.PVT ) THEN
! 195: *
! 196: * Pivot OK, so can now swap pivot rows and columns
! 197: *
! 198: A( PVT, PVT ) = A( J, J )
! 199: CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
! 200: IF( PVT.LT.N )
! 201: $ CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA,
! 202: $ A( PVT, PVT+1 ), LDA )
! 203: DO 140 I = J + 1, PVT - 1
! 204: ZTEMP = DCONJG( A( J, I ) )
! 205: A( J, I ) = DCONJG( A( I, PVT ) )
! 206: A( I, PVT ) = ZTEMP
! 207: 140 CONTINUE
! 208: A( J, PVT ) = DCONJG( A( J, PVT ) )
! 209: *
! 210: * Swap dot products and PIV
! 211: *
! 212: DTEMP = WORK( J )
! 213: WORK( J ) = WORK( PVT )
! 214: WORK( PVT ) = DTEMP
! 215: ITEMP = PIV( PVT )
! 216: PIV( PVT ) = PIV( J )
! 217: PIV( J ) = ITEMP
! 218: END IF
! 219: *
! 220: AJJ = SQRT( AJJ )
! 221: A( J, J ) = AJJ
! 222: *
! 223: * Compute elements J+1:N of row J
! 224: *
! 225: IF( J.LT.N ) THEN
! 226: CALL ZLACGV( J-1, A( 1, J ), 1 )
! 227: CALL ZGEMV( 'Trans', J-1, N-J, -CONE, A( 1, J+1 ), LDA,
! 228: $ A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
! 229: CALL ZLACGV( J-1, A( 1, J ), 1 )
! 230: CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
! 231: END IF
! 232: *
! 233: 150 CONTINUE
! 234: *
! 235: ELSE
! 236: *
! 237: * Compute the Cholesky factorization P' * A * P = L * L'
! 238: *
! 239: DO 180 J = 1, N
! 240: *
! 241: * Find pivot, test for exit, else swap rows and columns
! 242: * Update dot products, compute possible pivots which are
! 243: * stored in the second half of WORK
! 244: *
! 245: DO 160 I = J, N
! 246: *
! 247: IF( J.GT.1 ) THEN
! 248: WORK( I ) = WORK( I ) +
! 249: $ DBLE( DCONJG( A( I, J-1 ) )*
! 250: $ A( I, J-1 ) )
! 251: END IF
! 252: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
! 253: *
! 254: 160 CONTINUE
! 255: *
! 256: IF( J.GT.1 ) THEN
! 257: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
! 258: PVT = ITEMP + J - 1
! 259: AJJ = WORK( N+PVT )
! 260: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
! 261: A( J, J ) = AJJ
! 262: GO TO 190
! 263: END IF
! 264: END IF
! 265: *
! 266: IF( J.NE.PVT ) THEN
! 267: *
! 268: * Pivot OK, so can now swap pivot rows and columns
! 269: *
! 270: A( PVT, PVT ) = A( J, J )
! 271: CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
! 272: IF( PVT.LT.N )
! 273: $ CALL ZSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
! 274: $ 1 )
! 275: DO 170 I = J + 1, PVT - 1
! 276: ZTEMP = DCONJG( A( I, J ) )
! 277: A( I, J ) = DCONJG( A( PVT, I ) )
! 278: A( PVT, I ) = ZTEMP
! 279: 170 CONTINUE
! 280: A( PVT, J ) = DCONJG( A( PVT, J ) )
! 281: *
! 282: * Swap dot products and PIV
! 283: *
! 284: DTEMP = WORK( J )
! 285: WORK( J ) = WORK( PVT )
! 286: WORK( PVT ) = DTEMP
! 287: ITEMP = PIV( PVT )
! 288: PIV( PVT ) = PIV( J )
! 289: PIV( J ) = ITEMP
! 290: END IF
! 291: *
! 292: AJJ = SQRT( AJJ )
! 293: A( J, J ) = AJJ
! 294: *
! 295: * Compute elements J+1:N of column J
! 296: *
! 297: IF( J.LT.N ) THEN
! 298: CALL ZLACGV( J-1, A( J, 1 ), LDA )
! 299: CALL ZGEMV( 'No Trans', N-J, J-1, -CONE, A( J+1, 1 ),
! 300: $ LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
! 301: CALL ZLACGV( J-1, A( J, 1 ), LDA )
! 302: CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
! 303: END IF
! 304: *
! 305: 180 CONTINUE
! 306: *
! 307: END IF
! 308: *
! 309: * Ran to completion, A has full rank
! 310: *
! 311: RANK = N
! 312: *
! 313: GO TO 200
! 314: 190 CONTINUE
! 315: *
! 316: * Rank is number of steps completed. Set INFO = 1 to signal
! 317: * that the factorization cannot be used to solve a system.
! 318: *
! 319: RANK = J - 1
! 320: INFO = 1
! 321: *
! 322: 200 CONTINUE
! 323: RETURN
! 324: *
! 325: * End of ZPSTF2
! 326: *
! 327: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zpstf2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack