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