LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 hesvx
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 xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsycon(uplo, n, a, lda, ipiv, anorm, rcond, work, iwork, info)
DSYCON
Definition dsycon.f:130
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 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
subroutine dsytrf(uplo, n, a, lda, ipiv, work, lwork, info)
DSYTRF
Definition dsytrf.f:182
subroutine dsytrs(uplo, n, nrhs, a, lda, ipiv, b, ldb, info)
DSYTRS
Definition dsytrs.f:120
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