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