Annotation of rpl/lapack/lapack/zptsvx.f, revision 1.3
1.1 bertrand 1: SUBROUTINE ZPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
2: $ RCOND, FERR, BERR, WORK, RWORK, INFO )
3: *
4: * -- LAPACK 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 FACT
11: INTEGER INFO, LDB, LDX, N, NRHS
12: DOUBLE PRECISION RCOND
13: * ..
14: * .. Array Arguments ..
15: DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
16: $ RWORK( * )
17: COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
18: $ X( LDX, * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * ZPTSVX uses the factorization A = L*D*L**H to compute the solution
25: * to a complex system of linear equations A*X = B, where A is an
26: * N-by-N Hermitian positive definite tridiagonal matrix and X and B
27: * 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 = 'N', the matrix A is factored as A = L*D*L**H, where L
38: * is a unit lower bidiagonal matrix and D is diagonal. The
39: * factorization can also be regarded as having the form
40: * A = U**H*D*U.
41: *
42: * 2. If the leading i-by-i principal minor is not positive definite,
43: * then the routine returns with INFO = i. Otherwise, the factored
44: * form of A is used to estimate the condition number of the matrix
45: * A. If the reciprocal of the condition number is less than machine
46: * precision, INFO = N+1 is returned as a warning, but the routine
47: * still goes on to solve for X and compute error bounds as
48: * described below.
49: *
50: * 3. The system of equations is solved for X using the factored form
51: * of A.
52: *
53: * 4. Iterative refinement is applied to improve the computed solution
54: * matrix and calculate error bounds and backward error estimates
55: * for it.
56: *
57: * Arguments
58: * =========
59: *
60: * FACT (input) CHARACTER*1
61: * Specifies whether or not the factored form of the matrix
62: * A is supplied on entry.
63: * = 'F': On entry, DF and EF contain the factored form of A.
64: * D, E, DF, and EF will not be modified.
65: * = 'N': The matrix A will be copied to DF and EF and
66: * factored.
67: *
68: * N (input) INTEGER
69: * The order of the matrix A. N >= 0.
70: *
71: * NRHS (input) INTEGER
72: * The number of right hand sides, i.e., the number of columns
73: * of the matrices B and X. NRHS >= 0.
74: *
75: * D (input) DOUBLE PRECISION array, dimension (N)
76: * The n diagonal elements of the tridiagonal matrix A.
77: *
78: * E (input) COMPLEX*16 array, dimension (N-1)
79: * The (n-1) subdiagonal elements of the tridiagonal matrix A.
80: *
81: * DF (input or output) DOUBLE PRECISION array, dimension (N)
82: * If FACT = 'F', then DF is an input argument and on entry
83: * contains the n diagonal elements of the diagonal matrix D
84: * from the L*D*L**H factorization of A.
85: * If FACT = 'N', then DF is an output argument and on exit
86: * contains the n diagonal elements of the diagonal matrix D
87: * from the L*D*L**H factorization of A.
88: *
89: * EF (input or output) COMPLEX*16 array, dimension (N-1)
90: * If FACT = 'F', then EF is an input argument and on entry
91: * contains the (n-1) subdiagonal elements of the unit
92: * bidiagonal factor L from the L*D*L**H factorization of A.
93: * If FACT = 'N', then EF is an output argument and on exit
94: * contains the (n-1) subdiagonal elements of the unit
95: * bidiagonal factor L from the L*D*L**H factorization of A.
96: *
97: * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
98: * The N-by-NRHS right hand side matrix B.
99: *
100: * LDB (input) INTEGER
101: * The leading dimension of the array B. LDB >= max(1,N).
102: *
103: * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
104: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
105: *
106: * LDX (input) INTEGER
107: * The leading dimension of the array X. LDX >= max(1,N).
108: *
109: * RCOND (output) DOUBLE PRECISION
110: * The reciprocal condition number of the matrix A. If RCOND
111: * is less than the machine precision (in particular, if
112: * RCOND = 0), the matrix is singular to working precision.
113: * This condition is indicated by a return code of INFO > 0.
114: *
115: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
116: * The forward error bound for each solution vector
117: * X(j) (the j-th column of the solution matrix X).
118: * If XTRUE is the true solution corresponding to X(j), FERR(j)
119: * is an estimated upper bound for the magnitude of the largest
120: * element in (X(j) - XTRUE) divided by the magnitude of the
121: * largest element in X(j).
122: *
123: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
124: * The componentwise relative backward error of each solution
125: * vector X(j) (i.e., the smallest relative change in any
126: * element of A or B that makes X(j) an exact solution).
127: *
128: * WORK (workspace) COMPLEX*16 array, dimension (N)
129: *
130: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
131: *
132: * INFO (output) INTEGER
133: * = 0: successful exit
134: * < 0: if INFO = -i, the i-th argument had an illegal value
135: * > 0: if INFO = i, and i is
136: * <= N: the leading minor of order i of A is
137: * not positive definite, so the factorization
138: * could not be completed, and the solution has not
139: * been computed. RCOND = 0 is returned.
140: * = N+1: U is nonsingular, but RCOND is less than machine
141: * precision, meaning that the matrix is singular
142: * to working precision. Nevertheless, the
143: * solution and error bounds are computed because
144: * there are a number of situations where the
145: * computed solution can be more accurate than the
146: * value of RCOND would suggest.
147: *
148: * =====================================================================
149: *
150: * .. Parameters ..
151: DOUBLE PRECISION ZERO
152: PARAMETER ( ZERO = 0.0D+0 )
153: * ..
154: * .. Local Scalars ..
155: LOGICAL NOFACT
156: DOUBLE PRECISION ANORM
157: * ..
158: * .. External Functions ..
159: LOGICAL LSAME
160: DOUBLE PRECISION DLAMCH, ZLANHT
161: EXTERNAL LSAME, DLAMCH, ZLANHT
162: * ..
163: * .. External Subroutines ..
164: EXTERNAL DCOPY, XERBLA, ZCOPY, ZLACPY, ZPTCON, ZPTRFS,
165: $ ZPTTRF, ZPTTRS
166: * ..
167: * .. Intrinsic Functions ..
168: INTRINSIC MAX
169: * ..
170: * .. Executable Statements ..
171: *
172: * Test the input parameters.
173: *
174: INFO = 0
175: NOFACT = LSAME( FACT, 'N' )
176: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
177: INFO = -1
178: ELSE IF( N.LT.0 ) THEN
179: INFO = -2
180: ELSE IF( NRHS.LT.0 ) THEN
181: INFO = -3
182: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
183: INFO = -9
184: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
185: INFO = -11
186: END IF
187: IF( INFO.NE.0 ) THEN
188: CALL XERBLA( 'ZPTSVX', -INFO )
189: RETURN
190: END IF
191: *
192: IF( NOFACT ) THEN
193: *
194: * Compute the L*D*L' (or U'*D*U) factorization of A.
195: *
196: CALL DCOPY( N, D, 1, DF, 1 )
197: IF( N.GT.1 )
198: $ CALL ZCOPY( N-1, E, 1, EF, 1 )
199: CALL ZPTTRF( N, DF, EF, INFO )
200: *
201: * Return if INFO is non-zero.
202: *
203: IF( INFO.GT.0 )THEN
204: RCOND = ZERO
205: RETURN
206: END IF
207: END IF
208: *
209: * Compute the norm of the matrix A.
210: *
211: ANORM = ZLANHT( '1', N, D, E )
212: *
213: * Compute the reciprocal of the condition number of A.
214: *
215: CALL ZPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO )
216: *
217: * Compute the solution vectors X.
218: *
219: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
220: CALL ZPTTRS( 'Lower', N, NRHS, DF, EF, X, LDX, INFO )
221: *
222: * Use iterative refinement to improve the computed solutions and
223: * compute error bounds and backward error estimates for them.
224: *
225: CALL ZPTRFS( 'Lower', N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
226: $ BERR, WORK, RWORK, INFO )
227: *
228: * Set INFO = N+1 if the matrix is singular to working precision.
229: *
230: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
231: $ INFO = N + 1
232: *
233: RETURN
234: *
235: * End of ZPTSVX
236: *
237: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zptsvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack