Annotation of rpl/lapack/lapack/zcposv.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
! 2: + SWORK, RWORK, ITER, INFO )
! 3: *
! 4: * -- LAPACK PROTOTYPE driver routine (version 3.2.2) --
! 5: *
! 6: * -- June 2010 --
! 7: *
! 8: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 10: * ..
! 11: * .. Scalar Arguments ..
! 12: CHARACTER UPLO
! 13: INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
! 14: * ..
! 15: * .. Array Arguments ..
! 16: DOUBLE PRECISION RWORK( * )
! 17: COMPLEX SWORK( * )
! 18: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
! 19: + X( LDX, * )
! 20: * ..
! 21: *
! 22: * Purpose
! 23: * =======
! 24: *
! 25: * ZCPOSV computes the solution to a complex system of linear equations
! 26: * A * X = B,
! 27: * where A is an N-by-N Hermitian positive definite matrix and X and B
! 28: * are N-by-NRHS matrices.
! 29: *
! 30: * ZCPOSV first attempts to factorize the matrix in COMPLEX and use this
! 31: * factorization within an iterative refinement procedure to produce a
! 32: * solution with COMPLEX*16 normwise backward error quality (see below).
! 33: * If the approach fails the method switches to a COMPLEX*16
! 34: * factorization and solve.
! 35: *
! 36: * The iterative refinement is not going to be a winning strategy if
! 37: * the ratio COMPLEX performance over COMPLEX*16 performance is too
! 38: * small. A reasonable strategy should take the number of right-hand
! 39: * sides and the size of the matrix into account. This might be done
! 40: * with a call to ILAENV in the future. Up to now, we always try
! 41: * iterative refinement.
! 42: *
! 43: * The iterative refinement process is stopped if
! 44: * ITER > ITERMAX
! 45: * or for all the RHS we have:
! 46: * RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
! 47: * where
! 48: * o ITER is the number of the current iteration in the iterative
! 49: * refinement process
! 50: * o RNRM is the infinity-norm of the residual
! 51: * o XNRM is the infinity-norm of the solution
! 52: * o ANRM is the infinity-operator-norm of the matrix A
! 53: * o EPS is the machine epsilon returned by DLAMCH('Epsilon')
! 54: * The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
! 55: * respectively.
! 56: *
! 57: * Arguments
! 58: * =========
! 59: *
! 60: * UPLO (input) CHARACTER
! 61: * = 'U': Upper triangle of A is stored;
! 62: * = 'L': Lower triangle of A is stored.
! 63: *
! 64: * N (input) INTEGER
! 65: * The number of linear equations, i.e., the order of the
! 66: * matrix A. N >= 0.
! 67: *
! 68: * NRHS (input) INTEGER
! 69: * The number of right hand sides, i.e., the number of columns
! 70: * of the matrix B. NRHS >= 0.
! 71: *
! 72: * A (input/output) COMPLEX*16 array,
! 73: * dimension (LDA,N)
! 74: * On entry, the Hermitian matrix A. If UPLO = 'U', the leading
! 75: * N-by-N upper triangular part of A contains the upper
! 76: * triangular part of the matrix A, and the strictly lower
! 77: * triangular part of A is not referenced. If UPLO = 'L', the
! 78: * leading N-by-N lower triangular part of A contains the lower
! 79: * triangular part of the matrix A, and the strictly upper
! 80: * triangular part of A is not referenced.
! 81: *
! 82: * Note that the imaginary parts of the diagonal
! 83: * elements need not be set and are assumed to be zero.
! 84: *
! 85: * On exit, if iterative refinement has been successfully used
! 86: * (INFO.EQ.0 and ITER.GE.0, see description below), then A is
! 87: * unchanged, if double precision factorization has been used
! 88: * (INFO.EQ.0 and ITER.LT.0, see description below), then the
! 89: * array A contains the factor U or L from the Cholesky
! 90: * factorization A = U**H*U or A = L*L**H.
! 91: *
! 92: * LDA (input) INTEGER
! 93: * The leading dimension of the array A. LDA >= max(1,N).
! 94: *
! 95: * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
! 96: * The N-by-NRHS right hand side matrix B.
! 97: *
! 98: * LDB (input) INTEGER
! 99: * The leading dimension of the array B. LDB >= max(1,N).
! 100: *
! 101: * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
! 102: * If INFO = 0, the N-by-NRHS solution matrix X.
! 103: *
! 104: * LDX (input) INTEGER
! 105: * The leading dimension of the array X. LDX >= max(1,N).
! 106: *
! 107: * WORK (workspace) COMPLEX*16 array, dimension (N*NRHS)
! 108: * This array is used to hold the residual vectors.
! 109: *
! 110: * SWORK (workspace) COMPLEX array, dimension (N*(N+NRHS))
! 111: * This array is used to use the single precision matrix and the
! 112: * right-hand sides or solutions in single precision.
! 113: *
! 114: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
! 115: *
! 116: * ITER (output) INTEGER
! 117: * < 0: iterative refinement has failed, COMPLEX*16
! 118: * factorization has been performed
! 119: * -1 : the routine fell back to full precision for
! 120: * implementation- or machine-specific reasons
! 121: * -2 : narrowing the precision induced an overflow,
! 122: * the routine fell back to full precision
! 123: * -3 : failure of CPOTRF
! 124: * -31: stop the iterative refinement after the 30th
! 125: * iterations
! 126: * > 0: iterative refinement has been sucessfully used.
! 127: * Returns the number of iterations
! 128: *
! 129: * INFO (output) INTEGER
! 130: * = 0: successful exit
! 131: * < 0: if INFO = -i, the i-th argument had an illegal value
! 132: * > 0: if INFO = i, the leading minor of order i of
! 133: * (COMPLEX*16) A is not positive definite, so the
! 134: * factorization could not be completed, and the solution
! 135: * has not been computed.
! 136: *
! 137: * =========
! 138: *
! 139: * .. Parameters ..
! 140: LOGICAL DOITREF
! 141: PARAMETER ( DOITREF = .TRUE. )
! 142: *
! 143: INTEGER ITERMAX
! 144: PARAMETER ( ITERMAX = 30 )
! 145: *
! 146: DOUBLE PRECISION BWDMAX
! 147: PARAMETER ( BWDMAX = 1.0E+00 )
! 148: *
! 149: COMPLEX*16 NEGONE, ONE
! 150: PARAMETER ( NEGONE = ( -1.0D+00, 0.0D+00 ),
! 151: + ONE = ( 1.0D+00, 0.0D+00 ) )
! 152: *
! 153: * .. Local Scalars ..
! 154: INTEGER I, IITER, PTSA, PTSX
! 155: DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
! 156: COMPLEX*16 ZDUM
! 157: *
! 158: * .. External Subroutines ..
! 159: EXTERNAL ZAXPY, ZHEMM, ZLACPY, ZLAT2C, ZLAG2C, CLAG2Z,
! 160: + CPOTRF, CPOTRS, XERBLA
! 161: * ..
! 162: * .. External Functions ..
! 163: INTEGER IZAMAX
! 164: DOUBLE PRECISION DLAMCH, ZLANHE
! 165: LOGICAL LSAME
! 166: EXTERNAL IZAMAX, DLAMCH, ZLANHE, LSAME
! 167: * ..
! 168: * .. Intrinsic Functions ..
! 169: INTRINSIC ABS, DBLE, MAX, SQRT
! 170: * .. Statement Functions ..
! 171: DOUBLE PRECISION CABS1
! 172: * ..
! 173: * .. Statement Function definitions ..
! 174: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
! 175: * ..
! 176: * .. Executable Statements ..
! 177: *
! 178: INFO = 0
! 179: ITER = 0
! 180: *
! 181: * Test the input parameters.
! 182: *
! 183: IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 184: INFO = -1
! 185: ELSE IF( N.LT.0 ) THEN
! 186: INFO = -2
! 187: ELSE IF( NRHS.LT.0 ) THEN
! 188: INFO = -3
! 189: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 190: INFO = -5
! 191: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 192: INFO = -7
! 193: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
! 194: INFO = -9
! 195: END IF
! 196: IF( INFO.NE.0 ) THEN
! 197: CALL XERBLA( 'ZCPOSV', -INFO )
! 198: RETURN
! 199: END IF
! 200: *
! 201: * Quick return if (N.EQ.0).
! 202: *
! 203: IF( N.EQ.0 )
! 204: + RETURN
! 205: *
! 206: * Skip single precision iterative refinement if a priori slower
! 207: * than double precision factorization.
! 208: *
! 209: IF( .NOT.DOITREF ) THEN
! 210: ITER = -1
! 211: GO TO 40
! 212: END IF
! 213: *
! 214: * Compute some constants.
! 215: *
! 216: ANRM = ZLANHE( 'I', UPLO, N, A, LDA, RWORK )
! 217: EPS = DLAMCH( 'Epsilon' )
! 218: CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
! 219: *
! 220: * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
! 221: *
! 222: PTSA = 1
! 223: PTSX = PTSA + N*N
! 224: *
! 225: * Convert B from double precision to single precision and store the
! 226: * result in SX.
! 227: *
! 228: CALL ZLAG2C( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
! 229: *
! 230: IF( INFO.NE.0 ) THEN
! 231: ITER = -2
! 232: GO TO 40
! 233: END IF
! 234: *
! 235: * Convert A from double precision to single precision and store the
! 236: * result in SA.
! 237: *
! 238: CALL ZLAT2C( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO )
! 239: *
! 240: IF( INFO.NE.0 ) THEN
! 241: ITER = -2
! 242: GO TO 40
! 243: END IF
! 244: *
! 245: * Compute the Cholesky factorization of SA.
! 246: *
! 247: CALL CPOTRF( UPLO, N, SWORK( PTSA ), N, INFO )
! 248: *
! 249: IF( INFO.NE.0 ) THEN
! 250: ITER = -3
! 251: GO TO 40
! 252: END IF
! 253: *
! 254: * Solve the system SA*SX = SB.
! 255: *
! 256: CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
! 257: + INFO )
! 258: *
! 259: * Convert SX back to COMPLEX*16
! 260: *
! 261: CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
! 262: *
! 263: * Compute R = B - AX (R is WORK).
! 264: *
! 265: CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
! 266: *
! 267: CALL ZHEMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
! 268: + WORK, N )
! 269: *
! 270: * Check whether the NRHS normwise backward errors satisfy the
! 271: * stopping criterion. If yes, set ITER=0 and return.
! 272: *
! 273: DO I = 1, NRHS
! 274: XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
! 275: RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) )
! 276: IF( RNRM.GT.XNRM*CTE )
! 277: + GO TO 10
! 278: END DO
! 279: *
! 280: * If we are here, the NRHS normwise backward errors satisfy the
! 281: * stopping criterion. We are good to exit.
! 282: *
! 283: ITER = 0
! 284: RETURN
! 285: *
! 286: 10 CONTINUE
! 287: *
! 288: DO 30 IITER = 1, ITERMAX
! 289: *
! 290: * Convert R (in WORK) from double precision to single precision
! 291: * and store the result in SX.
! 292: *
! 293: CALL ZLAG2C( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
! 294: *
! 295: IF( INFO.NE.0 ) THEN
! 296: ITER = -2
! 297: GO TO 40
! 298: END IF
! 299: *
! 300: * Solve the system SA*SX = SR.
! 301: *
! 302: CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
! 303: + INFO )
! 304: *
! 305: * Convert SX back to double precision and update the current
! 306: * iterate.
! 307: *
! 308: CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
! 309: *
! 310: DO I = 1, NRHS
! 311: CALL ZAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
! 312: END DO
! 313: *
! 314: * Compute R = B - AX (R is WORK).
! 315: *
! 316: CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
! 317: *
! 318: CALL ZHEMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
! 319: + WORK, N )
! 320: *
! 321: * Check whether the NRHS normwise backward errors satisfy the
! 322: * stopping criterion. If yes, set ITER=IITER>0 and return.
! 323: *
! 324: DO I = 1, NRHS
! 325: XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
! 326: RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) )
! 327: IF( RNRM.GT.XNRM*CTE )
! 328: + GO TO 20
! 329: END DO
! 330: *
! 331: * If we are here, the NRHS normwise backward errors satisfy the
! 332: * stopping criterion, we are good to exit.
! 333: *
! 334: ITER = IITER
! 335: *
! 336: RETURN
! 337: *
! 338: 20 CONTINUE
! 339: *
! 340: 30 CONTINUE
! 341: *
! 342: * If we are at this place of the code, this is because we have
! 343: * performed ITER=ITERMAX iterations and never satisified the
! 344: * stopping criterion, set up the ITER flag accordingly and follow
! 345: * up on double precision routine.
! 346: *
! 347: ITER = -ITERMAX - 1
! 348: *
! 349: 40 CONTINUE
! 350: *
! 351: * Single-precision iterative refinement failed to converge to a
! 352: * satisfactory solution, so we resort to double precision.
! 353: *
! 354: CALL ZPOTRF( UPLO, N, A, LDA, INFO )
! 355: *
! 356: IF( INFO.NE.0 )
! 357: + RETURN
! 358: *
! 359: CALL ZLACPY( 'All', N, NRHS, B, LDB, X, LDX )
! 360: CALL ZPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )
! 361: *
! 362: RETURN
! 363: *
! 364: * End of ZCPOSV.
! 365: *
! 366: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zcposv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack