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