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