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