1 *> \brief <b> SSPSVX computes the solution to system of linear equations A * X = B for OTHER matrices</b>
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SSPSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sspsvx.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sspsvx.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sspsvx.f">
21 * SUBROUTINE SSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X,
22 * LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
24 * .. Scalar Arguments ..
25 * CHARACTER FACT, UPLO
26 * INTEGER INFO, LDB, LDX, N, NRHS
29 * .. Array Arguments ..
30 * INTEGER IPIV( * ), IWORK( * )
31 * REAL AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
32 * $ FERR( * ), WORK( * ), X( LDX, * )
41 *> SSPSVX 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.
46 *> Error bounds on the solution and a condition estimate are also
55 *> The following steps are performed:
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.
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.
71 *> 3. The system of equations is solved for X using the factored form
74 *> 4. Iterative refinement is applied to improve the computed solution
75 *> matrix and calculate error bounds and backward error estimates
84 *> FACT is CHARACTER*1
85 *> Specifies whether or not the factored form of A has been
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.
94 *> UPLO is CHARACTER*1
95 *> = 'U': Upper triangle of A is stored;
96 *> = 'L': Lower triangle of A is stored.
102 *> The number of linear equations, i.e., the order of the
109 *> The number of right hand sides, i.e., the number of columns
110 *> of the matrices B and X. NRHS >= 0.
115 *> AP is REAL 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.
124 *> \param[in,out] AFP
126 *> AFP is REAL 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 SSPTRF, stored as
131 *> a packed triangular matrix in the same storage format as A.
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 SSPTRF, stored as
137 *> a packed triangular matrix in the same storage format as A.
140 *> \param[in,out] IPIV
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 SSPTRF.
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.
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 SSPTRF.
161 *> B is REAL array, dimension (LDB,NRHS)
162 *> The N-by-NRHS right hand side matrix B.
168 *> The leading dimension of the array B. LDB >= max(1,N).
173 *> X is REAL array, dimension (LDX,NRHS)
174 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
180 *> The leading dimension of the array X. LDX >= max(1,N).
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
195 *> FERR is REAL 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.
208 *> BERR is REAL 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).
216 *> WORK is REAL array, dimension (3*N)
221 *> IWORK is INTEGER array, dimension (N)
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.
246 *> \author Univ. of Tennessee
247 *> \author Univ. of California Berkeley
248 *> \author Univ. of Colorado Denver
253 *> \ingroup realOTHERsolve
255 *> \par Further Details:
256 * =====================
260 *> The packed storage scheme is illustrated by the following example
261 *> when N = 4, UPLO = 'U':
263 *> Two-dimensional storage of the symmetric matrix A:
267 *> a33 a34 (aij = aji)
270 *> Packed storage of the upper triangle of A:
272 *> AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
275 * =====================================================================
276 SUBROUTINE SSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X,
277 $ LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
279 * -- LAPACK driver routine (version 3.7.0) --
280 * -- LAPACK is a software package provided by Univ. of Tennessee, --
281 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
284 * .. Scalar Arguments ..
286 INTEGER INFO, LDB, LDX, N, NRHS
289 * .. Array Arguments ..
290 INTEGER IPIV( * ), IWORK( * )
291 REAL AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
292 $ FERR( * ), WORK( * ), X( LDX, * )
295 * =====================================================================
299 PARAMETER ( ZERO = 0.0E+0 )
301 * .. Local Scalars ..
305 * .. External Functions ..
308 EXTERNAL LSAME, SLAMCH, SLANSP
310 * .. External Subroutines ..
311 EXTERNAL SCOPY, SLACPY, SSPCON, SSPRFS, SSPTRF, SSPTRS,
314 * .. Intrinsic Functions ..
317 * .. Executable Statements ..
319 * Test the input parameters.
322 NOFACT = LSAME( FACT, 'N' )
323 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
325 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
328 ELSE IF( N.LT.0 ) THEN
330 ELSE IF( NRHS.LT.0 ) THEN
332 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
334 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
338 CALL XERBLA( 'SSPSVX', -INFO )
344 * Compute the factorization A = U*D*U**T or A = L*D*L**T.
346 CALL SCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
347 CALL SSPTRF( UPLO, N, AFP, IPIV, INFO )
349 * Return if INFO is non-zero.
357 * Compute the norm of the matrix A.
359 ANORM = SLANSP( 'I', UPLO, N, AP, WORK )
361 * Compute the reciprocal of the condition number of A.
363 CALL SSPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, IWORK, INFO )
365 * Compute the solution vectors X.
367 CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
368 CALL SSPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO )
370 * Use iterative refinement to improve the computed solutions and
371 * compute error bounds and backward error estimates for them.
373 CALL SSPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR,
374 $ BERR, WORK, IWORK, INFO )
376 * Set INFO = N+1 if the matrix is singular to working precision.
378 IF( RCOND.LT.SLAMCH( 'Epsilon' ) )