Annotation of rpl/lapack/lapack/zgesvx.f, revision 1.7
1.1 bertrand 1: SUBROUTINE ZGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
2: $ EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR,
3: $ WORK, RWORK, INFO )
4: *
5: * -- LAPACK driver routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: CHARACTER EQUED, FACT, TRANS
12: INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
13: DOUBLE PRECISION RCOND
14: * ..
15: * .. Array Arguments ..
16: INTEGER IPIV( * )
17: DOUBLE PRECISION BERR( * ), C( * ), FERR( * ), R( * ),
18: $ RWORK( * )
19: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
20: $ WORK( * ), X( LDX, * )
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * ZGESVX uses the LU factorization to compute the solution to a complex
27: * system of linear equations
28: * A * X = B,
29: * where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
30: *
31: * Error bounds on the solution and a condition estimate are also
32: * provided.
33: *
34: * Description
35: * ===========
36: *
37: * The following steps are performed:
38: *
39: * 1. If FACT = 'E', real scaling factors are computed to equilibrate
40: * the system:
41: * TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
42: * TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
43: * TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
44: * Whether or not the system will be equilibrated depends on the
45: * scaling of the matrix A, but if equilibration is used, A is
46: * overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
47: * or diag(C)*B (if TRANS = 'T' or 'C').
48: *
49: * 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
50: * matrix A (after equilibration if FACT = 'E') as
51: * A = P * L * U,
52: * where P is a permutation matrix, L is a unit lower triangular
53: * matrix, and U is upper triangular.
54: *
55: * 3. If some U(i,i)=0, so that U is exactly singular, then the routine
56: * returns with INFO = i. Otherwise, the factored form of A is used
57: * to estimate the condition number of the matrix A. If the
58: * reciprocal of the condition number is less than machine precision,
59: * INFO = N+1 is returned as a warning, but the routine still goes on
60: * to solve for X and compute error bounds as described below.
61: *
62: * 4. The system of equations is solved for X using the factored form
63: * of A.
64: *
65: * 5. Iterative refinement is applied to improve the computed solution
66: * matrix and calculate error bounds and backward error estimates
67: * for it.
68: *
69: * 6. If equilibration was used, the matrix X is premultiplied by
70: * diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
71: * that it solves the original system before equilibration.
72: *
73: * Arguments
74: * =========
75: *
76: * FACT (input) CHARACTER*1
77: * Specifies whether or not the factored form of the matrix A is
78: * supplied on entry, and if not, whether the matrix A should be
79: * equilibrated before it is factored.
80: * = 'F': On entry, AF and IPIV contain the factored form of A.
81: * If EQUED is not 'N', the matrix A has been
82: * equilibrated with scaling factors given by R and C.
83: * A, AF, and IPIV are not modified.
84: * = 'N': The matrix A will be copied to AF and factored.
85: * = 'E': The matrix A will be equilibrated if necessary, then
86: * copied to AF and factored.
87: *
88: * TRANS (input) CHARACTER*1
89: * Specifies the form of the system of equations:
90: * = 'N': A * X = B (No transpose)
91: * = 'T': A**T * X = B (Transpose)
92: * = 'C': A**H * X = B (Conjugate transpose)
93: *
94: * N (input) INTEGER
95: * The number of linear equations, i.e., the order of the
96: * matrix A. N >= 0.
97: *
98: * NRHS (input) INTEGER
99: * The number of right hand sides, i.e., the number of columns
100: * of the matrices B and X. NRHS >= 0.
101: *
102: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
103: * On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is
104: * not 'N', then A must have been equilibrated by the scaling
105: * factors in R and/or C. A is not modified if FACT = 'F' or
106: * 'N', or if FACT = 'E' and EQUED = 'N' on exit.
107: *
108: * On exit, if EQUED .ne. 'N', A is scaled as follows:
109: * EQUED = 'R': A := diag(R) * A
110: * EQUED = 'C': A := A * diag(C)
111: * EQUED = 'B': A := diag(R) * A * diag(C).
112: *
113: * LDA (input) INTEGER
114: * The leading dimension of the array A. LDA >= max(1,N).
115: *
116: * AF (input or output) COMPLEX*16 array, dimension (LDAF,N)
117: * If FACT = 'F', then AF is an input argument and on entry
118: * contains the factors L and U from the factorization
119: * A = P*L*U as computed by ZGETRF. If EQUED .ne. 'N', then
120: * AF is the factored form of the equilibrated matrix A.
121: *
122: * If FACT = 'N', then AF is an output argument and on exit
123: * returns the factors L and U from the factorization A = P*L*U
124: * of the original matrix A.
125: *
126: * If FACT = 'E', then AF is an output argument and on exit
127: * returns the factors L and U from the factorization A = P*L*U
128: * of the equilibrated matrix A (see the description of A for
129: * the form of the equilibrated matrix).
130: *
131: * LDAF (input) INTEGER
132: * The leading dimension of the array AF. LDAF >= max(1,N).
133: *
134: * IPIV (input or output) INTEGER array, dimension (N)
135: * If FACT = 'F', then IPIV is an input argument and on entry
136: * contains the pivot indices from the factorization A = P*L*U
137: * as computed by ZGETRF; row i of the matrix was interchanged
138: * with row IPIV(i).
139: *
140: * If FACT = 'N', then IPIV is an output argument and on exit
141: * contains the pivot indices from the factorization A = P*L*U
142: * of the original matrix A.
143: *
144: * If FACT = 'E', then IPIV is an output argument and on exit
145: * contains the pivot indices from the factorization A = P*L*U
146: * of the equilibrated matrix A.
147: *
148: * EQUED (input or output) CHARACTER*1
149: * Specifies the form of equilibration that was done.
150: * = 'N': No equilibration (always true if FACT = 'N').
151: * = 'R': Row equilibration, i.e., A has been premultiplied by
152: * diag(R).
153: * = 'C': Column equilibration, i.e., A has been postmultiplied
154: * by diag(C).
155: * = 'B': Both row and column equilibration, i.e., A has been
156: * replaced by diag(R) * A * diag(C).
157: * EQUED is an input argument if FACT = 'F'; otherwise, it is an
158: * output argument.
159: *
160: * R (input or output) DOUBLE PRECISION array, dimension (N)
161: * The row scale factors for A. If EQUED = 'R' or 'B', A is
162: * multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
163: * is not accessed. R is an input argument if FACT = 'F';
164: * otherwise, R is an output argument. If FACT = 'F' and
165: * EQUED = 'R' or 'B', each element of R must be positive.
166: *
167: * C (input or output) DOUBLE PRECISION array, dimension (N)
168: * The column scale factors for A. If EQUED = 'C' or 'B', A is
169: * multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
170: * is not accessed. C is an input argument if FACT = 'F';
171: * otherwise, C is an output argument. If FACT = 'F' and
172: * EQUED = 'C' or 'B', each element of C must be positive.
173: *
174: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
175: * On entry, the N-by-NRHS right hand side matrix B.
176: * On exit,
177: * if EQUED = 'N', B is not modified;
178: * if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
179: * diag(R)*B;
180: * if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
181: * overwritten by diag(C)*B.
182: *
183: * LDB (input) INTEGER
184: * The leading dimension of the array B. LDB >= max(1,N).
185: *
186: * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
187: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X
188: * to the original system of equations. Note that A and B are
189: * modified on exit if EQUED .ne. 'N', and the solution to the
190: * equilibrated system is inv(diag(C))*X if TRANS = 'N' and
191: * EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
192: * and EQUED = 'R' or 'B'.
193: *
194: * LDX (input) INTEGER
195: * The leading dimension of the array X. LDX >= max(1,N).
196: *
197: * RCOND (output) DOUBLE PRECISION
198: * The estimate of the reciprocal condition number of the matrix
199: * A after equilibration (if done). If RCOND is less than the
200: * machine precision (in particular, if RCOND = 0), the matrix
201: * is singular to working precision. This condition is
202: * indicated by a return code of INFO > 0.
203: *
204: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
205: * The estimated forward error bound for each solution vector
206: * X(j) (the j-th column of the solution matrix X).
207: * If XTRUE is the true solution corresponding to X(j), FERR(j)
208: * is an estimated upper bound for the magnitude of the largest
209: * element in (X(j) - XTRUE) divided by the magnitude of the
210: * largest element in X(j). The estimate is as reliable as
211: * the estimate for RCOND, and is almost always a slight
212: * overestimate of the true error.
213: *
214: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
215: * The componentwise relative backward error of each solution
216: * vector X(j) (i.e., the smallest relative change in
217: * any element of A or B that makes X(j) an exact solution).
218: *
219: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
220: *
221: * RWORK (workspace/output) DOUBLE PRECISION array, dimension (2*N)
222: * On exit, RWORK(1) contains the reciprocal pivot growth
223: * factor norm(A)/norm(U). The "max absolute element" norm is
224: * used. If RWORK(1) is much less than 1, then the stability
225: * of the LU factorization of the (equilibrated) matrix A
226: * could be poor. This also means that the solution X, condition
227: * estimator RCOND, and forward error bound FERR could be
228: * unreliable. If factorization fails with 0<INFO<=N, then
229: * RWORK(1) contains the reciprocal pivot growth factor for the
230: * leading INFO columns of A.
231: *
232: * INFO (output) INTEGER
233: * = 0: successful exit
234: * < 0: if INFO = -i, the i-th argument had an illegal value
235: * > 0: if INFO = i, and i is
236: * <= N: U(i,i) is exactly zero. The factorization has
237: * been completed, but the factor U is exactly
238: * singular, so the solution and error bounds
239: * could not be computed. RCOND = 0 is returned.
240: * = N+1: U is nonsingular, but RCOND is less than machine
241: * precision, meaning that the matrix is singular
242: * to working precision. Nevertheless, the
243: * solution and error bounds are computed because
244: * there are a number of situations where the
245: * computed solution can be more accurate than the
246: * value of RCOND would suggest.
247: *
248: * =====================================================================
249: *
250: * .. Parameters ..
251: DOUBLE PRECISION ZERO, ONE
252: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
253: * ..
254: * .. Local Scalars ..
255: LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
256: CHARACTER NORM
257: INTEGER I, INFEQU, J
258: DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
259: $ ROWCND, RPVGRW, SMLNUM
260: * ..
261: * .. External Functions ..
262: LOGICAL LSAME
263: DOUBLE PRECISION DLAMCH, ZLANGE, ZLANTR
264: EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANTR
265: * ..
266: * .. External Subroutines ..
267: EXTERNAL XERBLA, ZGECON, ZGEEQU, ZGERFS, ZGETRF, ZGETRS,
268: $ ZLACPY, ZLAQGE
269: * ..
270: * .. Intrinsic Functions ..
271: INTRINSIC MAX, MIN
272: * ..
273: * .. Executable Statements ..
274: *
275: INFO = 0
276: NOFACT = LSAME( FACT, 'N' )
277: EQUIL = LSAME( FACT, 'E' )
278: NOTRAN = LSAME( TRANS, 'N' )
279: IF( NOFACT .OR. EQUIL ) THEN
280: EQUED = 'N'
281: ROWEQU = .FALSE.
282: COLEQU = .FALSE.
283: ELSE
284: ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
285: COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
286: SMLNUM = DLAMCH( 'Safe minimum' )
287: BIGNUM = ONE / SMLNUM
288: END IF
289: *
290: * Test the input parameters.
291: *
292: IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
293: $ THEN
294: INFO = -1
295: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
296: $ LSAME( TRANS, 'C' ) ) THEN
297: INFO = -2
298: ELSE IF( N.LT.0 ) THEN
299: INFO = -3
300: ELSE IF( NRHS.LT.0 ) THEN
301: INFO = -4
302: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
303: INFO = -6
304: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
305: INFO = -8
306: ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
307: $ ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
308: INFO = -10
309: ELSE
310: IF( ROWEQU ) THEN
311: RCMIN = BIGNUM
312: RCMAX = ZERO
313: DO 10 J = 1, N
314: RCMIN = MIN( RCMIN, R( J ) )
315: RCMAX = MAX( RCMAX, R( J ) )
316: 10 CONTINUE
317: IF( RCMIN.LE.ZERO ) THEN
318: INFO = -11
319: ELSE IF( N.GT.0 ) THEN
320: ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
321: ELSE
322: ROWCND = ONE
323: END IF
324: END IF
325: IF( COLEQU .AND. INFO.EQ.0 ) THEN
326: RCMIN = BIGNUM
327: RCMAX = ZERO
328: DO 20 J = 1, N
329: RCMIN = MIN( RCMIN, C( J ) )
330: RCMAX = MAX( RCMAX, C( J ) )
331: 20 CONTINUE
332: IF( RCMIN.LE.ZERO ) THEN
333: INFO = -12
334: ELSE IF( N.GT.0 ) THEN
335: COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
336: ELSE
337: COLCND = ONE
338: END IF
339: END IF
340: IF( INFO.EQ.0 ) THEN
341: IF( LDB.LT.MAX( 1, N ) ) THEN
342: INFO = -14
343: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
344: INFO = -16
345: END IF
346: END IF
347: END IF
348: *
349: IF( INFO.NE.0 ) THEN
350: CALL XERBLA( 'ZGESVX', -INFO )
351: RETURN
352: END IF
353: *
354: IF( EQUIL ) THEN
355: *
356: * Compute row and column scalings to equilibrate the matrix A.
357: *
358: CALL ZGEEQU( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFEQU )
359: IF( INFEQU.EQ.0 ) THEN
360: *
361: * Equilibrate the matrix.
362: *
363: CALL ZLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
364: $ EQUED )
365: ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
366: COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
367: END IF
368: END IF
369: *
370: * Scale the right hand side.
371: *
372: IF( NOTRAN ) THEN
373: IF( ROWEQU ) THEN
374: DO 40 J = 1, NRHS
375: DO 30 I = 1, N
376: B( I, J ) = R( I )*B( I, J )
377: 30 CONTINUE
378: 40 CONTINUE
379: END IF
380: ELSE IF( COLEQU ) THEN
381: DO 60 J = 1, NRHS
382: DO 50 I = 1, N
383: B( I, J ) = C( I )*B( I, J )
384: 50 CONTINUE
385: 60 CONTINUE
386: END IF
387: *
388: IF( NOFACT .OR. EQUIL ) THEN
389: *
390: * Compute the LU factorization of A.
391: *
392: CALL ZLACPY( 'Full', N, N, A, LDA, AF, LDAF )
393: CALL ZGETRF( N, N, AF, LDAF, IPIV, INFO )
394: *
395: * Return if INFO is non-zero.
396: *
397: IF( INFO.GT.0 ) THEN
398: *
399: * Compute the reciprocal pivot growth factor of the
400: * leading rank-deficient INFO columns of A.
401: *
402: RPVGRW = ZLANTR( 'M', 'U', 'N', INFO, INFO, AF, LDAF,
403: $ RWORK )
404: IF( RPVGRW.EQ.ZERO ) THEN
405: RPVGRW = ONE
406: ELSE
407: RPVGRW = ZLANGE( 'M', N, INFO, A, LDA, RWORK ) /
408: $ RPVGRW
409: END IF
410: RWORK( 1 ) = RPVGRW
411: RCOND = ZERO
412: RETURN
413: END IF
414: END IF
415: *
416: * Compute the norm of the matrix A and the
417: * reciprocal pivot growth factor RPVGRW.
418: *
419: IF( NOTRAN ) THEN
420: NORM = '1'
421: ELSE
422: NORM = 'I'
423: END IF
424: ANORM = ZLANGE( NORM, N, N, A, LDA, RWORK )
425: RPVGRW = ZLANTR( 'M', 'U', 'N', N, N, AF, LDAF, RWORK )
426: IF( RPVGRW.EQ.ZERO ) THEN
427: RPVGRW = ONE
428: ELSE
429: RPVGRW = ZLANGE( 'M', N, N, A, LDA, RWORK ) / RPVGRW
430: END IF
431: *
432: * Compute the reciprocal of the condition number of A.
433: *
434: CALL ZGECON( NORM, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO )
435: *
436: * Compute the solution matrix X.
437: *
438: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
439: CALL ZGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
440: *
441: * Use iterative refinement to improve the computed solution and
442: * compute error bounds and backward error estimates for it.
443: *
444: CALL ZGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
445: $ LDX, FERR, BERR, WORK, RWORK, INFO )
446: *
447: * Transform the solution matrix X to a solution of the original
448: * system.
449: *
450: IF( NOTRAN ) THEN
451: IF( COLEQU ) THEN
452: DO 80 J = 1, NRHS
453: DO 70 I = 1, N
454: X( I, J ) = C( I )*X( I, J )
455: 70 CONTINUE
456: 80 CONTINUE
457: DO 90 J = 1, NRHS
458: FERR( J ) = FERR( J ) / COLCND
459: 90 CONTINUE
460: END IF
461: ELSE IF( ROWEQU ) THEN
462: DO 110 J = 1, NRHS
463: DO 100 I = 1, N
464: X( I, J ) = R( I )*X( I, J )
465: 100 CONTINUE
466: 110 CONTINUE
467: DO 120 J = 1, NRHS
468: FERR( J ) = FERR( J ) / ROWCND
469: 120 CONTINUE
470: END IF
471: *
472: * Set INFO = N+1 if the matrix is singular to working precision.
473: *
474: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
475: $ INFO = N + 1
476: *
477: RWORK( 1 ) = RPVGRW
478: RETURN
479: *
480: * End of ZGESVX
481: *
482: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgesvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack