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