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.3.1) --
5 *
6 * -- April 2011 --
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*1
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
2 $ SWORK, RWORK, ITER, INFO )
3 *
4 * -- LAPACK PROTOTYPE driver routine (version 3.3.1) --
5 *
6 * -- April 2011 --
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*1
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