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