OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
real Other Solve Routines
LAPACK » Other Solve Routines

Functions

subroutine sgglse (m, n, p, a, lda, b, ldb, c, d, x, work, lwork, info)
  SGGLSE solves overdetermined or underdetermined systems for OTHER matrices
subroutine spbsv (uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
  SPBSV computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine spbsvx (fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info)
  SPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine sppsv (uplo, n, nrhs, ap, b, ldb, info)
  SPPSV computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine sppsvx (fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info)
  SPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine sspsv (uplo, n, nrhs, ap, ipiv, b, ldb, info)
  SSPSV computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine sspsvx (fact, uplo, n, nrhs, ap, afp, ipiv, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info)
  SSPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

Detailed Description

This is the group of real Other Solve routines

Function Documentation

◆ sgglse()

subroutine sgglse ( integer m,
integer n,
integer p,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) c,
real, dimension( * ) d,
real, dimension( * ) x,
real, dimension( * ) work,
integer lwork,
integer info )

SGGLSE solves overdetermined or underdetermined systems for OTHER matrices

Download SGGLSE + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SGGLSE solves the linear equality-constrained least squares (LSE)
!> problem:
!>
!>         minimize || c - A*x ||_2   subject to   B*x = d
!>
!> where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
!> M-vector, and d is a given P-vector. It is assumed that
!> P <= N <= M+P, and
!>
!>          rank(B) = P and  rank( (A) ) = N.
!>                               ( (B) )
!>
!> These conditions ensure that the LSE problem has a unique solution,
!> which is obtained using a generalized RQ factorization of the
!> matrices (B, A) given by
!>
!>    B = (0 R)*Q,   A = Z*T*Q.
!>
Parameters
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrices A and B. N >= 0.
!>
[in]P
!>          P is INTEGER
!>          The number of rows of the matrix B. 0 <= P <= N <= M+P.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit, the elements on and above the diagonal of the array
!>          contain the min(M,N)-by-N upper trapezoidal matrix T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,M).
!>
[in,out]B
!>          B is REAL array, dimension (LDB,N)
!>          On entry, the P-by-N matrix B.
!>          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
!>          contains the P-by-P upper triangular matrix R.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,P).
!>
[in,out]C
!>          C is REAL array, dimension (M)
!>          On entry, C contains the right hand side vector for the
!>          least squares part of the LSE problem.
!>          On exit, the residual sum of squares for the solution
!>          is given by the sum of squares of elements N-P+1 to M of
!>          vector C.
!>
[in,out]D
!>          D is REAL array, dimension (P)
!>          On entry, D contains the right hand side vector for the
!>          constrained equation.
!>          On exit, D is destroyed.
!>
[out]X
!>          X is REAL array, dimension (N)
!>          On exit, X is the solution of the LSE problem.
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= max(1,M+N+P).
!>          For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
!>          where NB is an upper bound for the optimal blocksizes for
!>          SGEQRF, SGERQF, SORMQR and SORMRQ.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1:  the upper triangular factor R associated with B in the
!>                generalized RQ factorization of the pair (B, A) is
!>                singular, so that rank(B) < P; the least squares
!>                solution could not be computed.
!>          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor
!>                T associated with A in the generalized RQ factorization
!>                of the pair (B, A) is singular, so that
!>                rank( (A) ) < N; the least squares solution could not
!>                    ( (B) )
!>                be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 178 of file sgglse.f.

180*
181* -- LAPACK driver routine --
182* -- LAPACK is a software package provided by Univ. of Tennessee, --
183* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184*
185* .. Scalar Arguments ..
186 INTEGER INFO, LDA, LDB, LWORK, M, N, P
187* ..
188* .. Array Arguments ..
189 REAL A( LDA, * ), B( LDB, * ), C( * ), D( * ),
190 $ WORK( * ), X( * )
191* ..
192*
193* =====================================================================
194*
195* .. Parameters ..
196 REAL ONE
197 parameter( one = 1.0e+0 )
198* ..
199* .. Local Scalars ..
200 LOGICAL LQUERY
201 INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
202 $ NB4, NR
203* ..
204* .. External Subroutines ..
205 EXTERNAL saxpy, scopy, sgemv, sggrqf, sormqr, sormrq,
206 $ strmv, strtrs, xerbla
207* ..
208* .. External Functions ..
209 INTEGER ILAENV
210 EXTERNAL ilaenv
211* ..
212* .. Intrinsic Functions ..
213 INTRINSIC int, max, min
214* ..
215* .. Executable Statements ..
216*
217* Test the input parameters
218*
219 info = 0
220 mn = min( m, n )
221 lquery = ( lwork.EQ.-1 )
222 IF( m.LT.0 ) THEN
223 info = -1
224 ELSE IF( n.LT.0 ) THEN
225 info = -2
226 ELSE IF( p.LT.0 .OR. p.GT.n .OR. p.LT.n-m ) THEN
227 info = -3
228 ELSE IF( lda.LT.max( 1, m ) ) THEN
229 info = -5
230 ELSE IF( ldb.LT.max( 1, p ) ) THEN
231 info = -7
232 END IF
233*
234* Calculate workspace
235*
236 IF( info.EQ.0) THEN
237 IF( n.EQ.0 ) THEN
238 lwkmin = 1
239 lwkopt = 1
240 ELSE
241 nb1 = ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
242 nb2 = ilaenv( 1, 'SGERQF', ' ', m, n, -1, -1 )
243 nb3 = ilaenv( 1, 'SORMQR', ' ', m, n, p, -1 )
244 nb4 = ilaenv( 1, 'SORMRQ', ' ', m, n, p, -1 )
245 nb = max( nb1, nb2, nb3, nb4 )
246 lwkmin = m + n + p
247 lwkopt = p + mn + max( m, n )*nb
248 END IF
249 work( 1 ) = lwkopt
250*
251 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
252 info = -12
253 END IF
254 END IF
255*
256 IF( info.NE.0 ) THEN
257 CALL xerbla( 'SGGLSE', -info )
258 RETURN
259 ELSE IF( lquery ) THEN
260 RETURN
261 END IF
262*
263* Quick return if possible
264*
265 IF( n.EQ.0 )
266 $ RETURN
267*
268* Compute the GRQ factorization of matrices B and A:
269*
270* B*Q**T = ( 0 T12 ) P Z**T*A*Q**T = ( R11 R12 ) N-P
271* N-P P ( 0 R22 ) M+P-N
272* N-P P
273*
274* where T12 and R11 are upper triangular, and Q and Z are
275* orthogonal.
276*
277 CALL sggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),
278 $ work( p+mn+1 ), lwork-p-mn, info )
279 lopt = work( p+mn+1 )
280*
281* Update c = Z**T *c = ( c1 ) N-P
282* ( c2 ) M+P-N
283*
284 CALL sormqr( 'Left', 'Transpose', m, 1, mn, a, lda, work( p+1 ),
285 $ c, max( 1, m ), work( p+mn+1 ), lwork-p-mn, info )
286 lopt = max( lopt, int( work( p+mn+1 ) ) )
287*
288* Solve T12*x2 = d for x2
289*
290 IF( p.GT.0 ) THEN
291 CALL strtrs( 'Upper', 'No transpose', 'Non-unit', p, 1,
292 $ b( 1, n-p+1 ), ldb, d, p, info )
293*
294 IF( info.GT.0 ) THEN
295 info = 1
296 RETURN
297 END IF
298*
299* Put the solution in X
300*
301 CALL scopy( p, d, 1, x( n-p+1 ), 1 )
302*
303* Update c1
304*
305 CALL sgemv( 'No transpose', n-p, p, -one, a( 1, n-p+1 ), lda,
306 $ d, 1, one, c, 1 )
307 END IF
308*
309* Solve R11*x1 = c1 for x1
310*
311 IF( n.GT.p ) THEN
312 CALL strtrs( 'Upper', 'No transpose', 'Non-unit', n-p, 1,
313 $ a, lda, c, n-p, info )
314*
315 IF( info.GT.0 ) THEN
316 info = 2
317 RETURN
318 END IF
319*
320* Put the solutions in X
321*
322 CALL scopy( n-p, c, 1, x, 1 )
323 END IF
324*
325* Compute the residual vector:
326*
327 IF( m.LT.n ) THEN
328 nr = m + p - n
329 IF( nr.GT.0 )
330 $ CALL sgemv( 'No transpose', nr, n-m, -one, a( n-p+1, m+1 ),
331 $ lda, d( nr+1 ), 1, one, c( n-p+1 ), 1 )
332 ELSE
333 nr = p
334 END IF
335 IF( nr.GT.0 ) THEN
336 CALL strmv( 'Upper', 'No transpose', 'Non unit', nr,
337 $ a( n-p+1, n-p+1 ), lda, d, 1 )
338 CALL saxpy( nr, -one, d, 1, c( n-p+1 ), 1 )
339 END IF
340*
341* Backward transformation x = Q**T*x
342*
343 CALL sormrq( 'Left', 'Transpose', n, 1, p, b, ldb, work( 1 ), x,
344 $ n, work( p+mn+1 ), lwork-p-mn, info )
345 work( 1 ) = p + mn + max( lopt, int( work( p+mn+1 ) ) )
346*
347 RETURN
348*
349* End of SGGLSE
350*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
strtrs
subroutine strtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
STRTRS
Definition strtrs.f:140
sormrq
subroutine sormrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMRQ
Definition sormrq.f:168
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
sggrqf
subroutine sggrqf(m, p, n, a, lda, taua, b, ldb, taub, work, lwork, info)
SGGRQF
Definition sggrqf.f:214
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:156
strmv
subroutine strmv(uplo, trans, diag, n, a, lda, x, incx)
STRMV
Definition strmv.f:147
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ spbsv()

subroutine spbsv ( character uplo,
integer n,
integer kd,
integer nrhs,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( ldb, * ) b,
integer ldb,
integer info )

SPBSV computes the solution to system of linear equations A * X = B for OTHER matrices

Download SPBSV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SPBSV computes the solution to a real system of linear equations
!>    A * X = B,
!> where A is an N-by-N symmetric positive definite band matrix and X
!> and B are N-by-NRHS matrices.
!>
!> The Cholesky decomposition is used to factor A as
!>    A = U**T * U,  if UPLO = 'U', or
!>    A = L * L**T,  if UPLO = 'L',
!> where U is an upper triangular band matrix, and L is a lower
!> triangular band matrix, with the same number of superdiagonals or
!> subdiagonals as A.  The factored form of A is then used to solve the
!> system of equations A * X = B.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]AB
!>          AB is REAL array, dimension (LDAB,N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(N,j+KD).
!>          See below for further details.
!>
!>          On exit, if INFO = 0, the triangular factor U or L from the
!>          Cholesky factorization A = U**T*U or A = L*L**T of the band
!>          matrix A, in the same storage format as A.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD+1.
!>
[in,out]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the leading minor of order i of A is not
!>                positive definite, so the factorization could not be
!>                completed, and the solution has not been computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The band storage scheme is illustrated by the following example, when
!>  N = 6, KD = 2, and UPLO = 'U':
!>
!>  On entry:                       On exit:
!>
!>      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
!>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
!>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
!>
!>  Similarly, if UPLO = 'L' the format of A is as follows:
!>
!>  On entry:                       On exit:
!>
!>     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
!>     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
!>     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
!>
!>  Array elements marked * are not used by the routine.
!>

Definition at line 163 of file spbsv.f.

164*
165* -- LAPACK driver routine --
166* -- LAPACK is a software package provided by Univ. of Tennessee, --
167* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168*
169* .. Scalar Arguments ..
170 CHARACTER UPLO
171 INTEGER INFO, KD, LDAB, LDB, N, NRHS
172* ..
173* .. Array Arguments ..
174 REAL AB( LDAB, * ), B( LDB, * )
175* ..
176*
177* =====================================================================
178*
179* .. External Functions ..
180 LOGICAL LSAME
181 EXTERNAL lsame
182* ..
183* .. External Subroutines ..
184 EXTERNAL spbtrf, spbtrs, xerbla
185* ..
186* .. Intrinsic Functions ..
187 INTRINSIC max
188* ..
189* .. Executable Statements ..
190*
191* Test the input parameters.
192*
193 info = 0
194 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
195 info = -1
196 ELSE IF( n.LT.0 ) THEN
197 info = -2
198 ELSE IF( kd.LT.0 ) THEN
199 info = -3
200 ELSE IF( nrhs.LT.0 ) THEN
201 info = -4
202 ELSE IF( ldab.LT.kd+1 ) THEN
203 info = -6
204 ELSE IF( ldb.LT.max( 1, n ) ) THEN
205 info = -8
206 END IF
207 IF( info.NE.0 ) THEN
208 CALL xerbla( 'SPBSV ', -info )
209 RETURN
210 END IF
211*
212* Compute the Cholesky factorization A = U**T*U or A = L*L**T.
213*
214 CALL spbtrf( uplo, n, kd, ab, ldab, info )
215 IF( info.EQ.0 ) THEN
216*
217* Solve the system A*X = B, overwriting B with X.
218*
219 CALL spbtrs( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
220*
221 END IF
222 RETURN
223*
224* End of SPBSV
225*
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
spbtrs
subroutine spbtrs(uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
SPBTRS
Definition spbtrs.f:121
spbtrf
subroutine spbtrf(uplo, n, kd, ab, ldab, info)
SPBTRF
Definition spbtrf.f:142

◆ spbsvx()

subroutine spbsvx ( character fact,
character uplo,
integer n,
integer kd,
integer nrhs,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( ldafb, * ) afb,
integer ldafb,
character equed,
real, dimension( * ) s,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldx, * ) x,
integer ldx,
real rcond,
real, dimension( * ) ferr,
real, dimension( * ) berr,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices

Download SPBSVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SPBSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
!> compute the solution to a real system of linear equations
!>    A * X = B,
!> where A is an N-by-N symmetric positive definite band matrix and X
!> and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'E', real scaling factors are computed to equilibrate
!>    the system:
!>       diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
!>    Whether or not the system will be equilibrated depends on the
!>    scaling of the matrix A, but if equilibration is used, A is
!>    overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
!>
!> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
!>    factor the matrix A (after equilibration if FACT = 'E') as
!>       A = U**T * U,  if UPLO = 'U', or
!>       A = L * L**T,  if UPLO = 'L',
!>    where U is an upper triangular band matrix, and L is a lower
!>    triangular band matrix.
!>
!> 3. If the leading i-by-i principal minor is not positive definite,
!>    then the routine returns with INFO = i. Otherwise, the factored
!>    form of A is used to estimate the condition number of the matrix
!>    A.  If the reciprocal of the condition number is less than machine
!>    precision, INFO = N+1 is returned as a warning, but the routine
!>    still goes on to solve for X and compute error bounds as
!>    described below.
!>
!> 4. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 5. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
!> 6. If equilibration was used, the matrix X is premultiplied by
!>    diag(S) so that it solves the original system before
!>    equilibration.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of the matrix A is
!>          supplied on entry, and if not, whether the matrix A should be
!>          equilibrated before it is factored.
!>          = 'F':  On entry, AFB contains the factored form of A.
!>                  If EQUED = 'Y', the matrix A has been equilibrated
!>                  with scaling factors given by S.  AB and AFB will not
!>                  be modified.
!>          = 'N':  The matrix A will be copied to AFB and factored.
!>          = 'E':  The matrix A will be equilibrated if necessary, then
!>                  copied to AFB and factored.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right-hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in,out]AB
!>          AB is REAL array, dimension (LDAB,N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix A, stored in the first KD+1 rows of the array, except
!>          if FACT = 'F' and EQUED = 'Y', then A must contain the
!>          equilibrated matrix diag(S)*A*diag(S).  The j-th column of A
!>          is stored in the j-th column of the array AB as follows:
!>          if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(N,j+KD).
!>          See below for further details.
!>
!>          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
!>          diag(S)*A*diag(S).
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array A.  LDAB >= KD+1.
!>
[in,out]AFB
!>          AFB is REAL array, dimension (LDAFB,N)
!>          If FACT = 'F', then AFB is an input argument and on entry
!>          contains the triangular factor U or L from the Cholesky
!>          factorization A = U**T*U or A = L*L**T of the band matrix
!>          A, in the same storage format as A (see AB).  If EQUED = 'Y',
!>          then AFB is the factored form of the equilibrated matrix A.
!>
!>          If FACT = 'N', then AFB is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**T*U or A = L*L**T.
!>
!>          If FACT = 'E', then AFB is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**T*U or A = L*L**T of the equilibrated
!>          matrix A (see the description of A for the form of the
!>          equilibrated matrix).
!>
[in]LDAFB
!>          LDAFB is INTEGER
!>          The leading dimension of the array AFB.  LDAFB >= KD+1.
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>          Specifies the form of equilibration that was done.
!>          = 'N':  No equilibration (always true if FACT = 'N').
!>          = 'Y':  Equilibration was done, i.e., A has been replaced by
!>                  diag(S) * A * diag(S).
!>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>          output argument.
!>
[in,out]S
!>          S is REAL array, dimension (N)
!>          The scale factors for A; not accessed if EQUED = 'N'.  S is
!>          an input argument if FACT = 'F'; otherwise, S is an output
!>          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
!>          must be positive.
!>
[in,out]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
!>          B is overwritten by diag(S) * B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is REAL array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
!>          the original system of equations.  Note that if EQUED = 'Y',
!>          A and B are modified on exit, and the solution to the
!>          equilibrated system is inv(diag(S))*X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is REAL
!>          The estimate of the reciprocal condition number of the matrix
!>          A after equilibration (if done).  If RCOND is less than the
!>          machine precision (in particular, if RCOND = 0), the matrix
!>          is singular to working precision.  This condition is
!>          indicated by a return code of INFO > 0.
!>
[out]FERR
!>          FERR is REAL array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is REAL array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is REAL array, dimension (3*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= N:  the leading minor of order i of A is
!>                       not positive definite, so the factorization
!>                       could not be completed, and the solution has not
!>                       been computed. RCOND = 0 is returned.
!>                = N+1: U is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The band storage scheme is illustrated by the following example, when
!>  N = 6, KD = 2, and UPLO = 'U':
!>
!>  Two-dimensional storage of the symmetric matrix A:
!>
!>     a11  a12  a13
!>          a22  a23  a24
!>               a33  a34  a35
!>                    a44  a45  a46
!>                         a55  a56
!>     (aij=conjg(aji))         a66
!>
!>  Band storage of the upper triangle of A:
!>
!>      *    *   a13  a24  a35  a46
!>      *   a12  a23  a34  a45  a56
!>     a11  a22  a33  a44  a55  a66
!>
!>  Similarly, if UPLO = 'L' the format of A is as follows:
!>
!>     a11  a22  a33  a44  a55  a66
!>     a21  a32  a43  a54  a65   *
!>     a31  a42  a53  a64   *    *
!>
!>  Array elements marked * are not used by the routine.
!>

Definition at line 340 of file spbsvx.f.

343*
344* -- LAPACK driver routine --
345* -- LAPACK is a software package provided by Univ. of Tennessee, --
346* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
347*
348* .. Scalar Arguments ..
349 CHARACTER EQUED, FACT, UPLO
350 INTEGER INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
351 REAL RCOND
352* ..
353* .. Array Arguments ..
354 INTEGER IWORK( * )
355 REAL AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
356 $ BERR( * ), FERR( * ), S( * ), WORK( * ),
357 $ X( LDX, * )
358* ..
359*
360* =====================================================================
361*
362* .. Parameters ..
363 REAL ZERO, ONE
364 parameter( zero = 0.0e+0, one = 1.0e+0 )
365* ..
366* .. Local Scalars ..
367 LOGICAL EQUIL, NOFACT, RCEQU, UPPER
368 INTEGER I, INFEQU, J, J1, J2
369 REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
370* ..
371* .. External Functions ..
372 LOGICAL LSAME
373 REAL SLAMCH, SLANSB
374 EXTERNAL lsame, slamch, slansb
375* ..
376* .. External Subroutines ..
377 EXTERNAL scopy, slacpy, slaqsb, spbcon, spbequ, spbrfs,
378 $ spbtrf, spbtrs, xerbla
379* ..
380* .. Intrinsic Functions ..
381 INTRINSIC max, min
382* ..
383* .. Executable Statements ..
384*
385 info = 0
386 nofact = lsame( fact, 'N' )
387 equil = lsame( fact, 'E' )
388 upper = lsame( uplo, 'U' )
389 IF( nofact .OR. equil ) THEN
390 equed = 'N'
391 rcequ = .false.
392 ELSE
393 rcequ = lsame( equed, 'Y' )
394 smlnum = slamch( 'Safe minimum' )
395 bignum = one / smlnum
396 END IF
397*
398* Test the input parameters.
399*
400 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
401 $ THEN
402 info = -1
403 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
404 info = -2
405 ELSE IF( n.LT.0 ) THEN
406 info = -3
407 ELSE IF( kd.LT.0 ) THEN
408 info = -4
409 ELSE IF( nrhs.LT.0 ) THEN
410 info = -5
411 ELSE IF( ldab.LT.kd+1 ) THEN
412 info = -7
413 ELSE IF( ldafb.LT.kd+1 ) THEN
414 info = -9
415 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
416 $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
417 info = -10
418 ELSE
419 IF( rcequ ) THEN
420 smin = bignum
421 smax = zero
422 DO 10 j = 1, n
423 smin = min( smin, s( j ) )
424 smax = max( smax, s( j ) )
425 10 CONTINUE
426 IF( smin.LE.zero ) THEN
427 info = -11
428 ELSE IF( n.GT.0 ) THEN
429 scond = max( smin, smlnum ) / min( smax, bignum )
430 ELSE
431 scond = one
432 END IF
433 END IF
434 IF( info.EQ.0 ) THEN
435 IF( ldb.LT.max( 1, n ) ) THEN
436 info = -13
437 ELSE IF( ldx.LT.max( 1, n ) ) THEN
438 info = -15
439 END IF
440 END IF
441 END IF
442*
443 IF( info.NE.0 ) THEN
444 CALL xerbla( 'SPBSVX', -info )
445 RETURN
446 END IF
447*
448 IF( equil ) THEN
449*
450* Compute row and column scalings to equilibrate the matrix A.
451*
452 CALL spbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
453 IF( infequ.EQ.0 ) THEN
454*
455* Equilibrate the matrix.
456*
457 CALL slaqsb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
458 rcequ = lsame( equed, 'Y' )
459 END IF
460 END IF
461*
462* Scale the right-hand side.
463*
464 IF( rcequ ) THEN
465 DO 30 j = 1, nrhs
466 DO 20 i = 1, n
467 b( i, j ) = s( i )*b( i, j )
468 20 CONTINUE
469 30 CONTINUE
470 END IF
471*
472 IF( nofact .OR. equil ) THEN
473*
474* Compute the Cholesky factorization A = U**T *U or A = L*L**T.
475*
476 IF( upper ) THEN
477 DO 40 j = 1, n
478 j1 = max( j-kd, 1 )
479 CALL scopy( j-j1+1, ab( kd+1-j+j1, j ), 1,
480 $ afb( kd+1-j+j1, j ), 1 )
481 40 CONTINUE
482 ELSE
483 DO 50 j = 1, n
484 j2 = min( j+kd, n )
485 CALL scopy( j2-j+1, ab( 1, j ), 1, afb( 1, j ), 1 )
486 50 CONTINUE
487 END IF
488*
489 CALL spbtrf( uplo, n, kd, afb, ldafb, info )
490*
491* Return if INFO is non-zero.
492*
493 IF( info.GT.0 )THEN
494 rcond = zero
495 RETURN
496 END IF
497 END IF
498*
499* Compute the norm of the matrix A.
500*
501 anorm = slansb( '1', uplo, n, kd, ab, ldab, work )
502*
503* Compute the reciprocal of the condition number of A.
504*
505 CALL spbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, iwork,
506 $ info )
507*
508* Compute the solution matrix X.
509*
510 CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldx )
511 CALL spbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
512*
513* Use iterative refinement to improve the computed solution and
514* compute error bounds and backward error estimates for it.
515*
516 CALL spbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,
517 $ ldx, ferr, berr, work, iwork, info )
518*
519* Transform the solution matrix X to a solution of the original
520* system.
521*
522 IF( rcequ ) THEN
523 DO 70 j = 1, nrhs
524 DO 60 i = 1, n
525 x( i, j ) = s( i )*x( i, j )
526 60 CONTINUE
527 70 CONTINUE
528 DO 80 j = 1, nrhs
529 ferr( j ) = ferr( j ) / scond
530 80 CONTINUE
531 END IF
532*
533* Set INFO = N+1 if the matrix is singular to working precision.
534*
535 IF( rcond.LT.slamch( 'Epsilon' ) )
536 $ info = n + 1
537*
538 RETURN
539*
540* End of SPBSVX
541*
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
slansb
real function slansb(norm, uplo, n, k, ab, ldab, work)
SLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansb.f:129
slaqsb
subroutine slaqsb(uplo, n, kd, ab, ldab, s, scond, amax, equed)
SLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ.
Definition slaqsb.f:140
spbequ
subroutine spbequ(uplo, n, kd, ab, ldab, s, scond, amax, info)
SPBEQU
Definition spbequ.f:129
spbcon
subroutine spbcon(uplo, n, kd, ab, ldab, anorm, rcond, work, iwork, info)
SPBCON
Definition spbcon.f:132
spbrfs
subroutine spbrfs(uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x, ldx, ferr, berr, work, iwork, info)
SPBRFS
Definition spbrfs.f:189
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68

◆ sppsv()

subroutine sppsv ( character uplo,
integer n,
integer nrhs,
real, dimension( * ) ap,
real, dimension( ldb, * ) b,
integer ldb,
integer info )

SPPSV computes the solution to system of linear equations A * X = B for OTHER matrices

Download SPPSV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SPPSV computes the solution to a real system of linear equations
!>    A * X = B,
!> where A is an N-by-N symmetric positive definite matrix stored in
!> packed format and X and B are N-by-NRHS matrices.
!>
!> The Cholesky decomposition is used to factor A as
!>    A = U**T* U,  if UPLO = 'U', or
!>    A = L * L**T,  if UPLO = 'L',
!> where U is an upper triangular matrix and L is a lower triangular
!> matrix.  The factored form of A is then used to solve the system of
!> equations A * X = B.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.
!>
!>          On exit, if INFO = 0, the factor U or L from the Cholesky
!>          factorization A = U**T*U or A = L*L**T, in the same storage
!>          format as A.
!>
[in,out]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the leading minor of order i of A is not
!>                positive definite, so the factorization could not be
!>                completed, and the solution has not been computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the symmetric matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = conjg(aji))
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 143 of file sppsv.f.

144*
145* -- LAPACK driver routine --
146* -- LAPACK is a software package provided by Univ. of Tennessee, --
147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148*
149* .. Scalar Arguments ..
150 CHARACTER UPLO
151 INTEGER INFO, LDB, N, NRHS
152* ..
153* .. Array Arguments ..
154 REAL AP( * ), B( LDB, * )
155* ..
156*
157* =====================================================================
158*
159* .. External Functions ..
160 LOGICAL LSAME
161 EXTERNAL lsame
162* ..
163* .. External Subroutines ..
164 EXTERNAL spptrf, spptrs, xerbla
165* ..
166* .. Intrinsic Functions ..
167 INTRINSIC max
168* ..
169* .. Executable Statements ..
170*
171* Test the input parameters.
172*
173 info = 0
174 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
175 info = -1
176 ELSE IF( n.LT.0 ) THEN
177 info = -2
178 ELSE IF( nrhs.LT.0 ) THEN
179 info = -3
180 ELSE IF( ldb.LT.max( 1, n ) ) THEN
181 info = -6
182 END IF
183 IF( info.NE.0 ) THEN
184 CALL xerbla( 'SPPSV ', -info )
185 RETURN
186 END IF
187*
188* Compute the Cholesky factorization A = U**T*U or A = L*L**T.
189*
190 CALL spptrf( uplo, n, ap, info )
191 IF( info.EQ.0 ) THEN
192*
193* Solve the system A*X = B, overwriting B with X.
194*
195 CALL spptrs( uplo, n, nrhs, ap, b, ldb, info )
196*
197 END IF
198 RETURN
199*
200* End of SPPSV
201*
spptrf
subroutine spptrf(uplo, n, ap, info)
SPPTRF
Definition spptrf.f:119
spptrs
subroutine spptrs(uplo, n, nrhs, ap, b, ldb, info)
SPPTRS
Definition spptrs.f:108

◆ sppsvx()

subroutine sppsvx ( character fact,
character uplo,
integer n,
integer nrhs,
real, dimension( * ) ap,
real, dimension( * ) afp,
character equed,
real, dimension( * ) s,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldx, * ) x,
integer ldx,
real rcond,
real, dimension( * ) ferr,
real, dimension( * ) berr,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

Download SPPSVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
!> compute the solution to a real system of linear equations
!>    A * X = B,
!> where A is an N-by-N symmetric positive definite matrix stored in
!> packed format and X and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'E', real scaling factors are computed to equilibrate
!>    the system:
!>       diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
!>    Whether or not the system will be equilibrated depends on the
!>    scaling of the matrix A, but if equilibration is used, A is
!>    overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
!>
!> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
!>    factor the matrix A (after equilibration if FACT = 'E') as
!>       A = U**T* U,  if UPLO = 'U', or
!>       A = L * L**T,  if UPLO = 'L',
!>    where U is an upper triangular matrix and L is a lower triangular
!>    matrix.
!>
!> 3. If the leading i-by-i principal minor is not positive definite,
!>    then the routine returns with INFO = i. Otherwise, the factored
!>    form of A is used to estimate the condition number of the matrix
!>    A.  If the reciprocal of the condition number is less than machine
!>    precision, INFO = N+1 is returned as a warning, but the routine
!>    still goes on to solve for X and compute error bounds as
!>    described below.
!>
!> 4. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 5. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
!> 6. If equilibration was used, the matrix X is premultiplied by
!>    diag(S) so that it solves the original system before
!>    equilibration.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of the matrix A is
!>          supplied on entry, and if not, whether the matrix A should be
!>          equilibrated before it is factored.
!>          = 'F':  On entry, AFP contains the factored form of A.
!>                  If EQUED = 'Y', the matrix A has been equilibrated
!>                  with scaling factors given by S.  AP and AFP will not
!>                  be modified.
!>          = 'N':  The matrix A will be copied to AFP and factored.
!>          = 'E':  The matrix A will be equilibrated if necessary, then
!>                  copied to AFP and factored.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in,out]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array, except if FACT = 'F'
!>          and EQUED = 'Y', then A must contain the equilibrated matrix
!>          diag(S)*A*diag(S).  The j-th column of A is stored in the
!>          array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.  A is not modified if
!>          FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
!>
!>          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
!>          diag(S)*A*diag(S).
!>
[in,out]AFP
!>          AFP is REAL array, dimension (N*(N+1)/2)
!>          If FACT = 'F', then AFP is an input argument and on entry
!>          contains the triangular factor U or L from the Cholesky
!>          factorization A = U**T*U or A = L*L**T, in the same storage
!>          format as A.  If EQUED .ne. 'N', then AFP is the factored
!>          form of the equilibrated matrix A.
!>
!>          If FACT = 'N', then AFP is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**T * U or A = L * L**T of the original
!>          matrix A.
!>
!>          If FACT = 'E', then AFP is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**T * U or A = L * L**T of the equilibrated
!>          matrix A (see the description of AP for the form of the
!>          equilibrated matrix).
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>          Specifies the form of equilibration that was done.
!>          = 'N':  No equilibration (always true if FACT = 'N').
!>          = 'Y':  Equilibration was done, i.e., A has been replaced by
!>                  diag(S) * A * diag(S).
!>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>          output argument.
!>
[in,out]S
!>          S is REAL array, dimension (N)
!>          The scale factors for A; not accessed if EQUED = 'N'.  S is
!>          an input argument if FACT = 'F'; otherwise, S is an output
!>          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
!>          must be positive.
!>
[in,out]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
!>          B is overwritten by diag(S) * B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is REAL array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
!>          the original system of equations.  Note that if EQUED = 'Y',
!>          A and B are modified on exit, and the solution to the
!>          equilibrated system is inv(diag(S))*X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is REAL
!>          The estimate of the reciprocal condition number of the matrix
!>          A after equilibration (if done).  If RCOND is less than the
!>          machine precision (in particular, if RCOND = 0), the matrix
!>          is singular to working precision.  This condition is
!>          indicated by a return code of INFO > 0.
!>
[out]FERR
!>          FERR is REAL array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is REAL array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is REAL array, dimension (3*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= N:  the leading minor of order i of A is
!>                       not positive definite, so the factorization
!>                       could not be completed, and the solution has not
!>                       been computed. RCOND = 0 is returned.
!>                = N+1: U is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the symmetric matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = conjg(aji))
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 309 of file sppsvx.f.

311*
312* -- LAPACK driver routine --
313* -- LAPACK is a software package provided by Univ. of Tennessee, --
314* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
315*
316* .. Scalar Arguments ..
317 CHARACTER EQUED, FACT, UPLO
318 INTEGER INFO, LDB, LDX, N, NRHS
319 REAL RCOND
320* ..
321* .. Array Arguments ..
322 INTEGER IWORK( * )
323 REAL AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
324 $ FERR( * ), S( * ), WORK( * ), X( LDX, * )
325* ..
326*
327* =====================================================================
328*
329* .. Parameters ..
330 REAL ZERO, ONE
331 parameter( zero = 0.0e+0, one = 1.0e+0 )
332* ..
333* .. Local Scalars ..
334 LOGICAL EQUIL, NOFACT, RCEQU
335 INTEGER I, INFEQU, J
336 REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
337* ..
338* .. External Functions ..
339 LOGICAL LSAME
340 REAL SLAMCH, SLANSP
341 EXTERNAL lsame, slamch, slansp
342* ..
343* .. External Subroutines ..
344 EXTERNAL scopy, slacpy, slaqsp, sppcon, sppequ, spprfs,
345 $ spptrf, spptrs, xerbla
346* ..
347* .. Intrinsic Functions ..
348 INTRINSIC max, min
349* ..
350* .. Executable Statements ..
351*
352 info = 0
353 nofact = lsame( fact, 'N' )
354 equil = lsame( fact, 'E' )
355 IF( nofact .OR. equil ) THEN
356 equed = 'N'
357 rcequ = .false.
358 ELSE
359 rcequ = lsame( equed, 'Y' )
360 smlnum = slamch( 'Safe minimum' )
361 bignum = one / smlnum
362 END IF
363*
364* Test the input parameters.
365*
366 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
367 $ THEN
368 info = -1
369 ELSE IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) )
370 $ THEN
371 info = -2
372 ELSE IF( n.LT.0 ) THEN
373 info = -3
374 ELSE IF( nrhs.LT.0 ) THEN
375 info = -4
376 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
377 $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
378 info = -7
379 ELSE
380 IF( rcequ ) THEN
381 smin = bignum
382 smax = zero
383 DO 10 j = 1, n
384 smin = min( smin, s( j ) )
385 smax = max( smax, s( j ) )
386 10 CONTINUE
387 IF( smin.LE.zero ) THEN
388 info = -8
389 ELSE IF( n.GT.0 ) THEN
390 scond = max( smin, smlnum ) / min( smax, bignum )
391 ELSE
392 scond = one
393 END IF
394 END IF
395 IF( info.EQ.0 ) THEN
396 IF( ldb.LT.max( 1, n ) ) THEN
397 info = -10
398 ELSE IF( ldx.LT.max( 1, n ) ) THEN
399 info = -12
400 END IF
401 END IF
402 END IF
403*
404 IF( info.NE.0 ) THEN
405 CALL xerbla( 'SPPSVX', -info )
406 RETURN
407 END IF
408*
409 IF( equil ) THEN
410*
411* Compute row and column scalings to equilibrate the matrix A.
412*
413 CALL sppequ( uplo, n, ap, s, scond, amax, infequ )
414 IF( infequ.EQ.0 ) THEN
415*
416* Equilibrate the matrix.
417*
418 CALL slaqsp( uplo, n, ap, s, scond, amax, equed )
419 rcequ = lsame( equed, 'Y' )
420 END IF
421 END IF
422*
423* Scale the right-hand side.
424*
425 IF( rcequ ) THEN
426 DO 30 j = 1, nrhs
427 DO 20 i = 1, n
428 b( i, j ) = s( i )*b( i, j )
429 20 CONTINUE
430 30 CONTINUE
431 END IF
432*
433 IF( nofact .OR. equil ) THEN
434*
435* Compute the Cholesky factorization A = U**T * U or A = L * L**T.
436*
437 CALL scopy( n*( n+1 ) / 2, ap, 1, afp, 1 )
438 CALL spptrf( uplo, n, afp, info )
439*
440* Return if INFO is non-zero.
441*
442 IF( info.GT.0 )THEN
443 rcond = zero
444 RETURN
445 END IF
446 END IF
447*
448* Compute the norm of the matrix A.
449*
450 anorm = slansp( 'I', uplo, n, ap, work )
451*
452* Compute the reciprocal of the condition number of A.
453*
454 CALL sppcon( uplo, n, afp, anorm, rcond, work, iwork, info )
455*
456* Compute the solution matrix X.
457*
458 CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldx )
459 CALL spptrs( uplo, n, nrhs, afp, x, ldx, info )
460*
461* Use iterative refinement to improve the computed solution and
462* compute error bounds and backward error estimates for it.
463*
464 CALL spprfs( uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr,
465 $ work, iwork, info )
466*
467* Transform the solution matrix X to a solution of the original
468* system.
469*
470 IF( rcequ ) THEN
471 DO 50 j = 1, nrhs
472 DO 40 i = 1, n
473 x( i, j ) = s( i )*x( i, j )
474 40 CONTINUE
475 50 CONTINUE
476 DO 60 j = 1, nrhs
477 ferr( j ) = ferr( j ) / scond
478 60 CONTINUE
479 END IF
480*
481* Set INFO = N+1 if the matrix is singular to working precision.
482*
483 IF( rcond.LT.slamch( 'Epsilon' ) )
484 $ info = n + 1
485*
486 RETURN
487*
488* End of SPPSVX
489*
slaqsp
subroutine slaqsp(uplo, n, ap, s, scond, amax, equed)
SLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppeq...
Definition slaqsp.f:125
slansp
real function slansp(norm, uplo, n, ap, work)
SLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansp.f:114
sppequ
subroutine sppequ(uplo, n, ap, s, scond, amax, info)
SPPEQU
Definition sppequ.f:116
sppcon
subroutine sppcon(uplo, n, ap, anorm, rcond, work, iwork, info)
SPPCON
Definition sppcon.f:118
spprfs
subroutine spprfs(uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr, work, iwork, info)
SPPRFS
Definition spprfs.f:171

◆ sspsv()

subroutine sspsv ( character uplo,
integer n,
integer nrhs,
real, dimension( * ) ap,
integer, dimension( * ) ipiv,
real, dimension( ldb, * ) b,
integer ldb,
integer info )

SSPSV computes the solution to system of linear equations A * X = B for OTHER matrices

Download SSPSV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SSPSV computes the solution to a real system of linear equations
!>    A * X = B,
!> where A is an N-by-N symmetric matrix stored in packed format and X
!> and B are N-by-NRHS matrices.
!>
!> The diagonal pivoting method is used to factor A as
!>    A = U * D * U**T,  if UPLO = 'U', or
!>    A = L * D * L**T,  if UPLO = 'L',
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, D is symmetric and block diagonal with 1-by-1
!> and 2-by-2 diagonal blocks.  The factored form of A is then used to
!> solve the system of equations A * X = B.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.
!>
!>          On exit, the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          Details of the interchanges and the block structure of D, as
!>          determined by SSPTRF.  If IPIV(k) > 0, then rows and columns
!>          k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1
!>          diagonal block.  If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0,
!>          then rows and columns k-1 and -IPIV(k) were interchanged and
!>          D(k-1:k,k-1:k) is a 2-by-2 diagonal block.  If UPLO = 'L' and
!>          IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and
!>          -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2
!>          diagonal block.
!>
[in,out]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
!>                has been completed, but the block diagonal matrix D is
!>                exactly singular, so the solution could not be
!>                computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the symmetric matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = aji)
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 161 of file sspsv.f.

162*
163* -- LAPACK driver routine --
164* -- LAPACK is a software package provided by Univ. of Tennessee, --
165* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*
167* .. Scalar Arguments ..
168 CHARACTER UPLO
169 INTEGER INFO, LDB, N, NRHS
170* ..
171* .. Array Arguments ..
172 INTEGER IPIV( * )
173 REAL AP( * ), B( LDB, * )
174* ..
175*
176* =====================================================================
177*
178* .. External Functions ..
179 LOGICAL LSAME
180 EXTERNAL lsame
181* ..
182* .. External Subroutines ..
183 EXTERNAL ssptrf, ssptrs, xerbla
184* ..
185* .. Intrinsic Functions ..
186 INTRINSIC max
187* ..
188* .. Executable Statements ..
189*
190* Test the input parameters.
191*
192 info = 0
193 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
194 info = -1
195 ELSE IF( n.LT.0 ) THEN
196 info = -2
197 ELSE IF( nrhs.LT.0 ) THEN
198 info = -3
199 ELSE IF( ldb.LT.max( 1, n ) ) THEN
200 info = -7
201 END IF
202 IF( info.NE.0 ) THEN
203 CALL xerbla( 'SSPSV ', -info )
204 RETURN
205 END IF
206*
207* Compute the factorization A = U*D*U**T or A = L*D*L**T.
208*
209 CALL ssptrf( uplo, n, ap, ipiv, info )
210 IF( info.EQ.0 ) THEN
211*
212* Solve the system A*X = B, overwriting B with X.
213*
214 CALL ssptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
215*
216 END IF
217 RETURN
218*
219* End of SSPSV
220*
ssptrf
subroutine ssptrf(uplo, n, ap, ipiv, info)
SSPTRF
Definition ssptrf.f:157
ssptrs
subroutine ssptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
SSPTRS
Definition ssptrs.f:115

◆ sspsvx()

subroutine sspsvx ( character fact,
character uplo,
integer n,
integer nrhs,
real, dimension( * ) ap,
real, dimension( * ) afp,
integer, dimension( * ) ipiv,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldx, * ) x,
integer ldx,
real rcond,
real, dimension( * ) ferr,
real, dimension( * ) berr,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SSPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

Download SSPSVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SSPSVX uses the diagonal pivoting factorization A = U*D*U**T or
!> A = L*D*L**T to compute the solution to a real system of linear
!> equations A * X = B, where A is an N-by-N symmetric matrix stored
!> in packed format and X and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'N', the diagonal pivoting method is used to factor A as
!>       A = U * D * U**T,  if UPLO = 'U', or
!>       A = L * D * L**T,  if UPLO = 'L',
!>    where U (or L) is a product of permutation and unit upper (lower)
!>    triangular matrices and D is symmetric and block diagonal with
!>    1-by-1 and 2-by-2 diagonal blocks.
!>
!> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
!>    returns with INFO = i. Otherwise, the factored form of A is used
!>    to estimate the condition number of the matrix A.  If the
!>    reciprocal of the condition number is less than machine precision,
!>    INFO = N+1 is returned as a warning, but the routine still goes on
!>    to solve for X and compute error bounds as described below.
!>
!> 3. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 4. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of A has been
!>          supplied on entry.
!>          = 'F':  On entry, AFP and IPIV contain the factored form of
!>                  A.  AP, AFP and IPIV will not be modified.
!>          = 'N':  The matrix A will be copied to AFP and factored.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          The upper or lower triangle of the symmetric matrix A, packed
!>          columnwise in a linear array.  The j-th column of A is stored
!>          in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.
!>
[in,out]AFP
!>          AFP is REAL array, dimension (N*(N+1)/2)
!>          If FACT = 'F', then AFP is an input argument and on entry
!>          contains the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
!>          If FACT = 'N', then AFP is an output argument and on exit
!>          contains the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
[in,out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          If FACT = 'F', then IPIV is an input argument and on entry
!>          contains details of the interchanges and the block structure
!>          of D, as determined by SSPTRF.
!>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
!>          interchanged and D(k,k) is a 1-by-1 diagonal block.
!>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
!>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
!>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
!>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
!>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
!>
!>          If FACT = 'N', then IPIV is an output argument and on exit
!>          contains details of the interchanges and the block structure
!>          of D, as determined by SSPTRF.
!>
[in]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          The N-by-NRHS right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is REAL array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is REAL
!>          The estimate of the reciprocal condition number of the matrix
!>          A.  If RCOND is less than the machine precision (in
!>          particular, if RCOND = 0), the matrix is singular to working
!>          precision.  This condition is indicated by a return code of
!>          INFO > 0.
!>
[out]FERR
!>          FERR is REAL array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is REAL array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is REAL array, dimension (3*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= N:  D(i,i) is exactly zero.  The factorization
!>                       has been completed but the factor D is exactly
!>                       singular, so the solution and error bounds could
!>                       not be computed. RCOND = 0 is returned.
!>                = N+1: D is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the symmetric matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = aji)
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 274 of file sspsvx.f.

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 REAL RCOND
285* ..
286* .. Array Arguments ..
287 INTEGER IPIV( * ), IWORK( * )
288 REAL AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
289 $ FERR( * ), WORK( * ), X( LDX, * )
290* ..
291*
292* =====================================================================
293*
294* .. Parameters ..
295 REAL ZERO
296 parameter( zero = 0.0e+0 )
297* ..
298* .. Local Scalars ..
299 LOGICAL NOFACT
300 REAL ANORM
301* ..
302* .. External Functions ..
303 LOGICAL LSAME
304 REAL SLAMCH, SLANSP
305 EXTERNAL lsame, slamch, slansp
306* ..
307* .. External Subroutines ..
308 EXTERNAL scopy, slacpy, sspcon, ssprfs, ssptrf, ssptrs,
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( 'SSPSVX', -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 scopy( n*( n+1 ) / 2, ap, 1, afp, 1 )
344 CALL ssptrf( 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 = slansp( 'I', uplo, n, ap, work )
357*
358* Compute the reciprocal of the condition number of A.
359*
360 CALL sspcon( uplo, n, afp, ipiv, anorm, rcond, work, iwork, info )
361*
362* Compute the solution vectors X.
363*
364 CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldx )
365 CALL ssptrs( 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 ssprfs( 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.slamch( 'Epsilon' ) )
376 $ info = n + 1
377*
378 RETURN
379*
380* End of SSPSVX
381*
ssprfs
subroutine ssprfs(uplo, n, nrhs, ap, afp, ipiv, b, ldb, x, ldx, ferr, berr, work, iwork, info)
SSPRFS
Definition ssprfs.f:179
sspcon
subroutine sspcon(uplo, n, ap, ipiv, anorm, rcond, work, iwork, info)
SSPCON
Definition sspcon.f:125