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