1: SUBROUTINE ZPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED,
2: $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
3: $ 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, UPLO
12: INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
13: DOUBLE PRECISION RCOND
14: * ..
15: * .. Array Arguments ..
16: DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * )
17: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
18: $ WORK( * ), X( LDX, * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * ZPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
25: * compute the solution to a complex system of linear equations
26: * A * X = B,
27: * where A is an N-by-N Hermitian positive definite matrix and X and B
28: * are N-by-NRHS matrices.
29: *
30: * Error bounds on the solution and a condition estimate are also
31: * provided.
32: *
33: * Description
34: * ===========
35: *
36: * The following steps are performed:
37: *
38: * 1. If FACT = 'E', real scaling factors are computed to equilibrate
39: * the system:
40: * diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
41: * Whether or not the system will be equilibrated depends on the
42: * scaling of the matrix A, but if equilibration is used, A is
43: * overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
44: *
45: * 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
46: * factor the matrix A (after equilibration if FACT = 'E') as
47: * A = U**H* U, if UPLO = 'U', or
48: * A = L * L**H, if UPLO = 'L',
49: * where U is an upper triangular matrix and L is a lower triangular
50: * matrix.
51: *
52: * 3. If the leading i-by-i principal minor is not positive definite,
53: * then the routine returns with INFO = i. Otherwise, the factored
54: * form of A is used to estimate the condition number of the matrix
55: * A. If the reciprocal of the condition number is less than machine
56: * precision, INFO = N+1 is returned as a warning, but the routine
57: * still goes on to solve for X and compute error bounds as
58: * described below.
59: *
60: * 4. The system of equations is solved for X using the factored form
61: * of A.
62: *
63: * 5. Iterative refinement is applied to improve the computed solution
64: * matrix and calculate error bounds and backward error estimates
65: * for it.
66: *
67: * 6. If equilibration was used, the matrix X is premultiplied by
68: * diag(S) so that it solves the original system before
69: * equilibration.
70: *
71: * Arguments
72: * =========
73: *
74: * FACT (input) CHARACTER*1
75: * Specifies whether or not the factored form of the matrix A is
76: * supplied on entry, and if not, whether the matrix A should be
77: * equilibrated before it is factored.
78: * = 'F': On entry, AF contains the factored form of A.
79: * If EQUED = 'Y', the matrix A has been equilibrated
80: * with scaling factors given by S. A and AF will not
81: * be modified.
82: * = 'N': The matrix A will be copied to AF and factored.
83: * = 'E': The matrix A will be equilibrated if necessary, then
84: * copied to AF and factored.
85: *
86: * UPLO (input) CHARACTER*1
87: * = 'U': Upper triangle of A is stored;
88: * = 'L': Lower triangle of A is stored.
89: *
90: * N (input) INTEGER
91: * The number of linear equations, i.e., the order of the
92: * matrix A. N >= 0.
93: *
94: * NRHS (input) INTEGER
95: * The number of right hand sides, i.e., the number of columns
96: * of the matrices B and X. NRHS >= 0.
97: *
98: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
99: * On entry, the Hermitian matrix A, except if FACT = 'F' and
100: * EQUED = 'Y', then A must contain the equilibrated matrix
101: * diag(S)*A*diag(S). If UPLO = 'U', the leading
102: * N-by-N upper triangular part of A contains the upper
103: * triangular part of the matrix A, and the strictly lower
104: * triangular part of A is not referenced. If UPLO = 'L', the
105: * leading N-by-N lower triangular part of A contains the lower
106: * triangular part of the matrix A, and the strictly upper
107: * triangular part of A is not referenced. A is not modified if
108: * FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
109: *
110: * On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
111: * diag(S)*A*diag(S).
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 triangular factor U or L from the Cholesky
119: * factorization A = U**H*U or A = L*L**H, in the same storage
120: * format as A. If EQUED .ne. 'N', then AF is the factored form
121: * of the equilibrated matrix diag(S)*A*diag(S).
122: *
123: * If FACT = 'N', then AF is an output argument and on exit
124: * returns the triangular factor U or L from the Cholesky
125: * factorization A = U**H*U or A = L*L**H of the original
126: * matrix A.
127: *
128: * If FACT = 'E', then AF is an output argument and on exit
129: * returns the triangular factor U or L from the Cholesky
130: * factorization A = U**H*U or A = L*L**H of the equilibrated
131: * matrix A (see the description of A for the form of the
132: * equilibrated matrix).
133: *
134: * LDAF (input) INTEGER
135: * The leading dimension of the array AF. LDAF >= max(1,N).
136: *
137: * EQUED (input or output) CHARACTER*1
138: * Specifies the form of equilibration that was done.
139: * = 'N': No equilibration (always true if FACT = 'N').
140: * = 'Y': Equilibration was done, i.e., A has been replaced by
141: * diag(S) * A * diag(S).
142: * EQUED is an input argument if FACT = 'F'; otherwise, it is an
143: * output argument.
144: *
145: * S (input or output) DOUBLE PRECISION array, dimension (N)
146: * The scale factors for A; not accessed if EQUED = 'N'. S is
147: * an input argument if FACT = 'F'; otherwise, S is an output
148: * argument. If FACT = 'F' and EQUED = 'Y', each element of S
149: * must be positive.
150: *
151: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
152: * On entry, the N-by-NRHS righthand side matrix B.
153: * On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
154: * B is overwritten by diag(S) * B.
155: *
156: * LDB (input) INTEGER
157: * The leading dimension of the array B. LDB >= max(1,N).
158: *
159: * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
160: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
161: * the original system of equations. Note that if EQUED = 'Y',
162: * A and B are modified on exit, and the solution to the
163: * equilibrated system is inv(diag(S))*X.
164: *
165: * LDX (input) INTEGER
166: * The leading dimension of the array X. LDX >= max(1,N).
167: *
168: * RCOND (output) DOUBLE PRECISION
169: * The estimate of the reciprocal condition number of the matrix
170: * A after equilibration (if done). If RCOND is less than the
171: * machine precision (in particular, if RCOND = 0), the matrix
172: * is singular to working precision. This condition is
173: * indicated by a return code of INFO > 0.
174: *
175: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
176: * The estimated forward error bound for each solution vector
177: * X(j) (the j-th column of the solution matrix X).
178: * If XTRUE is the true solution corresponding to X(j), FERR(j)
179: * is an estimated upper bound for the magnitude of the largest
180: * element in (X(j) - XTRUE) divided by the magnitude of the
181: * largest element in X(j). The estimate is as reliable as
182: * the estimate for RCOND, and is almost always a slight
183: * overestimate of the true error.
184: *
185: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
186: * The componentwise relative backward error of each solution
187: * vector X(j) (i.e., the smallest relative change in
188: * any element of A or B that makes X(j) an exact solution).
189: *
190: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
191: *
192: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
193: *
194: * INFO (output) INTEGER
195: * = 0: successful exit
196: * < 0: if INFO = -i, the i-th argument had an illegal value
197: * > 0: if INFO = i, and i is
198: * <= N: the leading minor of order i of A is
199: * not positive definite, so the factorization
200: * could not be completed, and the solution has not
201: * been computed. RCOND = 0 is returned.
202: * = N+1: U is nonsingular, but RCOND is less than machine
203: * precision, meaning that the matrix is singular
204: * to working precision. Nevertheless, the
205: * solution and error bounds are computed because
206: * there are a number of situations where the
207: * computed solution can be more accurate than the
208: * value of RCOND would suggest.
209: *
210: * =====================================================================
211: *
212: * .. Parameters ..
213: DOUBLE PRECISION ZERO, ONE
214: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
215: * ..
216: * .. Local Scalars ..
217: LOGICAL EQUIL, NOFACT, RCEQU
218: INTEGER I, INFEQU, J
219: DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
220: * ..
221: * .. External Functions ..
222: LOGICAL LSAME
223: DOUBLE PRECISION DLAMCH, ZLANHE
224: EXTERNAL LSAME, DLAMCH, ZLANHE
225: * ..
226: * .. External Subroutines ..
227: EXTERNAL XERBLA, ZLACPY, ZLAQHE, ZPOCON, ZPOEQU, ZPORFS,
228: $ ZPOTRF, ZPOTRS
229: * ..
230: * .. Intrinsic Functions ..
231: INTRINSIC MAX, MIN
232: * ..
233: * .. Executable Statements ..
234: *
235: INFO = 0
236: NOFACT = LSAME( FACT, 'N' )
237: EQUIL = LSAME( FACT, 'E' )
238: IF( NOFACT .OR. EQUIL ) THEN
239: EQUED = 'N'
240: RCEQU = .FALSE.
241: ELSE
242: RCEQU = LSAME( EQUED, 'Y' )
243: SMLNUM = DLAMCH( 'Safe minimum' )
244: BIGNUM = ONE / SMLNUM
245: END IF
246: *
247: * Test the input parameters.
248: *
249: IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
250: $ THEN
251: INFO = -1
252: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
253: $ THEN
254: INFO = -2
255: ELSE IF( N.LT.0 ) THEN
256: INFO = -3
257: ELSE IF( NRHS.LT.0 ) THEN
258: INFO = -4
259: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
260: INFO = -6
261: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
262: INFO = -8
263: ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
264: $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
265: INFO = -9
266: ELSE
267: IF( RCEQU ) THEN
268: SMIN = BIGNUM
269: SMAX = ZERO
270: DO 10 J = 1, N
271: SMIN = MIN( SMIN, S( J ) )
272: SMAX = MAX( SMAX, S( J ) )
273: 10 CONTINUE
274: IF( SMIN.LE.ZERO ) THEN
275: INFO = -10
276: ELSE IF( N.GT.0 ) THEN
277: SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
278: ELSE
279: SCOND = ONE
280: END IF
281: END IF
282: IF( INFO.EQ.0 ) THEN
283: IF( LDB.LT.MAX( 1, N ) ) THEN
284: INFO = -12
285: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
286: INFO = -14
287: END IF
288: END IF
289: END IF
290: *
291: IF( INFO.NE.0 ) THEN
292: CALL XERBLA( 'ZPOSVX', -INFO )
293: RETURN
294: END IF
295: *
296: IF( EQUIL ) THEN
297: *
298: * Compute row and column scalings to equilibrate the matrix A.
299: *
300: CALL ZPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU )
301: IF( INFEQU.EQ.0 ) THEN
302: *
303: * Equilibrate the matrix.
304: *
305: CALL ZLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED )
306: RCEQU = LSAME( EQUED, 'Y' )
307: END IF
308: END IF
309: *
310: * Scale the right hand side.
311: *
312: IF( RCEQU ) THEN
313: DO 30 J = 1, NRHS
314: DO 20 I = 1, N
315: B( I, J ) = S( I )*B( I, J )
316: 20 CONTINUE
317: 30 CONTINUE
318: END IF
319: *
320: IF( NOFACT .OR. EQUIL ) THEN
321: *
322: * Compute the Cholesky factorization A = U'*U or A = L*L'.
323: *
324: CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF )
325: CALL ZPOTRF( UPLO, N, AF, LDAF, INFO )
326: *
327: * Return if INFO is non-zero.
328: *
329: IF( INFO.GT.0 )THEN
330: RCOND = ZERO
331: RETURN
332: END IF
333: END IF
334: *
335: * Compute the norm of the matrix A.
336: *
337: ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
338: *
339: * Compute the reciprocal of the condition number of A.
340: *
341: CALL ZPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO )
342: *
343: * Compute the solution matrix X.
344: *
345: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
346: CALL ZPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO )
347: *
348: * Use iterative refinement to improve the computed solution and
349: * compute error bounds and backward error estimates for it.
350: *
351: CALL ZPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX,
352: $ FERR, BERR, WORK, RWORK, INFO )
353: *
354: * Transform the solution matrix X to a solution of the original
355: * system.
356: *
357: IF( RCEQU ) THEN
358: DO 50 J = 1, NRHS
359: DO 40 I = 1, N
360: X( I, J ) = S( I )*X( I, J )
361: 40 CONTINUE
362: 50 CONTINUE
363: DO 60 J = 1, NRHS
364: FERR( J ) = FERR( J ) / SCOND
365: 60 CONTINUE
366: END IF
367: *
368: * Set INFO = N+1 if the matrix is singular to working precision.
369: *
370: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
371: $ INFO = N + 1
372: *
373: RETURN
374: *
375: * End of ZPOSVX
376: *
377: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zposvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack