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