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