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