LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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