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