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.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 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**H (or U**H*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
2 $ RCOND, FERR, BERR, WORK, RWORK, 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 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**H (or U**H*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