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