Annotation of rpl/lapack/lapack/zcgesv.f, revision 1.8
1.1 bertrand 1: SUBROUTINE ZCGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
2: + SWORK, RWORK, ITER, INFO )
3: *
1.5 bertrand 4: * -- LAPACK PROTOTYPE driver routine (version 3.2.2) --
1.1 bertrand 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * January 2007
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: *
1.5 bertrand 66: * A (input/output) COMPLEX*16 array,
1.1 bertrand 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
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