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 = 0 and ITER >= 0, see description below), then A is
115: *> unchanged, if double precision factorization has been used
116: *> (INFO = 0 and ITER < 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 successfully 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: *> \ingroup complex16POsolve
205: *
206: * =====================================================================
207: SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
208: $ SWORK, RWORK, ITER, INFO )
209: *
210: * -- LAPACK driver routine --
211: * -- LAPACK is a software package provided by Univ. of Tennessee, --
212: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213: *
214: * .. Scalar Arguments ..
215: CHARACTER UPLO
216: INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
217: * ..
218: * .. Array Arguments ..
219: DOUBLE PRECISION RWORK( * )
220: COMPLEX SWORK( * )
221: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
222: $ X( LDX, * )
223: * ..
224: *
225: * =====================================================================
226: *
227: * .. Parameters ..
228: LOGICAL DOITREF
229: PARAMETER ( DOITREF = .TRUE. )
230: *
231: INTEGER ITERMAX
232: PARAMETER ( ITERMAX = 30 )
233: *
234: DOUBLE PRECISION BWDMAX
235: PARAMETER ( BWDMAX = 1.0E+00 )
236: *
237: COMPLEX*16 NEGONE, ONE
238: PARAMETER ( NEGONE = ( -1.0D+00, 0.0D+00 ),
239: $ ONE = ( 1.0D+00, 0.0D+00 ) )
240: *
241: * .. Local Scalars ..
242: INTEGER I, IITER, PTSA, PTSX
243: DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
244: COMPLEX*16 ZDUM
245: *
246: * .. External Subroutines ..
247: EXTERNAL ZAXPY, ZHEMM, ZLACPY, ZLAT2C, ZLAG2C, CLAG2Z,
248: $ CPOTRF, CPOTRS, XERBLA, ZPOTRF, ZPOTRS
249: * ..
250: * .. External Functions ..
251: INTEGER IZAMAX
252: DOUBLE PRECISION DLAMCH, ZLANHE
253: LOGICAL LSAME
254: EXTERNAL IZAMAX, DLAMCH, ZLANHE, LSAME
255: * ..
256: * .. Intrinsic Functions ..
257: INTRINSIC ABS, DBLE, MAX, SQRT
258: * .. Statement Functions ..
259: DOUBLE PRECISION CABS1
260: * ..
261: * .. Statement Function definitions ..
262: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
263: * ..
264: * .. Executable Statements ..
265: *
266: INFO = 0
267: ITER = 0
268: *
269: * Test the input parameters.
270: *
271: IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
272: INFO = -1
273: ELSE IF( N.LT.0 ) THEN
274: INFO = -2
275: ELSE IF( NRHS.LT.0 ) THEN
276: INFO = -3
277: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
278: INFO = -5
279: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
280: INFO = -7
281: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
282: INFO = -9
283: END IF
284: IF( INFO.NE.0 ) THEN
285: CALL XERBLA( 'ZCPOSV', -INFO )
286: RETURN
287: END IF
288: *
289: * Quick return if (N.EQ.0).
290: *
291: IF( N.EQ.0 )
292: $ RETURN
293: *
294: * Skip single precision iterative refinement if a priori slower
295: * than double precision factorization.
296: *
297: IF( .NOT.DOITREF ) THEN
298: ITER = -1
299: GO TO 40
300: END IF
301: *
302: * Compute some constants.
303: *
304: ANRM = ZLANHE( 'I', UPLO, N, A, LDA, RWORK )
305: EPS = DLAMCH( 'Epsilon' )
306: CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
307: *
308: * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
309: *
310: PTSA = 1
311: PTSX = PTSA + N*N
312: *
313: * Convert B from double precision to single precision and store the
314: * result in SX.
315: *
316: CALL ZLAG2C( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
317: *
318: IF( INFO.NE.0 ) THEN
319: ITER = -2
320: GO TO 40
321: END IF
322: *
323: * Convert A from double precision to single precision and store the
324: * result in SA.
325: *
326: CALL ZLAT2C( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO )
327: *
328: IF( INFO.NE.0 ) THEN
329: ITER = -2
330: GO TO 40
331: END IF
332: *
333: * Compute the Cholesky factorization of SA.
334: *
335: CALL CPOTRF( UPLO, N, SWORK( PTSA ), N, INFO )
336: *
337: IF( INFO.NE.0 ) THEN
338: ITER = -3
339: GO TO 40
340: END IF
341: *
342: * Solve the system SA*SX = SB.
343: *
344: CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
345: $ INFO )
346: *
347: * Convert SX back to COMPLEX*16
348: *
349: CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
350: *
351: * Compute R = B - AX (R is WORK).
352: *
353: CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
354: *
355: CALL ZHEMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
356: $ WORK, N )
357: *
358: * Check whether the NRHS normwise backward errors satisfy the
359: * stopping criterion. If yes, set ITER=0 and return.
360: *
361: DO I = 1, NRHS
362: XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
363: RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) )
364: IF( RNRM.GT.XNRM*CTE )
365: $ GO TO 10
366: END DO
367: *
368: * If we are here, the NRHS normwise backward errors satisfy the
369: * stopping criterion. We are good to exit.
370: *
371: ITER = 0
372: RETURN
373: *
374: 10 CONTINUE
375: *
376: DO 30 IITER = 1, ITERMAX
377: *
378: * Convert R (in WORK) from double precision to single precision
379: * and store the result in SX.
380: *
381: CALL ZLAG2C( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
382: *
383: IF( INFO.NE.0 ) THEN
384: ITER = -2
385: GO TO 40
386: END IF
387: *
388: * Solve the system SA*SX = SR.
389: *
390: CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
391: $ INFO )
392: *
393: * Convert SX back to double precision and update the current
394: * iterate.
395: *
396: CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
397: *
398: DO I = 1, NRHS
399: CALL ZAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
400: END DO
401: *
402: * Compute R = B - AX (R is WORK).
403: *
404: CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
405: *
406: CALL ZHEMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
407: $ WORK, N )
408: *
409: * Check whether the NRHS normwise backward errors satisfy the
410: * stopping criterion. If yes, set ITER=IITER>0 and return.
411: *
412: DO I = 1, NRHS
413: XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
414: RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) )
415: IF( RNRM.GT.XNRM*CTE )
416: $ GO TO 20
417: END DO
418: *
419: * If we are here, the NRHS normwise backward errors satisfy the
420: * stopping criterion, we are good to exit.
421: *
422: ITER = IITER
423: *
424: RETURN
425: *
426: 20 CONTINUE
427: *
428: 30 CONTINUE
429: *
430: * If we are at this place of the code, this is because we have
431: * performed ITER=ITERMAX iterations and never satisfied the
432: * stopping criterion, set up the ITER flag accordingly and follow
433: * up on double precision routine.
434: *
435: ITER = -ITERMAX - 1
436: *
437: 40 CONTINUE
438: *
439: * Single-precision iterative refinement failed to converge to a
440: * satisfactory solution, so we resort to double precision.
441: *
442: CALL ZPOTRF( UPLO, N, A, LDA, INFO )
443: *
444: IF( INFO.NE.0 )
445: $ RETURN
446: *
447: CALL ZLACPY( 'All', N, NRHS, B, LDB, X, LDX )
448: CALL ZPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )
449: *
450: RETURN
451: *
452: * End of ZCPOSV
453: *
454: 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