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