LAPACK  3.6.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 *> \date April 2012
279 *
280 *> \ingroup doubleSYsolve
281 *
282 * =====================================================================
283  SUBROUTINE dsysvx( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
284  $ ldb, x, ldx, rcond, ferr, berr, work, lwork,
285  $ iwork, info )
286 *
287 * -- LAPACK driver routine (version 3.4.1) --
288 * -- LAPACK is a software package provided by Univ. of Tennessee, --
289 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
290 * April 2012
291 *
292 * .. Scalar Arguments ..
293  CHARACTER FACT, UPLO
294  INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
295  DOUBLE PRECISION RCOND
296 * ..
297 * .. Array Arguments ..
298  INTEGER IPIV( * ), IWORK( * )
299  DOUBLE PRECISION A( lda, * ), AF( ldaf, * ), B( ldb, * ),
300  $ berr( * ), ferr( * ), work( * ), x( ldx, * )
301 * ..
302 *
303 * =====================================================================
304 *
305 * .. Parameters ..
306  DOUBLE PRECISION ZERO
307  parameter ( zero = 0.0d+0 )
308 * ..
309 * .. Local Scalars ..
310  LOGICAL LQUERY, NOFACT
311  INTEGER LWKOPT, NB
312  DOUBLE PRECISION ANORM
313 * ..
314 * .. External Functions ..
315  LOGICAL LSAME
316  INTEGER ILAENV
317  DOUBLE PRECISION DLAMCH, DLANSY
318  EXTERNAL lsame, ilaenv, dlamch, dlansy
319 * ..
320 * .. External Subroutines ..
321  EXTERNAL dlacpy, dsycon, dsyrfs, dsytrf, dsytrs, xerbla
322 * ..
323 * .. Intrinsic Functions ..
324  INTRINSIC max
325 * ..
326 * .. Executable Statements ..
327 *
328 * Test the input parameters.
329 *
330  info = 0
331  nofact = lsame( fact, 'N' )
332  lquery = ( lwork.EQ.-1 )
333  IF( .NOT.nofact .AND. .NOT.lsame( fact, 'F' ) ) THEN
334  info = -1
335  ELSE IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) )
336  $ THEN
337  info = -2
338  ELSE IF( n.LT.0 ) THEN
339  info = -3
340  ELSE IF( nrhs.LT.0 ) THEN
341  info = -4
342  ELSE IF( lda.LT.max( 1, n ) ) THEN
343  info = -6
344  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
345  info = -8
346  ELSE IF( ldb.LT.max( 1, n ) ) THEN
347  info = -11
348  ELSE IF( ldx.LT.max( 1, n ) ) THEN
349  info = -13
350  ELSE IF( lwork.LT.max( 1, 3*n ) .AND. .NOT.lquery ) THEN
351  info = -18
352  END IF
353 *
354  IF( info.EQ.0 ) THEN
355  lwkopt = max( 1, 3*n )
356  IF( nofact ) THEN
357  nb = ilaenv( 1, 'DSYTRF', uplo, n, -1, -1, -1 )
358  lwkopt = max( lwkopt, n*nb )
359  END IF
360  work( 1 ) = lwkopt
361  END IF
362 *
363  IF( info.NE.0 ) THEN
364  CALL xerbla( 'DSYSVX', -info )
365  RETURN
366  ELSE IF( lquery ) THEN
367  RETURN
368  END IF
369 *
370  IF( nofact ) THEN
371 *
372 * Compute the factorization A = U*D*U**T or A = L*D*L**T.
373 *
374  CALL dlacpy( uplo, n, n, a, lda, af, ldaf )
375  CALL dsytrf( uplo, n, af, ldaf, ipiv, work, lwork, info )
376 *
377 * Return if INFO is non-zero.
378 *
379  IF( info.GT.0 )THEN
380  rcond = zero
381  RETURN
382  END IF
383  END IF
384 *
385 * Compute the norm of the matrix A.
386 *
387  anorm = dlansy( 'I', uplo, n, a, lda, work )
388 *
389 * Compute the reciprocal of the condition number of A.
390 *
391  CALL dsycon( uplo, n, af, ldaf, ipiv, anorm, rcond, work, iwork,
392  $ info )
393 *
394 * Compute the solution vectors X.
395 *
396  CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
397  CALL dsytrs( uplo, n, nrhs, af, ldaf, ipiv, x, ldx, info )
398 *
399 * Use iterative refinement to improve the computed solutions and
400 * compute error bounds and backward error estimates for them.
401 *
402  CALL dsyrfs( uplo, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,
403  $ ldx, ferr, berr, work, iwork, info )
404 *
405 * Set INFO = N+1 if the matrix is singular to working precision.
406 *
407  IF( rcond.LT.dlamch( 'Epsilon' ) )
408  $ info = n + 1
409 *
410  work( 1 ) = lwkopt
411 *
412  RETURN
413 *
414 * End of DSYSVX
415 *
416  END
subroutine dsyrfs(UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
DSYRFS
Definition: dsyrfs.f:193
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dsycon(UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, IWORK, INFO)
DSYCON
Definition: dsycon.f:132
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dsytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
DSYTRS
Definition: dsytrs.f:122
subroutine dsytrf(UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
DSYTRF
Definition: dsytrf.f:184
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:286