LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
dsysvx.f
Go to the documentation of this file.
1 *> \brief <b> DSYSVX computes the solution to system of linear equations A * X = B for SY matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DSYSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsysvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsysvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsysvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
22 * LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
23 * IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER FACT, UPLO
27 * INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
28 * DOUBLE PRECISION RCOND
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IPIV( * ), IWORK( * )
32 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
33 * $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> DSYSVX uses the diagonal pivoting factorization to compute the
43 *> solution to a real system of linear equations A * X = B,
44 *> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
45 *> matrices.
46 *>
47 *> Error bounds on the solution and a condition estimate are also
48 *> provided.
49 *> \endverbatim
50 *
51 *> \par Description:
52 * =================
53 *>
54 *> \verbatim
55 *>
56 *> The following steps are performed:
57 *>
58 *> 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
59 *> The form of the factorization is
60 *> A = U * D * U**T, if UPLO = 'U', or
61 *> A = L * D * L**T, if UPLO = 'L',
62 *> where U (or L) is a product of permutation and unit upper (lower)
63 *> triangular matrices, and D is symmetric and block diagonal with
64 *> 1-by-1 and 2-by-2 diagonal blocks.
65 *>
66 *> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
67 *> returns with INFO = i. Otherwise, the factored form of A is used
68 *> to estimate the condition number of the matrix A. If the
69 *> reciprocal of the condition number is less than machine precision,
70 *> INFO = N+1 is returned as a warning, but the routine still goes on
71 *> to solve for X and compute error bounds as described below.
72 *>
73 *> 3. The system of equations is solved for X using the factored form
74 *> of A.
75 *>
76 *> 4. Iterative refinement is applied to improve the computed solution
77 *> matrix and calculate error bounds and backward error estimates
78 *> for it.
79 *> \endverbatim
80 *
81 * Arguments:
82 * ==========
83 *
84 *> \param[in] FACT
85 *> \verbatim
86 *> FACT is CHARACTER*1
87 *> Specifies whether or not the factored form of A has been
88 *> supplied on entry.
89 *> = 'F': On entry, AF and IPIV contain the factored form of
90 *> A. AF and IPIV will not be modified.
91 *> = 'N': The matrix A will be copied to AF and factored.
92 *> \endverbatim
93 *>
94 *> \param[in] UPLO
95 *> \verbatim
96 *> UPLO is CHARACTER*1
97 *> = 'U': Upper triangle of A is stored;
98 *> = 'L': Lower triangle of A is stored.
99 *> \endverbatim
100 *>
101 *> \param[in] N
102 *> \verbatim
103 *> N is INTEGER
104 *> The number of linear equations, i.e., the order of the
105 *> matrix A. N >= 0.
106 *> \endverbatim
107 *>
108 *> \param[in] NRHS
109 *> \verbatim
110 *> NRHS is INTEGER
111 *> The number of right hand sides, i.e., the number of columns
112 *> of the matrices B and X. NRHS >= 0.
113 *> \endverbatim
114 *>
115 *> \param[in] A
116 *> \verbatim
117 *> A is DOUBLE PRECISION array, dimension (LDA,N)
118 *> The symmetric matrix A. If UPLO = 'U', the leading N-by-N
119 *> upper triangular part of A contains the upper triangular part
120 *> of the matrix A, and the strictly lower triangular part of A
121 *> is not referenced. If UPLO = 'L', the leading N-by-N lower
122 *> triangular part of A contains the lower triangular part of
123 *> the matrix A, and the strictly upper triangular part of A is
124 *> not referenced.
125 *> \endverbatim
126 *>
127 *> \param[in] LDA
128 *> \verbatim
129 *> LDA is INTEGER
130 *> The leading dimension of the array A. LDA >= max(1,N).
131 *> \endverbatim
132 *>
133 *> \param[in,out] AF
134 *> \verbatim
135 *> AF is DOUBLE PRECISION array, dimension (LDAF,N)
136 *> If FACT = 'F', then AF is an input argument and on entry
137 *> contains the block diagonal matrix D and the multipliers used
138 *> to obtain the factor U or L from the factorization
139 *> A = U*D*U**T or A = L*D*L**T as computed by DSYTRF.
140 *>
141 *> If FACT = 'N', then AF is an output argument and on exit
142 *> returns the block diagonal matrix D and the multipliers used
143 *> to obtain the factor U or L from the factorization
144 *> A = U*D*U**T or A = L*D*L**T.
145 *> \endverbatim
146 *>
147 *> \param[in] LDAF
148 *> \verbatim
149 *> LDAF is INTEGER
150 *> The leading dimension of the array AF. LDAF >= max(1,N).
151 *> \endverbatim
152 *>
153 *> \param[in,out] IPIV
154 *> \verbatim
155 *> IPIV is INTEGER array, dimension (N)
156 *> If FACT = 'F', then IPIV is an input argument and on entry
157 *> contains details of the interchanges and the block structure
158 *> of D, as determined by DSYTRF.
159 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
160 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
161 *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
162 *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
163 *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
164 *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
165 *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
166 *>
167 *> If FACT = 'N', then IPIV is an output argument and on exit
168 *> contains details of the interchanges and the block structure
169 *> of D, as determined by DSYTRF.
170 *> \endverbatim
171 *>
172 *> \param[in] B
173 *> \verbatim
174 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
175 *> The N-by-NRHS right hand side matrix B.
176 *> \endverbatim
177 *>
178 *> \param[in] LDB
179 *> \verbatim
180 *> LDB is INTEGER
181 *> The leading dimension of the array B. LDB >= max(1,N).
182 *> \endverbatim
183 *>
184 *> \param[out] X
185 *> \verbatim
186 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
187 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
188 *> \endverbatim
189 *>
190 *> \param[in] LDX
191 *> \verbatim
192 *> LDX is INTEGER
193 *> The leading dimension of the array X. LDX >= max(1,N).
194 *> \endverbatim
195 *>
196 *> \param[out] RCOND
197 *> \verbatim
198 *> RCOND is DOUBLE PRECISION
199 *> The estimate of the reciprocal condition number of the matrix
200 *> A. If RCOND is less than the machine precision (in
201 *> particular, if RCOND = 0), the matrix is singular to working
202 *> precision. This condition is indicated by a return code of
203 *> INFO > 0.
204 *> \endverbatim
205 *>
206 *> \param[out] FERR
207 *> \verbatim
208 *> FERR is DOUBLE PRECISION array, dimension (NRHS)
209 *> The estimated forward error bound for each solution vector
210 *> X(j) (the j-th column of the solution matrix X).
211 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
212 *> is an estimated upper bound for the magnitude of the largest
213 *> element in (X(j) - XTRUE) divided by the magnitude of the
214 *> largest element in X(j). The estimate is as reliable as
215 *> the estimate for RCOND, and is almost always a slight
216 *> overestimate of the true error.
217 *> \endverbatim
218 *>
219 *> \param[out] BERR
220 *> \verbatim
221 *> BERR is DOUBLE PRECISION array, dimension (NRHS)
222 *> The componentwise relative backward error of each solution
223 *> vector X(j) (i.e., the smallest relative change in
224 *> any element of A or B that makes X(j) an exact solution).
225 *> \endverbatim
226 *>
227 *> \param[out] WORK
228 *> \verbatim
229 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
230 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
231 *> \endverbatim
232 *>
233 *> \param[in] LWORK
234 *> \verbatim
235 *> LWORK is INTEGER
236 *> The length of WORK. LWORK >= max(1,3*N), and for best
237 *> performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where
238 *> NB is the optimal blocksize for DSYTRF.
239 *>
240 *> If LWORK = -1, then a workspace query is assumed; the routine
241 *> only calculates the optimal size of the WORK array, returns
242 *> this value as the first entry of the WORK array, and no error
243 *> message related to LWORK is issued by XERBLA.
244 *> \endverbatim
245 *>
246 *> \param[out] IWORK
247 *> \verbatim
248 *> IWORK is INTEGER array, dimension (N)
249 *> \endverbatim
250 *>
251 *> \param[out] INFO
252 *> \verbatim
253 *> INFO is INTEGER
254 *> = 0: successful exit
255 *> < 0: if INFO = -i, the i-th argument had an illegal value
256 *> > 0: if INFO = i, and i is
257 *> <= N: D(i,i) is exactly zero. The factorization
258 *> has been completed but the factor D is exactly
259 *> singular, so the solution and error bounds could
260 *> not be computed. RCOND = 0 is returned.
261 *> = N+1: D is nonsingular, but RCOND is less than machine
262 *> precision, meaning that the matrix is singular
263 *> to working precision. Nevertheless, the
264 *> solution and error bounds are computed because
265 *> there are a number of situations where the
266 *> computed solution can be more accurate than the
267 *> value of RCOND would suggest.
268 *> \endverbatim
269 *
270 * Authors:
271 * ========
272 *
273 *> \author Univ. of Tennessee
274 *> \author Univ. of California Berkeley
275 *> \author Univ. of Colorado Denver
276 *> \author NAG Ltd.
277 *
278 *> \ingroup doubleSYsolve
279 *
280 * =====================================================================
281  SUBROUTINE dsysvx( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
282  $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
283  $ IWORK, INFO )
284 *
285 * -- LAPACK driver routine --
286 * -- LAPACK is a software package provided by Univ. of Tennessee, --
287 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
288 *
289 * .. Scalar Arguments ..
290  CHARACTER FACT, UPLO
291  INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
292  DOUBLE PRECISION RCOND
293 * ..
294 * .. Array Arguments ..
295  INTEGER IPIV( * ), IWORK( * )
296  DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
297  $ berr( * ), ferr( * ), work( * ), x( ldx, * )
298 * ..
299 *
300 * =====================================================================
301 *
302 * .. Parameters ..
303  DOUBLE PRECISION ZERO
304  PARAMETER ( ZERO = 0.0d+0 )
305 * ..
306 * .. Local Scalars ..
307  LOGICAL LQUERY, NOFACT
308  INTEGER LWKOPT, NB
309  DOUBLE PRECISION ANORM
310 * ..
311 * .. External Functions ..
312  LOGICAL LSAME
313  INTEGER ILAENV
314  DOUBLE PRECISION DLAMCH, DLANSY
315  EXTERNAL lsame, ilaenv, dlamch, dlansy
316 * ..
317 * .. External Subroutines ..
318  EXTERNAL dlacpy, dsycon, dsyrfs, dsytrf, dsytrs, xerbla
319 * ..
320 * .. Intrinsic Functions ..
321  INTRINSIC max
322 * ..
323 * .. Executable Statements ..
324 *
325 * Test the input parameters.
326 *
327  info = 0
328  nofact = lsame( fact, 'N' )
329  lquery = ( lwork.EQ.-1 )
330  IF( .NOT.nofact .AND. .NOT.lsame( fact, 'F' ) ) THEN
331  info = -1
332  ELSE IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) )
333  $ THEN
334  info = -2
335  ELSE IF( n.LT.0 ) THEN
336  info = -3
337  ELSE IF( nrhs.LT.0 ) THEN
338  info = -4
339  ELSE IF( lda.LT.max( 1, n ) ) THEN
340  info = -6
341  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
342  info = -8
343  ELSE IF( ldb.LT.max( 1, n ) ) THEN
344  info = -11
345  ELSE IF( ldx.LT.max( 1, n ) ) THEN
346  info = -13
347  ELSE IF( lwork.LT.max( 1, 3*n ) .AND. .NOT.lquery ) THEN
348  info = -18
349  END IF
350 *
351  IF( info.EQ.0 ) THEN
352  lwkopt = max( 1, 3*n )
353  IF( nofact ) THEN
354  nb = ilaenv( 1, 'DSYTRF', uplo, n, -1, -1, -1 )
355  lwkopt = max( lwkopt, n*nb )
356  END IF
357  work( 1 ) = lwkopt
358  END IF
359 *
360  IF( info.NE.0 ) THEN
361  CALL xerbla( 'DSYSVX', -info )
362  RETURN
363  ELSE IF( lquery ) THEN
364  RETURN
365  END IF
366 *
367  IF( nofact ) THEN
368 *
369 * Compute the factorization A = U*D*U**T or A = L*D*L**T.
370 *
371  CALL dlacpy( uplo, n, n, a, lda, af, ldaf )
372  CALL dsytrf( uplo, n, af, ldaf, ipiv, work, lwork, info )
373 *
374 * Return if INFO is non-zero.
375 *
376  IF( info.GT.0 )THEN
377  rcond = zero
378  RETURN
379  END IF
380  END IF
381 *
382 * Compute the norm of the matrix A.
383 *
384  anorm = dlansy( 'I', uplo, n, a, lda, work )
385 *
386 * Compute the reciprocal of the condition number of A.
387 *
388  CALL dsycon( uplo, n, af, ldaf, ipiv, anorm, rcond, work, iwork,
389  $ info )
390 *
391 * Compute the solution vectors X.
392 *
393  CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
394  CALL dsytrs( uplo, n, nrhs, af, ldaf, ipiv, x, ldx, info )
395 *
396 * Use iterative refinement to improve the computed solutions and
397 * compute error bounds and backward error estimates for them.
398 *
399  CALL dsyrfs( uplo, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,
400  $ ldx, ferr, berr, work, iwork, info )
401 *
402 * Set INFO = N+1 if the matrix is singular to working precision.
403 *
404  IF( rcond.LT.dlamch( 'Epsilon' ) )
405  $ info = n + 1
406 *
407  work( 1 ) = lwkopt
408 *
409  RETURN
410 *
411 * End of DSYSVX
412 *
413  END
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dsycon(UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, IWORK, INFO)
DSYCON
Definition: dsycon.f:130
subroutine dsytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
DSYTRS
Definition: dsytrs.f:120
subroutine dsyrfs(UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
DSYRFS
Definition: dsyrfs.f:191
subroutine dsytrf(UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
DSYTRF
Definition: dsytrf.f:182
subroutine dsysvx(FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, IWORK, INFO)
DSYSVX computes the solution to system of linear equations A * X = B for SY matrices
Definition: dsysvx.f:284