Annotation of rpl/lapack/lapack/dpstf2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DPSTF2( 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: DOUBLE PRECISION A( LDA, * ), WORK( 2*N )
! 14: INTEGER PIV( N )
! 15: * ..
! 16: *
! 17: * Purpose
! 18: * =======
! 19: *
! 20: * DPSTF2 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 2 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: * PIV (output) INTEGER array, dimension (N)
! 57: * PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
! 58: *
! 59: * RANK (output) INTEGER
! 60: * The rank of A given by the number of steps the algorithm
! 61: * completed.
! 62: *
! 63: * TOL (input) DOUBLE PRECISION
! 64: * User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
! 65: * will be used. The algorithm terminates at the (K-1)st step
! 66: * if the pivot <= TOL.
! 67: *
! 68: * LDA (input) INTEGER
! 69: * The leading dimension of the array A. LDA >= max(1,N).
! 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, PVT
! 90: LOGICAL UPPER
! 91: * ..
! 92: * .. External Functions ..
! 93: DOUBLE PRECISION DLAMCH
! 94: LOGICAL LSAME, DISNAN
! 95: EXTERNAL DLAMCH, LSAME, DISNAN
! 96: * ..
! 97: * .. External Subroutines ..
! 98: EXTERNAL DGEMV, DSCAL, DSWAP, XERBLA
! 99: * ..
! 100: * .. Intrinsic Functions ..
! 101: INTRINSIC MAX, SQRT, MAXLOC
! 102: * ..
! 103: * .. Executable Statements ..
! 104: *
! 105: * Test the input parameters
! 106: *
! 107: INFO = 0
! 108: UPPER = LSAME( UPLO, 'U' )
! 109: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 110: INFO = -1
! 111: ELSE IF( N.LT.0 ) THEN
! 112: INFO = -2
! 113: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 114: INFO = -4
! 115: END IF
! 116: IF( INFO.NE.0 ) THEN
! 117: CALL XERBLA( 'DPSTF2', -INFO )
! 118: RETURN
! 119: END IF
! 120: *
! 121: * Quick return if possible
! 122: *
! 123: IF( N.EQ.0 )
! 124: $ RETURN
! 125: *
! 126: * Initialize PIV
! 127: *
! 128: DO 100 I = 1, N
! 129: PIV( I ) = I
! 130: 100 CONTINUE
! 131: *
! 132: * Compute stopping value
! 133: *
! 134: PVT = 1
! 135: AJJ = A( PVT, PVT )
! 136: DO I = 2, N
! 137: IF( A( I, I ).GT.AJJ ) THEN
! 138: PVT = I
! 139: AJJ = A( PVT, PVT )
! 140: END IF
! 141: END DO
! 142: IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
! 143: RANK = 0
! 144: INFO = 1
! 145: GO TO 170
! 146: END IF
! 147: *
! 148: * Compute stopping value if not supplied
! 149: *
! 150: IF( TOL.LT.ZERO ) THEN
! 151: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
! 152: ELSE
! 153: DSTOP = TOL
! 154: END IF
! 155: *
! 156: * Set first half of WORK to zero, holds dot products
! 157: *
! 158: DO 110 I = 1, N
! 159: WORK( I ) = 0
! 160: 110 CONTINUE
! 161: *
! 162: IF( UPPER ) THEN
! 163: *
! 164: * Compute the Cholesky factorization P' * A * P = U' * U
! 165: *
! 166: DO 130 J = 1, N
! 167: *
! 168: * Find pivot, test for exit, else swap rows and columns
! 169: * Update dot products, compute possible pivots which are
! 170: * stored in the second half of WORK
! 171: *
! 172: DO 120 I = J, N
! 173: *
! 174: IF( J.GT.1 ) THEN
! 175: WORK( I ) = WORK( I ) + A( J-1, I )**2
! 176: END IF
! 177: WORK( N+I ) = A( I, I ) - WORK( I )
! 178: *
! 179: 120 CONTINUE
! 180: *
! 181: IF( J.GT.1 ) THEN
! 182: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
! 183: PVT = ITEMP + J - 1
! 184: AJJ = WORK( N+PVT )
! 185: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
! 186: A( J, J ) = AJJ
! 187: GO TO 160
! 188: END IF
! 189: END IF
! 190: *
! 191: IF( J.NE.PVT ) THEN
! 192: *
! 193: * Pivot OK, so can now swap pivot rows and columns
! 194: *
! 195: A( PVT, PVT ) = A( J, J )
! 196: CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
! 197: IF( PVT.LT.N )
! 198: $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA,
! 199: $ A( PVT, PVT+1 ), LDA )
! 200: CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA, A( J+1, PVT ), 1 )
! 201: *
! 202: * Swap dot products and PIV
! 203: *
! 204: DTEMP = WORK( J )
! 205: WORK( J ) = WORK( PVT )
! 206: WORK( PVT ) = DTEMP
! 207: ITEMP = PIV( PVT )
! 208: PIV( PVT ) = PIV( J )
! 209: PIV( J ) = ITEMP
! 210: END IF
! 211: *
! 212: AJJ = SQRT( AJJ )
! 213: A( J, J ) = AJJ
! 214: *
! 215: * Compute elements J+1:N of row J
! 216: *
! 217: IF( J.LT.N ) THEN
! 218: CALL DGEMV( 'Trans', J-1, N-J, -ONE, A( 1, J+1 ), LDA,
! 219: $ A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
! 220: CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
! 221: END IF
! 222: *
! 223: 130 CONTINUE
! 224: *
! 225: ELSE
! 226: *
! 227: * Compute the Cholesky factorization P' * A * P = L * L'
! 228: *
! 229: DO 150 J = 1, N
! 230: *
! 231: * Find pivot, test for exit, else swap rows and columns
! 232: * Update dot products, compute possible pivots which are
! 233: * stored in the second half of WORK
! 234: *
! 235: DO 140 I = J, N
! 236: *
! 237: IF( J.GT.1 ) THEN
! 238: WORK( I ) = WORK( I ) + A( I, J-1 )**2
! 239: END IF
! 240: WORK( N+I ) = A( I, I ) - WORK( I )
! 241: *
! 242: 140 CONTINUE
! 243: *
! 244: IF( J.GT.1 ) THEN
! 245: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
! 246: PVT = ITEMP + J - 1
! 247: AJJ = WORK( N+PVT )
! 248: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
! 249: A( J, J ) = AJJ
! 250: GO TO 160
! 251: END IF
! 252: END IF
! 253: *
! 254: IF( J.NE.PVT ) THEN
! 255: *
! 256: * Pivot OK, so can now swap pivot rows and columns
! 257: *
! 258: A( PVT, PVT ) = A( J, J )
! 259: CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
! 260: IF( PVT.LT.N )
! 261: $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
! 262: $ 1 )
! 263: CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ), LDA )
! 264: *
! 265: * Swap dot products and PIV
! 266: *
! 267: DTEMP = WORK( J )
! 268: WORK( J ) = WORK( PVT )
! 269: WORK( PVT ) = DTEMP
! 270: ITEMP = PIV( PVT )
! 271: PIV( PVT ) = PIV( J )
! 272: PIV( J ) = ITEMP
! 273: END IF
! 274: *
! 275: AJJ = SQRT( AJJ )
! 276: A( J, J ) = AJJ
! 277: *
! 278: * Compute elements J+1:N of column J
! 279: *
! 280: IF( J.LT.N ) THEN
! 281: CALL DGEMV( 'No Trans', N-J, J-1, -ONE, A( J+1, 1 ), LDA,
! 282: $ A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
! 283: CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
! 284: END IF
! 285: *
! 286: 150 CONTINUE
! 287: *
! 288: END IF
! 289: *
! 290: * Ran to completion, A has full rank
! 291: *
! 292: RANK = N
! 293: *
! 294: GO TO 170
! 295: 160 CONTINUE
! 296: *
! 297: * Rank is number of steps completed. Set INFO = 1 to signal
! 298: * that the factorization cannot be used to solve a system.
! 299: *
! 300: RANK = J - 1
! 301: INFO = 1
! 302: *
! 303: 170 CONTINUE
! 304: RETURN
! 305: *
! 306: * End of DPSTF2
! 307: *
! 308: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dpstf2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack