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.3.1) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * -- April 2011 --
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**T*U or A = L*L**T, 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**T * U or A = L * L**T of the original
122: * matrix A.
123: *
124: * If FACT = 'E', then AFP is an output argument and on exit
125: * returns the triangular factor U or L from the Cholesky
126: * factorization A = U**T * U or A = L * L**T of the equilibrated
127: * matrix A (see the description of AP for the form of the
128: * equilibrated matrix).
129: *
130: * EQUED (input or output) CHARACTER*1
131: * Specifies the form of equilibration that was done.
132: * = 'N': No equilibration (always true if FACT = 'N').
133: * = 'Y': Equilibration was done, i.e., A has been replaced by
134: * diag(S) * A * diag(S).
135: * EQUED is an input argument if FACT = 'F'; otherwise, it is an
136: * output argument.
137: *
138: * S (input or output) DOUBLE PRECISION array, dimension (N)
139: * The scale factors for A; not accessed if EQUED = 'N'. S is
140: * an input argument if FACT = 'F'; otherwise, S is an output
141: * argument. If FACT = 'F' and EQUED = 'Y', each element of S
142: * must be positive.
143: *
144: * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
145: * On entry, the N-by-NRHS right hand side matrix B.
146: * On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
147: * B is overwritten by diag(S) * B.
148: *
149: * LDB (input) INTEGER
150: * The leading dimension of the array B. LDB >= max(1,N).
151: *
152: * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
153: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
154: * the original system of equations. Note that if EQUED = 'Y',
155: * A and B are modified on exit, and the solution to the
156: * equilibrated system is inv(diag(S))*X.
157: *
158: * LDX (input) INTEGER
159: * The leading dimension of the array X. LDX >= max(1,N).
160: *
161: * RCOND (output) DOUBLE PRECISION
162: * The estimate of the reciprocal condition number of the matrix
163: * A after equilibration (if done). If RCOND is less than the
164: * machine precision (in particular, if RCOND = 0), the matrix
165: * is singular to working precision. This condition is
166: * indicated by a return code of INFO > 0.
167: *
168: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
169: * The estimated forward error bound for each solution vector
170: * X(j) (the j-th column of the solution matrix X).
171: * If XTRUE is the true solution corresponding to X(j), FERR(j)
172: * is an estimated upper bound for the magnitude of the largest
173: * element in (X(j) - XTRUE) divided by the magnitude of the
174: * largest element in X(j). The estimate is as reliable as
175: * the estimate for RCOND, and is almost always a slight
176: * overestimate of the true error.
177: *
178: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
179: * The componentwise relative backward error of each solution
180: * vector X(j) (i.e., the smallest relative change in
181: * any element of A or B that makes X(j) an exact solution).
182: *
183: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
184: *
185: * IWORK (workspace) INTEGER array, dimension (N)
186: *
187: * INFO (output) INTEGER
188: * = 0: successful exit
189: * < 0: if INFO = -i, the i-th argument had an illegal value
190: * > 0: if INFO = i, and i is
191: * <= N: the leading minor of order i of A is
192: * not positive definite, so the factorization
193: * could not be completed, and the solution has not
194: * been computed. RCOND = 0 is returned.
195: * = N+1: U is nonsingular, but RCOND is less than machine
196: * precision, meaning that the matrix is singular
197: * to working precision. Nevertheless, the
198: * solution and error bounds are computed because
199: * there are a number of situations where the
200: * computed solution can be more accurate than the
201: * value of RCOND would suggest.
202: *
203: * Further Details
204: * ===============
205: *
206: * The packed storage scheme is illustrated by the following example
207: * when N = 4, UPLO = 'U':
208: *
209: * Two-dimensional storage of the symmetric matrix A:
210: *
211: * a11 a12 a13 a14
212: * a22 a23 a24
213: * a33 a34 (aij = conjg(aji))
214: * a44
215: *
216: * Packed storage of the upper triangle of A:
217: *
218: * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
219: *
220: * =====================================================================
221: *
222: * .. Parameters ..
223: DOUBLE PRECISION ZERO, ONE
224: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
225: * ..
226: * .. Local Scalars ..
227: LOGICAL EQUIL, NOFACT, RCEQU
228: INTEGER I, INFEQU, J
229: DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
230: * ..
231: * .. External Functions ..
232: LOGICAL LSAME
233: DOUBLE PRECISION DLAMCH, DLANSP
234: EXTERNAL LSAME, DLAMCH, DLANSP
235: * ..
236: * .. External Subroutines ..
237: EXTERNAL DCOPY, DLACPY, DLAQSP, DPPCON, DPPEQU, DPPRFS,
238: $ DPPTRF, DPPTRS, XERBLA
239: * ..
240: * .. Intrinsic Functions ..
241: INTRINSIC MAX, MIN
242: * ..
243: * .. Executable Statements ..
244: *
245: INFO = 0
246: NOFACT = LSAME( FACT, 'N' )
247: EQUIL = LSAME( FACT, 'E' )
248: IF( NOFACT .OR. EQUIL ) THEN
249: EQUED = 'N'
250: RCEQU = .FALSE.
251: ELSE
252: RCEQU = LSAME( EQUED, 'Y' )
253: SMLNUM = DLAMCH( 'Safe minimum' )
254: BIGNUM = ONE / SMLNUM
255: END IF
256: *
257: * Test the input parameters.
258: *
259: IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
260: $ THEN
261: INFO = -1
262: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
263: $ THEN
264: INFO = -2
265: ELSE IF( N.LT.0 ) THEN
266: INFO = -3
267: ELSE IF( NRHS.LT.0 ) THEN
268: INFO = -4
269: ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
270: $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
271: INFO = -7
272: ELSE
273: IF( RCEQU ) THEN
274: SMIN = BIGNUM
275: SMAX = ZERO
276: DO 10 J = 1, N
277: SMIN = MIN( SMIN, S( J ) )
278: SMAX = MAX( SMAX, S( J ) )
279: 10 CONTINUE
280: IF( SMIN.LE.ZERO ) THEN
281: INFO = -8
282: ELSE IF( N.GT.0 ) THEN
283: SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
284: ELSE
285: SCOND = ONE
286: END IF
287: END IF
288: IF( INFO.EQ.0 ) THEN
289: IF( LDB.LT.MAX( 1, N ) ) THEN
290: INFO = -10
291: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
292: INFO = -12
293: END IF
294: END IF
295: END IF
296: *
297: IF( INFO.NE.0 ) THEN
298: CALL XERBLA( 'DPPSVX', -INFO )
299: RETURN
300: END IF
301: *
302: IF( EQUIL ) THEN
303: *
304: * Compute row and column scalings to equilibrate the matrix A.
305: *
306: CALL DPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFEQU )
307: IF( INFEQU.EQ.0 ) THEN
308: *
309: * Equilibrate the matrix.
310: *
311: CALL DLAQSP( UPLO, N, AP, S, SCOND, AMAX, EQUED )
312: RCEQU = LSAME( EQUED, 'Y' )
313: END IF
314: END IF
315: *
316: * Scale the right-hand side.
317: *
318: IF( RCEQU ) THEN
319: DO 30 J = 1, NRHS
320: DO 20 I = 1, N
321: B( I, J ) = S( I )*B( I, J )
322: 20 CONTINUE
323: 30 CONTINUE
324: END IF
325: *
326: IF( NOFACT .OR. EQUIL ) THEN
327: *
328: * Compute the Cholesky factorization A = U**T * U or A = L * L**T.
329: *
330: CALL DCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
331: CALL DPPTRF( UPLO, N, AFP, INFO )
332: *
333: * Return if INFO is non-zero.
334: *
335: IF( INFO.GT.0 )THEN
336: RCOND = ZERO
337: RETURN
338: END IF
339: END IF
340: *
341: * Compute the norm of the matrix A.
342: *
343: ANORM = DLANSP( 'I', UPLO, N, AP, WORK )
344: *
345: * Compute the reciprocal of the condition number of A.
346: *
347: CALL DPPCON( UPLO, N, AFP, ANORM, RCOND, WORK, IWORK, INFO )
348: *
349: * Compute the solution matrix X.
350: *
351: CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
352: CALL DPPTRS( UPLO, N, NRHS, AFP, X, LDX, INFO )
353: *
354: * Use iterative refinement to improve the computed solution and
355: * compute error bounds and backward error estimates for it.
356: *
357: CALL DPPRFS( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR, BERR,
358: $ WORK, IWORK, INFO )
359: *
360: * Transform the solution matrix X to a solution of the original
361: * system.
362: *
363: IF( RCEQU ) THEN
364: DO 50 J = 1, NRHS
365: DO 40 I = 1, N
366: X( I, J ) = S( I )*X( I, J )
367: 40 CONTINUE
368: 50 CONTINUE
369: DO 60 J = 1, NRHS
370: FERR( J ) = FERR( J ) / SCOND
371: 60 CONTINUE
372: END IF
373: *
374: * Set INFO = N+1 if the matrix is singular to working precision.
375: *
376: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
377: $ INFO = N + 1
378: *
379: RETURN
380: *
381: * End of DPPSVX
382: *
383: 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