OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
real
LAPACK » Positive Definite Matrix » Computational routines

Functions

real function sla_porcond (uplo, n, a, lda, af, ldaf, cmode, c, info, work, iwork)
 SLA_PORCOND estimates the Skeel condition number for a symmetric positive-definite matrix.
subroutine sla_porfsx_extended (prec_type, uplo, n, nrhs, a, lda, af, ldaf, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)
 SLA_PORFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric or Hermitian positive-definite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.
real function sla_porpvgrw (uplo, ncols, a, lda, af, ldaf, work)
 SLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian positive-definite matrix.
subroutine spocon (uplo, n, a, lda, anorm, rcond, work, iwork, info)
 SPOCON
subroutine spoequ (n, a, lda, s, scond, amax, info)
 SPOEQU
subroutine spoequb (n, a, lda, s, scond, amax, info)
 SPOEQUB
subroutine sporfs (uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx, ferr, berr, work, iwork, info)
 SPORFS
subroutine sporfsx (uplo, equed, n, nrhs, a, lda, af, ldaf, s, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
 SPORFSX
subroutine spotf2 (uplo, n, a, lda, info)
 SPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm).
subroutine spotrf (uplo, n, a, lda, info)
 SPOTRF
recursive subroutine spotrf2 (uplo, n, a, lda, info)
 SPOTRF2
subroutine spotri (uplo, n, a, lda, info)
 SPOTRI
subroutine spotrs (uplo, n, nrhs, a, lda, b, ldb, info)
 SPOTRS

Detailed Description

This is the group of real computational functions for PO matrices

Function Documentation

◆ sla_porcond()

real function sla_porcond ( character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldaf, * ) af,
integer ldaf,
integer cmode,
real, dimension( * ) c,
integer info,
real, dimension( * ) work,
integer, dimension( * ) iwork )

SLA_PORCOND estimates the Skeel condition number for a symmetric positive-definite matrix.

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

Purpose:
!>
!>    SLA_PORCOND Estimates the Skeel condition number of  op(A) * op2(C)
!>    where op2 is determined by CMODE as follows
!>    CMODE =  1    op2(C) = C
!>    CMODE =  0    op2(C) = I
!>    CMODE = -1    op2(C) = inv(C)
!>    The Skeel condition number  cond(A) = norminf( |inv(A)||A| )
!>    is computed by computing scaling factors R such that
!>    diag(R)*A*op2(C) is row equilibrated and computing the standard
!>    infinity-norm condition number.
!>
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]A
!>          A is REAL array, dimension (LDA,N)
!>     On entry, the N-by-N matrix A.
!>
[in]LDA
!>          LDA is INTEGER
!>     The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]AF
!>          AF is REAL array, dimension (LDAF,N)
!>     The triangular factor U or L from the Cholesky factorization
!>     A = U**T*U or A = L*L**T, as computed by SPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]CMODE
!>          CMODE is INTEGER
!>     Determines op2(C) in the formula op(A) * op2(C) as follows:
!>     CMODE =  1    op2(C) = C
!>     CMODE =  0    op2(C) = I
!>     CMODE = -1    op2(C) = inv(C)
!>
[in]C
!>          C is REAL array, dimension (N)
!>     The vector C in the formula op(A) * op2(C).
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>     i > 0:  The ith argument is invalid.
!>
[out]WORK
!>          WORK is REAL array, dimension (3*N).
!>     Workspace.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 138 of file sla_porcond.f.

140*
141* -- LAPACK computational routine --
142* -- LAPACK is a software package provided by Univ. of Tennessee, --
143* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144*
145* .. Scalar Arguments ..
146 CHARACTER UPLO
147 INTEGER N, LDA, LDAF, INFO, CMODE
148 REAL A( LDA, * ), AF( LDAF, * ), WORK( * ),
149 $ C( * )
150* ..
151* .. Array Arguments ..
152 INTEGER IWORK( * )
153* ..
154*
155* =====================================================================
156*
157* .. Local Scalars ..
158 INTEGER KASE, I, J
159 REAL AINVNM, TMP
160 LOGICAL UP
161* ..
162* .. Array Arguments ..
163 INTEGER ISAVE( 3 )
164* ..
165* .. External Functions ..
166 LOGICAL LSAME
167 EXTERNAL lsame
168* ..
169* .. External Subroutines ..
170 EXTERNAL slacn2, spotrs, xerbla
171* ..
172* .. Intrinsic Functions ..
173 INTRINSIC abs, max
174* ..
175* .. Executable Statements ..
176*
177 sla_porcond = 0.0
178*
179 info = 0
180 IF( n.LT.0 ) THEN
181 info = -2
182 END IF
183 IF( info.NE.0 ) THEN
184 CALL xerbla( 'SLA_PORCOND', -info )
185 RETURN
186 END IF
187
188 IF( n.EQ.0 ) THEN
189 sla_porcond = 1.0
190 RETURN
191 END IF
192 up = .false.
193 IF ( lsame( uplo, 'U' ) ) up = .true.
194*
195* Compute the equilibration matrix R such that
196* inv(R)*A*C has unit 1-norm.
197*
198 IF ( up ) THEN
199 DO i = 1, n
200 tmp = 0.0
201 IF ( cmode .EQ. 1 ) THEN
202 DO j = 1, i
203 tmp = tmp + abs( a( j, i ) * c( j ) )
204 END DO
205 DO j = i+1, n
206 tmp = tmp + abs( a( i, j ) * c( j ) )
207 END DO
208 ELSE IF ( cmode .EQ. 0 ) THEN
209 DO j = 1, i
210 tmp = tmp + abs( a( j, i ) )
211 END DO
212 DO j = i+1, n
213 tmp = tmp + abs( a( i, j ) )
214 END DO
215 ELSE
216 DO j = 1, i
217 tmp = tmp + abs( a( j ,i ) / c( j ) )
218 END DO
219 DO j = i+1, n
220 tmp = tmp + abs( a( i, j ) / c( j ) )
221 END DO
222 END IF
223 work( 2*n+i ) = tmp
224 END DO
225 ELSE
226 DO i = 1, n
227 tmp = 0.0
228 IF ( cmode .EQ. 1 ) THEN
229 DO j = 1, i
230 tmp = tmp + abs( a( i, j ) * c( j ) )
231 END DO
232 DO j = i+1, n
233 tmp = tmp + abs( a( j, i ) * c( j ) )
234 END DO
235 ELSE IF ( cmode .EQ. 0 ) THEN
236 DO j = 1, i
237 tmp = tmp + abs( a( i, j ) )
238 END DO
239 DO j = i+1, n
240 tmp = tmp + abs( a( j, i ) )
241 END DO
242 ELSE
243 DO j = 1, i
244 tmp = tmp + abs( a( i, j ) / c( j ) )
245 END DO
246 DO j = i+1, n
247 tmp = tmp + abs( a( j, i ) / c( j ) )
248 END DO
249 END IF
250 work( 2*n+i ) = tmp
251 END DO
252 ENDIF
253*
254* Estimate the norm of inv(op(A)).
255*
256 ainvnm = 0.0
257
258 kase = 0
259 10 CONTINUE
260 CALL slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
261 IF( kase.NE.0 ) THEN
262 IF( kase.EQ.2 ) THEN
263*
264* Multiply by R.
265*
266 DO i = 1, n
267 work( i ) = work( i ) * work( 2*n+i )
268 END DO
269
270 IF (up) THEN
271 CALL spotrs( 'Upper', n, 1, af, ldaf, work, n, info )
272 ELSE
273 CALL spotrs( 'Lower', n, 1, af, ldaf, work, n, info )
274 ENDIF
275*
276* Multiply by inv(C).
277*
278 IF ( cmode .EQ. 1 ) THEN
279 DO i = 1, n
280 work( i ) = work( i ) / c( i )
281 END DO
282 ELSE IF ( cmode .EQ. -1 ) THEN
283 DO i = 1, n
284 work( i ) = work( i ) * c( i )
285 END DO
286 END IF
287 ELSE
288*
289* Multiply by inv(C**T).
290*
291 IF ( cmode .EQ. 1 ) THEN
292 DO i = 1, n
293 work( i ) = work( i ) / c( i )
294 END DO
295 ELSE IF ( cmode .EQ. -1 ) THEN
296 DO i = 1, n
297 work( i ) = work( i ) * c( i )
298 END DO
299 END IF
300
301 IF ( up ) THEN
302 CALL spotrs( 'Upper', n, 1, af, ldaf, work, n, info )
303 ELSE
304 CALL spotrs( 'Lower', n, 1, af, ldaf, work, n, info )
305 ENDIF
306*
307* Multiply by R.
308*
309 DO i = 1, n
310 work( i ) = work( i ) * work( 2*n+i )
311 END DO
312 END IF
313 GO TO 10
314 END IF
315*
316* Compute the estimate of the reciprocal condition number.
317*
318 IF( ainvnm .NE. 0.0 )
319 $ sla_porcond = ( 1.0 / ainvnm )
320*
321 RETURN
322*
323* End of SLA_PORCOND
324*
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
slacn2
subroutine slacn2(n, v, x, isgn, est, kase, isave)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition slacn2.f:136
spotrs
subroutine spotrs(uplo, n, nrhs, a, lda, b, ldb, info)
SPOTRS
Definition spotrs.f:110
sla_porcond
real function sla_porcond(uplo, n, a, lda, af, ldaf, cmode, c, info, work, iwork)
SLA_PORCOND estimates the Skeel condition number for a symmetric positive-definite matrix.
Definition sla_porcond.f:140
max
#define max(a, b)
Definition macros.h:21

◆ sla_porfsx_extended()

subroutine sla_porfsx_extended ( integer prec_type,
character uplo,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldaf, * ) af,
integer ldaf,
logical colequ,
real, dimension( * ) c,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldy, * ) y,
integer ldy,
real, dimension( * ) berr_out,
integer n_norms,
real, dimension( nrhs, * ) err_bnds_norm,
real, dimension( nrhs, * ) err_bnds_comp,
real, dimension( * ) res,
real, dimension(*) ayb,
real, dimension( * ) dy,
real, dimension( * ) y_tail,
real rcond,
integer ithresh,
real rthresh,
real dz_ub,
logical ignore_cwise,
integer info )

SLA_PORFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric or Hermitian positive-definite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

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

Purpose:
!>
!> SLA_PORFSX_EXTENDED improves the computed solution to a system of
!> linear equations by performing extra-precise iterative refinement
!> and provides error bounds and backward error estimates for the solution.
!> This subroutine is called by SPORFSX to perform iterative refinement.
!> In addition to normwise error bound, the code provides maximum
!> componentwise error bound if possible. See comments for ERR_BNDS_NORM
!> and ERR_BNDS_COMP for details of the error bounds. Note that this
!> subroutine is only responsible for setting the second fields of
!> ERR_BNDS_NORM and ERR_BNDS_COMP.
!>
Parameters
[in]PREC_TYPE
!>          PREC_TYPE is INTEGER
!>     Specifies the intermediate precision to be used in refinement.
!>     The value is defined by ILAPREC(P) where P is a CHARACTER and P
!>          = 'S':  Single
!>          = 'D':  Double
!>          = 'I':  Indigenous
!>          = 'X' or 'E':  Extra
!>
[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.
!>
[in]A
!>          A is REAL array, dimension (LDA,N)
!>     On entry, the N-by-N matrix A.
!>
[in]LDA
!>          LDA is INTEGER
!>     The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]AF
!>          AF is REAL array, dimension (LDAF,N)
!>     The triangular factor U or L from the Cholesky factorization
!>     A = U**T*U or A = L*L**T, as computed by SPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]COLEQU
!>          COLEQU is LOGICAL
!>     If .TRUE. then column equilibration was done to A before calling
!>     this routine. This is needed to compute the solution and error
!>     bounds correctly.
!>
[in]C
!>          C is REAL array, dimension (N)
!>     The column scale factors for A. If COLEQU = .FALSE., C
!>     is not accessed. If C is input, each element of C should be a power
!>     of the radix to ensure a reliable solution and error estimates.
!>     Scaling by powers of the radix does not cause rounding errors unless
!>     the result underflows or overflows. Rounding errors during scaling
!>     lead to refining with a matrix that is not equivalent to the
!>     input matrix, producing error estimates that may not be
!>     reliable.
!>
[in]B
!>          B is REAL array, dimension (LDB,NRHS)
!>     The right-hand-side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>     The leading dimension of the array B.  LDB >= max(1,N).
!>
[in,out]Y
!>          Y is REAL array, dimension (LDY,NRHS)
!>     On entry, the solution matrix X, as computed by SPOTRS.
!>     On exit, the improved solution matrix Y.
!>
[in]LDY
!>          LDY is INTEGER
!>     The leading dimension of the array Y.  LDY >= max(1,N).
!>
[out]BERR_OUT
!>          BERR_OUT is REAL array, dimension (NRHS)
!>     On exit, BERR_OUT(j) contains the componentwise relative backward
!>     error for right-hand-side j from the formula
!>         max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
!>     where abs(Z) is the componentwise absolute value of the matrix
!>     or vector Z. This is computed by SLA_LIN_BERR.
!>
[in]N_NORMS
!>          N_NORMS is INTEGER
!>     Determines which error bounds to return (see ERR_BNDS_NORM
!>     and ERR_BNDS_COMP).
!>     If N_NORMS >= 1 return normwise error bounds.
!>     If N_NORMS >= 2 return componentwise error bounds.
!>
[in,out]ERR_BNDS_NORM
!>          ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     normwise relative error, which is defined as follows:
!>
!>     Normwise relative error in the ith solution vector:
!>             max_j (abs(XTRUE(j,i) - X(j,i)))
!>            ------------------------------
!>                  max_j abs(X(j,i))
!>
!>     The array is indexed by the type of error information as described
!>     below. There currently are up to three pieces of information
!>     returned.
!>
!>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_NORM(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     This subroutine is only responsible for setting the second field
!>     above.
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in,out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     componentwise relative error, which is defined as follows:
!>
!>     Componentwise relative error in the ith solution vector:
!>                    abs(XTRUE(j,i) - X(j,i))
!>             max_j ----------------------
!>                         abs(X(j,i))
!>
!>     The array is indexed by the right-hand side i (on which the
!>     componentwise relative error depends), and the type of error
!>     information as described below. There currently are up to three
!>     pieces of information returned for each right-hand side. If
!>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
!>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     This subroutine is only responsible for setting the second field
!>     above.
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]RES
!>          RES is REAL array, dimension (N)
!>     Workspace to hold the intermediate residual.
!>
[in]AYB
!>          AYB is REAL array, dimension (N)
!>     Workspace. This can be the same workspace passed for Y_TAIL.
!>
[in]DY
!>          DY is REAL array, dimension (N)
!>     Workspace to hold the intermediate solution.
!>
[in]Y_TAIL
!>          Y_TAIL is REAL array, dimension (N)
!>     Workspace to hold the trailing bits of the intermediate solution.
!>
[in]RCOND
!>          RCOND is REAL
!>     Reciprocal scaled condition number.  This is an estimate of the
!>     reciprocal Skeel condition number of the matrix A after
!>     equilibration (if done).  If this is less than the machine
!>     precision (in particular, if it is zero), the matrix is singular
!>     to working precision.  Note that the error may still be small even
!>     if this number is very small and the matrix appears ill-
!>     conditioned.
!>
[in]ITHRESH
!>          ITHRESH is INTEGER
!>     The maximum number of residual computations allowed for
!>     refinement. The default is 10. For 'aggressive' set to 100 to
!>     permit convergence using approximate factorizations or
!>     factorizations other than LU. If the factorization uses a
!>     technique other than Gaussian elimination, the guarantees in
!>     ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
!>
[in]RTHRESH
!>          RTHRESH is REAL
!>     Determines when to stop refinement if the error estimate stops
!>     decreasing. Refinement will stop when the next solution no longer
!>     satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
!>     the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
!>     default value is 0.5. For 'aggressive' set to 0.9 to permit
!>     convergence on extremely ill-conditioned matrices. See LAWN 165
!>     for more details.
!>
[in]DZ_UB
!>          DZ_UB is REAL
!>     Determines when to start considering componentwise convergence.
!>     Componentwise convergence is only considered after each component
!>     of the solution Y is stable, which we define as the relative
!>     change in each component being less than DZ_UB. The default value
!>     is 0.25, requiring the first bit to be stable. See LAWN 165 for
!>     more details.
!>
[in]IGNORE_CWISE
!>          IGNORE_CWISE is LOGICAL
!>     If .TRUE. then ignore componentwise convergence. Default value
!>     is .FALSE..
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>       < 0:  if INFO = -i, the ith argument to SPOTRS had an illegal
!>             value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 380 of file sla_porfsx_extended.f.

387*
388* -- LAPACK computational routine --
389* -- LAPACK is a software package provided by Univ. of Tennessee, --
390* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
391*
392* .. Scalar Arguments ..
393 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
394 $ N_NORMS, ITHRESH
395 CHARACTER UPLO
396 LOGICAL COLEQU, IGNORE_CWISE
397 REAL RTHRESH, DZ_UB
398* ..
399* .. Array Arguments ..
400 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
401 $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
402 REAL C( * ), AYB(*), RCOND, BERR_OUT( * ),
403 $ ERR_BNDS_NORM( NRHS, * ),
404 $ ERR_BNDS_COMP( NRHS, * )
405* ..
406*
407* =====================================================================
408*
409* .. Local Scalars ..
410 INTEGER UPLO2, CNT, I, J, X_STATE, Z_STATE
411 REAL YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
412 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
413 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
414 $ EPS, HUGEVAL, INCR_THRESH
415 LOGICAL INCR_PREC
416* ..
417* .. Parameters ..
418 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
419 $ NOPROG_STATE, Y_PREC_STATE, BASE_RESIDUAL,
420 $ EXTRA_RESIDUAL, EXTRA_Y
421 parameter( unstable_state = 0, working_state = 1,
422 $ conv_state = 2, noprog_state = 3 )
423 parameter( base_residual = 0, extra_residual = 1,
424 $ extra_y = 2 )
425 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
426 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
427 INTEGER CMP_ERR_I, PIV_GROWTH_I
428 parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
429 $ berr_i = 3 )
430 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
431 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
432 $ piv_growth_i = 9 )
433 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
434 $ LA_LINRX_CWISE_I
435 parameter( la_linrx_itref_i = 1,
436 $ la_linrx_ithresh_i = 2 )
437 parameter( la_linrx_cwise_i = 3 )
438 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
439 $ LA_LINRX_RCOND_I
440 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
441 parameter( la_linrx_rcond_i = 3 )
442* ..
443* .. External Functions ..
444 LOGICAL LSAME
445 EXTERNAL ilauplo
446 INTEGER ILAUPLO
447* ..
448* .. External Subroutines ..
449 EXTERNAL saxpy, scopy, spotrs, ssymv, blas_ssymv_x,
450 $ blas_ssymv2_x, sla_syamv, sla_wwaddw,
451 $ sla_lin_berr
452 REAL SLAMCH
453* ..
454* .. Intrinsic Functions ..
455 INTRINSIC abs, max, min
456* ..
457* .. Executable Statements ..
458*
459 IF (info.NE.0) RETURN
460 eps = slamch( 'Epsilon' )
461 hugeval = slamch( 'Overflow' )
462* Force HUGEVAL to Inf
463 hugeval = hugeval * hugeval
464* Using HUGEVAL may lead to spurious underflows.
465 incr_thresh = real( n ) * eps
466
467 IF ( lsame( uplo, 'L' ) ) THEN
468 uplo2 = ilauplo( 'L' )
469 ELSE
470 uplo2 = ilauplo( 'U' )
471 ENDIF
472
473 DO j = 1, nrhs
474 y_prec_state = extra_residual
475 IF ( y_prec_state .EQ. extra_y ) THEN
476 DO i = 1, n
477 y_tail( i ) = 0.0
478 END DO
479 END IF
480
481 dxrat = 0.0
482 dxratmax = 0.0
483 dzrat = 0.0
484 dzratmax = 0.0
485 final_dx_x = hugeval
486 final_dz_z = hugeval
487 prevnormdx = hugeval
488 prev_dz_z = hugeval
489 dz_z = hugeval
490 dx_x = hugeval
491
492 x_state = working_state
493 z_state = unstable_state
494 incr_prec = .false.
495
496 DO cnt = 1, ithresh
497*
498* Compute residual RES = B_s - op(A_s) * Y,
499* op(A) = A, A**T, or A**H depending on TRANS (and type).
500*
501 CALL scopy( n, b( 1, j ), 1, res, 1 )
502 IF ( y_prec_state .EQ. base_residual ) THEN
503 CALL ssymv( uplo, n, -1.0, a, lda, y(1,j), 1,
504 $ 1.0, res, 1 )
505 ELSE IF ( y_prec_state .EQ. extra_residual ) THEN
506 CALL blas_ssymv_x( uplo2, n, -1.0, a, lda,
507 $ y( 1, j ), 1, 1.0, res, 1, prec_type )
508 ELSE
509 CALL blas_ssymv2_x(uplo2, n, -1.0, a, lda,
510 $ y(1, j), y_tail, 1, 1.0, res, 1, prec_type)
511 END IF
512
513! XXX: RES is no longer needed.
514 CALL scopy( n, res, 1, dy, 1 )
515 CALL spotrs( uplo, n, 1, af, ldaf, dy, n, info )
516*
517* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
518*
519 normx = 0.0
520 normy = 0.0
521 normdx = 0.0
522 dz_z = 0.0
523 ymin = hugeval
524
525 DO i = 1, n
526 yk = abs( y( i, j ) )
527 dyk = abs( dy( i ) )
528
529 IF ( yk .NE. 0.0 ) THEN
530 dz_z = max( dz_z, dyk / yk )
531 ELSE IF ( dyk .NE. 0.0 ) THEN
532 dz_z = hugeval
533 END IF
534
535 ymin = min( ymin, yk )
536
537 normy = max( normy, yk )
538
539 IF ( colequ ) THEN
540 normx = max( normx, yk * c( i ) )
541 normdx = max( normdx, dyk * c( i ) )
542 ELSE
543 normx = normy
544 normdx = max( normdx, dyk )
545 END IF
546 END DO
547
548 IF ( normx .NE. 0.0 ) THEN
549 dx_x = normdx / normx
550 ELSE IF ( normdx .EQ. 0.0 ) THEN
551 dx_x = 0.0
552 ELSE
553 dx_x = hugeval
554 END IF
555
556 dxrat = normdx / prevnormdx
557 dzrat = dz_z / prev_dz_z
558*
559* Check termination criteria.
560*
561 IF ( ymin*rcond .LT. incr_thresh*normy
562 $ .AND. y_prec_state .LT. extra_y )
563 $ incr_prec = .true.
564
565 IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
566 $ x_state = working_state
567 IF ( x_state .EQ. working_state ) THEN
568 IF ( dx_x .LE. eps ) THEN
569 x_state = conv_state
570 ELSE IF ( dxrat .GT. rthresh ) THEN
571 IF ( y_prec_state .NE. extra_y ) THEN
572 incr_prec = .true.
573 ELSE
574 x_state = noprog_state
575 END IF
576 ELSE
577 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
578 END IF
579 IF ( x_state .GT. working_state ) final_dx_x = dx_x
580 END IF
581
582 IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
583 $ z_state = working_state
584 IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
585 $ z_state = working_state
586 IF ( z_state .EQ. working_state ) THEN
587 IF ( dz_z .LE. eps ) THEN
588 z_state = conv_state
589 ELSE IF ( dz_z .GT. dz_ub ) THEN
590 z_state = unstable_state
591 dzratmax = 0.0
592 final_dz_z = hugeval
593 ELSE IF ( dzrat .GT. rthresh ) THEN
594 IF ( y_prec_state .NE. extra_y ) THEN
595 incr_prec = .true.
596 ELSE
597 z_state = noprog_state
598 END IF
599 ELSE
600 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
601 END IF
602 IF ( z_state .GT. working_state ) final_dz_z = dz_z
603 END IF
604
605 IF ( x_state.NE.working_state.AND.
606 $ ( ignore_cwise.OR.z_state.NE.working_state ) )
607 $ GOTO 666
608
609 IF ( incr_prec ) THEN
610 incr_prec = .false.
611 y_prec_state = y_prec_state + 1
612 DO i = 1, n
613 y_tail( i ) = 0.0
614 END DO
615 END IF
616
617 prevnormdx = normdx
618 prev_dz_z = dz_z
619*
620* Update soluton.
621*
622 IF (y_prec_state .LT. extra_y) THEN
623 CALL saxpy( n, 1.0, dy, 1, y(1,j), 1 )
624 ELSE
625 CALL sla_wwaddw( n, y( 1, j ), y_tail, dy )
626 END IF
627
628 END DO
629* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
630 666 CONTINUE
631*
632* Set final_* when cnt hits ithresh.
633*
634 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
635 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
636*
637* Compute error bounds.
638*
639 IF ( n_norms .GE. 1 ) THEN
640 err_bnds_norm( j, la_linrx_err_i ) =
641 $ final_dx_x / (1 - dxratmax)
642 END IF
643 IF ( n_norms .GE. 2 ) THEN
644 err_bnds_comp( j, la_linrx_err_i ) =
645 $ final_dz_z / (1 - dzratmax)
646 END IF
647*
648* Compute componentwise relative backward error from formula
649* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
650* where abs(Z) is the componentwise absolute value of the matrix
651* or vector Z.
652*
653* Compute residual RES = B_s - op(A_s) * Y,
654* op(A) = A, A**T, or A**H depending on TRANS (and type).
655*
656 CALL scopy( n, b( 1, j ), 1, res, 1 )
657 CALL ssymv( uplo, n, -1.0, a, lda, y(1,j), 1, 1.0, res, 1 )
658
659 DO i = 1, n
660 ayb( i ) = abs( b( i, j ) )
661 END DO
662*
663* Compute abs(op(A_s))*abs(Y) + abs(B_s).
664*
665 CALL sla_syamv( uplo2, n, 1.0,
666 $ a, lda, y(1, j), 1, 1.0, ayb, 1 )
667
668 CALL sla_lin_berr( n, n, 1, res, ayb, berr_out( j ) )
669*
670* End of loop for each RHS.
671*
672 END DO
673*
674 RETURN
675*
676* End of SLA_PORFSX_EXTENDED
677*
ilauplo
integer function ilauplo(uplo)
ILAUPLO
Definition ilauplo.f:58
sla_lin_berr
subroutine sla_lin_berr(n, nz, nrhs, res, ayb, berr)
SLA_LIN_BERR computes a component-wise relative backward error.
Definition sla_lin_berr.f:101
sla_wwaddw
subroutine sla_wwaddw(n, x, y, w)
SLA_WWADDW adds a vector into a doubled-single vector.
Definition sla_wwaddw.f:81
sla_syamv
subroutine sla_syamv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
SLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bou...
Definition sla_syamv.f:177
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
ssymv
subroutine ssymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
SSYMV
Definition ssymv.f:152
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
min
#define min(a, b)
Definition macros.h:20

◆ sla_porpvgrw()

real function sla_porpvgrw ( character*1 uplo,
integer ncols,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldaf, * ) af,
integer ldaf,
real, dimension( * ) work )

SLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian positive-definite matrix.

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

Purpose:
!>
!>
!> SLA_PORPVGRW computes the reciprocal pivot growth factor
!> norm(A)/norm(U). The  norm is used. If this is
!> much less than 1, the stability of the LU factorization of the
!> (equilibrated) matrix A could be poor. This also means that the
!> solution X, estimated condition numbers, and error bounds could be
!> unreliable.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>       = 'U':  Upper triangle of A is stored;
!>       = 'L':  Lower triangle of A is stored.
!>
[in]NCOLS
!>          NCOLS is INTEGER
!>     The number of columns of the matrix A. NCOLS >= 0.
!>
[in]A
!>          A is REAL array, dimension (LDA,N)
!>     On entry, the N-by-N matrix A.
!>
[in]LDA
!>          LDA is INTEGER
!>     The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]AF
!>          AF is REAL array, dimension (LDAF,N)
!>     The triangular factor U or L from the Cholesky factorization
!>     A = U**T*U or A = L*L**T, as computed by SPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[out]WORK
!>          WORK is REAL array, dimension (2*N)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 103 of file sla_porpvgrw.f.

104*
105* -- LAPACK computational routine --
106* -- LAPACK is a software package provided by Univ. of Tennessee, --
107* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
108*
109* .. Scalar Arguments ..
110 CHARACTER*1 UPLO
111 INTEGER NCOLS, LDA, LDAF
112* ..
113* .. Array Arguments ..
114 REAL A( LDA, * ), AF( LDAF, * ), WORK( * )
115* ..
116*
117* =====================================================================
118*
119* .. Local Scalars ..
120 INTEGER I, J
121 REAL AMAX, UMAX, RPVGRW
122 LOGICAL UPPER
123* ..
124* .. Intrinsic Functions ..
125 INTRINSIC abs, max, min
126* ..
127* .. External Functions ..
128 EXTERNAL lsame
129 LOGICAL LSAME
130* ..
131* .. Executable Statements ..
132*
133 upper = lsame( 'Upper', uplo )
134*
135* SPOTRF will have factored only the NCOLSxNCOLS leading minor, so
136* we restrict the growth search to that minor and use only the first
137* 2*NCOLS workspace entries.
138*
139 rpvgrw = 1.0
140 DO i = 1, 2*ncols
141 work( i ) = 0.0
142 END DO
143*
144* Find the max magnitude entry of each column.
145*
146 IF ( upper ) THEN
147 DO j = 1, ncols
148 DO i = 1, j
149 work( ncols+j ) =
150 $ max( abs( a( i, j ) ), work( ncols+j ) )
151 END DO
152 END DO
153 ELSE
154 DO j = 1, ncols
155 DO i = j, ncols
156 work( ncols+j ) =
157 $ max( abs( a( i, j ) ), work( ncols+j ) )
158 END DO
159 END DO
160 END IF
161*
162* Now find the max magnitude entry of each column of the factor in
163* AF. No pivoting, so no permutations.
164*
165 IF ( lsame( 'Upper', uplo ) ) THEN
166 DO j = 1, ncols
167 DO i = 1, j
168 work( j ) = max( abs( af( i, j ) ), work( j ) )
169 END DO
170 END DO
171 ELSE
172 DO j = 1, ncols
173 DO i = j, ncols
174 work( j ) = max( abs( af( i, j ) ), work( j ) )
175 END DO
176 END DO
177 END IF
178*
179* Compute the *inverse* of the max element growth factor. Dividing
180* by zero would imply the largest entry of the factor's column is
181* zero. Than can happen when either the column of A is zero or
182* massive pivots made the factor underflow to zero. Neither counts
183* as growth in itself, so simply ignore terms with zero
184* denominators.
185*
186 IF ( lsame( 'Upper', uplo ) ) THEN
187 DO i = 1, ncols
188 umax = work( i )
189 amax = work( ncols+i )
190 IF ( umax /= 0.0 ) THEN
191 rpvgrw = min( amax / umax, rpvgrw )
192 END IF
193 END DO
194 ELSE
195 DO i = 1, ncols
196 umax = work( i )
197 amax = work( ncols+i )
198 IF ( umax /= 0.0 ) THEN
199 rpvgrw = min( amax / umax, rpvgrw )
200 END IF
201 END DO
202 END IF
203
204 sla_porpvgrw = rpvgrw
205*
206* End of SLA_PORPVGRW
207*
sla_porpvgrw
real function sla_porpvgrw(uplo, ncols, a, lda, af, ldaf, work)
SLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian...
Definition sla_porpvgrw.f:104

◆ spocon()

subroutine spocon ( character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real anorm,
real rcond,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SPOCON

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

Purpose:
!>
!> SPOCON estimates the reciprocal of the condition number (in the
!> 1-norm) of a real symmetric positive definite matrix using the
!> Cholesky factorization A = U**T*U or A = L*L**T computed by SPOTRF.
!>
!> An estimate is obtained for norm(inv(A)), and the reciprocal of the
!> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
!>
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 order of the matrix A.  N >= 0.
!>
[in]A
!>          A is REAL array, dimension (LDA,N)
!>          The triangular factor U or L from the Cholesky factorization
!>          A = U**T*U or A = L*L**T, as computed by SPOTRF.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]ANORM
!>          ANORM is REAL
!>          The 1-norm (or infinity-norm) of the symmetric matrix A.
!>
[out]RCOND
!>          RCOND is REAL
!>          The reciprocal of the condition number of the matrix A,
!>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
!>          estimate of the 1-norm of inv(A) computed in this routine.
!>
[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
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 119 of file spocon.f.

121*
122* -- LAPACK computational routine --
123* -- LAPACK is a software package provided by Univ. of Tennessee, --
124* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
125*
126* .. Scalar Arguments ..
127 CHARACTER UPLO
128 INTEGER INFO, LDA, N
129 REAL ANORM, RCOND
130* ..
131* .. Array Arguments ..
132 INTEGER IWORK( * )
133 REAL A( LDA, * ), WORK( * )
134* ..
135*
136* =====================================================================
137*
138* .. Parameters ..
139 REAL ONE, ZERO
140 parameter( one = 1.0e+0, zero = 0.0e+0 )
141* ..
142* .. Local Scalars ..
143 LOGICAL UPPER
144 CHARACTER NORMIN
145 INTEGER IX, KASE
146 REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
147* ..
148* .. Local Arrays ..
149 INTEGER ISAVE( 3 )
150* ..
151* .. External Functions ..
152 LOGICAL LSAME
153 INTEGER ISAMAX
154 REAL SLAMCH
155 EXTERNAL lsame, isamax, slamch
156* ..
157* .. External Subroutines ..
158 EXTERNAL slacn2, slatrs, srscl, xerbla
159* ..
160* .. Intrinsic Functions ..
161 INTRINSIC abs, max
162* ..
163* .. Executable Statements ..
164*
165* Test the input parameters.
166*
167 info = 0
168 upper = lsame( uplo, 'U' )
169 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
170 info = -1
171 ELSE IF( n.LT.0 ) THEN
172 info = -2
173 ELSE IF( lda.LT.max( 1, n ) ) THEN
174 info = -4
175 ELSE IF( anorm.LT.zero ) THEN
176 info = -5
177 END IF
178 IF( info.NE.0 ) THEN
179 CALL xerbla( 'SPOCON', -info )
180 RETURN
181 END IF
182*
183* Quick return if possible
184*
185 rcond = zero
186 IF( n.EQ.0 ) THEN
187 rcond = one
188 RETURN
189 ELSE IF( anorm.EQ.zero ) THEN
190 RETURN
191 END IF
192*
193 smlnum = slamch( 'Safe minimum' )
194*
195* Estimate the 1-norm of inv(A).
196*
197 kase = 0
198 normin = 'N'
199 10 CONTINUE
200 CALL slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
201 IF( kase.NE.0 ) THEN
202 IF( upper ) THEN
203*
204* Multiply by inv(U**T).
205*
206 CALL slatrs( 'Upper', 'Transpose', 'Non-unit', normin, n, a,
207 $ lda, work, scalel, work( 2*n+1 ), info )
208 normin = 'Y'
209*
210* Multiply by inv(U).
211*
212 CALL slatrs( 'Upper', 'No transpose', 'Non-unit', normin, n,
213 $ a, lda, work, scaleu, work( 2*n+1 ), info )
214 ELSE
215*
216* Multiply by inv(L).
217*
218 CALL slatrs( 'Lower', 'No transpose', 'Non-unit', normin, n,
219 $ a, lda, work, scalel, work( 2*n+1 ), info )
220 normin = 'Y'
221*
222* Multiply by inv(L**T).
223*
224 CALL slatrs( 'Lower', 'Transpose', 'Non-unit', normin, n, a,
225 $ lda, work, scaleu, work( 2*n+1 ), info )
226 END IF
227*
228* Multiply by 1/SCALE if doing so will not cause overflow.
229*
230 scale = scalel*scaleu
231 IF( scale.NE.one ) THEN
232 ix = isamax( n, work, 1 )
233 IF( scale.LT.abs( work( ix ) )*smlnum .OR. scale.EQ.zero )
234 $ GO TO 20
235 CALL srscl( n, scale, work, 1 )
236 END IF
237 GO TO 10
238 END IF
239*
240* Compute the estimate of the reciprocal condition number.
241*
242 IF( ainvnm.NE.zero )
243 $ rcond = ( one / ainvnm ) / anorm
244*
245 20 CONTINUE
246 RETURN
247*
248* End of SPOCON
249*
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slatrs
subroutine slatrs(uplo, trans, diag, normin, n, a, lda, x, scale, cnorm, info)
SLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition slatrs.f:238
srscl
subroutine srscl(n, sa, sx, incx)
SRSCL multiplies a vector by the reciprocal of a real scalar.
Definition srscl.f:84

◆ spoequ()

subroutine spoequ ( integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) s,
real scond,
real amax,
integer info )

SPOEQU

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

Purpose:
!>
!> SPOEQU computes row and column scalings intended to equilibrate a
!> symmetric positive definite matrix A and reduce its condition number
!> (with respect to the two-norm).  S contains the scale factors,
!> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
!> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
!> choice of S puts the condition number of B within a factor N of the
!> smallest possible condition number over all possible diagonal
!> scalings.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]A
!>          A is REAL array, dimension (LDA,N)
!>          The N-by-N symmetric positive definite matrix whose scaling
!>          factors are to be computed.  Only the diagonal elements of A
!>          are referenced.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]S
!>          S is REAL array, dimension (N)
!>          If INFO = 0, S contains the scale factors for A.
!>
[out]SCOND
!>          SCOND is REAL
!>          If INFO = 0, S contains the ratio of the smallest S(i) to
!>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
!>          large nor too small, it is not worth scaling by S.
!>
[out]AMAX
!>          AMAX is REAL
!>          Absolute value of largest matrix element.  If AMAX is very
!>          close to overflow or very close to underflow, the matrix
!>          should be scaled.
!>
[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 i-th diagonal element is nonpositive.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 111 of file spoequ.f.

112*
113* -- LAPACK computational routine --
114* -- LAPACK is a software package provided by Univ. of Tennessee, --
115* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
116*
117* .. Scalar Arguments ..
118 INTEGER INFO, LDA, N
119 REAL AMAX, SCOND
120* ..
121* .. Array Arguments ..
122 REAL A( LDA, * ), S( * )
123* ..
124*
125* =====================================================================
126*
127* .. Parameters ..
128 REAL ZERO, ONE
129 parameter( zero = 0.0e+0, one = 1.0e+0 )
130* ..
131* .. Local Scalars ..
132 INTEGER I
133 REAL SMIN
134* ..
135* .. External Subroutines ..
136 EXTERNAL xerbla
137* ..
138* .. Intrinsic Functions ..
139 INTRINSIC max, min, sqrt
140* ..
141* .. Executable Statements ..
142*
143* Test the input parameters.
144*
145 info = 0
146 IF( n.LT.0 ) THEN
147 info = -1
148 ELSE IF( lda.LT.max( 1, n ) ) THEN
149 info = -3
150 END IF
151 IF( info.NE.0 ) THEN
152 CALL xerbla( 'SPOEQU', -info )
153 RETURN
154 END IF
155*
156* Quick return if possible
157*
158 IF( n.EQ.0 ) THEN
159 scond = one
160 amax = zero
161 RETURN
162 END IF
163*
164* Find the minimum and maximum diagonal elements.
165*
166 s( 1 ) = a( 1, 1 )
167 smin = s( 1 )
168 amax = s( 1 )
169 DO 10 i = 2, n
170 s( i ) = a( i, i )
171 smin = min( smin, s( i ) )
172 amax = max( amax, s( i ) )
173 10 CONTINUE
174*
175 IF( smin.LE.zero ) THEN
176*
177* Find the first non-positive diagonal element and return.
178*
179 DO 20 i = 1, n
180 IF( s( i ).LE.zero ) THEN
181 info = i
182 RETURN
183 END IF
184 20 CONTINUE
185 ELSE
186*
187* Set the scale factors to the reciprocals
188* of the diagonal elements.
189*
190 DO 30 i = 1, n
191 s( i ) = one / sqrt( s( i ) )
192 30 CONTINUE
193*
194* Compute SCOND = min(S(I)) / max(S(I))
195*
196 scond = sqrt( smin ) / sqrt( amax )
197 END IF
198 RETURN
199*
200* End of SPOEQU
201*

◆ spoequb()

subroutine spoequb ( integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) s,
real scond,
real amax,
integer info )

SPOEQUB

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

Purpose:
!>
!> SPOEQUB computes row and column scalings intended to equilibrate a
!> symmetric positive definite matrix A and reduce its condition number
!> (with respect to the two-norm).  S contains the scale factors,
!> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
!> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
!> choice of S puts the condition number of B within a factor N of the
!> smallest possible condition number over all possible diagonal
!> scalings.
!>
!> This routine differs from SPOEQU by restricting the scaling factors
!> to a power of the radix.  Barring over- and underflow, scaling by
!> these factors introduces no additional rounding errors.  However, the
!> scaled diagonal entries are no longer approximately 1 but lie
!> between sqrt(radix) and 1/sqrt(radix).
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]A
!>          A is REAL array, dimension (LDA,N)
!>          The N-by-N symmetric positive definite matrix whose scaling
!>          factors are to be computed.  Only the diagonal elements of A
!>          are referenced.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]S
!>          S is REAL array, dimension (N)
!>          If INFO = 0, S contains the scale factors for A.
!>
[out]SCOND
!>          SCOND is REAL
!>          If INFO = 0, S contains the ratio of the smallest S(i) to
!>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
!>          large nor too small, it is not worth scaling by S.
!>
[out]AMAX
!>          AMAX is REAL
!>          Absolute value of largest matrix element.  If AMAX is very
!>          close to overflow or very close to underflow, the matrix
!>          should be scaled.
!>
[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 i-th diagonal element is nonpositive.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 117 of file spoequb.f.

118*
119* -- LAPACK computational routine --
120* -- LAPACK is a software package provided by Univ. of Tennessee, --
121* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
122*
123* .. Scalar Arguments ..
124 INTEGER INFO, LDA, N
125 REAL AMAX, SCOND
126* ..
127* .. Array Arguments ..
128 REAL A( LDA, * ), S( * )
129* ..
130*
131* =====================================================================
132*
133* .. Parameters ..
134 REAL ZERO, ONE
135 parameter( zero = 0.0e+0, one = 1.0e+0 )
136* ..
137* .. Local Scalars ..
138 INTEGER I
139 REAL SMIN, BASE, TMP
140* ..
141* .. External Functions ..
142 REAL SLAMCH
143 EXTERNAL slamch
144* ..
145* .. External Subroutines ..
146 EXTERNAL xerbla
147* ..
148* .. Intrinsic Functions ..
149 INTRINSIC max, min, sqrt, log, int
150* ..
151* .. Executable Statements ..
152*
153* Test the input parameters.
154*
155* Positive definite only performs 1 pass of equilibration.
156*
157 info = 0
158 IF( n.LT.0 ) THEN
159 info = -1
160 ELSE IF( lda.LT.max( 1, n ) ) THEN
161 info = -3
162 END IF
163 IF( info.NE.0 ) THEN
164 CALL xerbla( 'SPOEQUB', -info )
165 RETURN
166 END IF
167*
168* Quick return if possible.
169*
170 IF( n.EQ.0 ) THEN
171 scond = one
172 amax = zero
173 RETURN
174 END IF
175
176 base = slamch( 'B' )
177 tmp = -0.5 / log( base )
178*
179* Find the minimum and maximum diagonal elements.
180*
181 s( 1 ) = a( 1, 1 )
182 smin = s( 1 )
183 amax = s( 1 )
184 DO 10 i = 2, n
185 s( i ) = a( i, i )
186 smin = min( smin, s( i ) )
187 amax = max( amax, s( i ) )
188 10 CONTINUE
189*
190 IF( smin.LE.zero ) THEN
191*
192* Find the first non-positive diagonal element and return.
193*
194 DO 20 i = 1, n
195 IF( s( i ).LE.zero ) THEN
196 info = i
197 RETURN
198 END IF
199 20 CONTINUE
200 ELSE
201*
202* Set the scale factors to the reciprocals
203* of the diagonal elements.
204*
205 DO 30 i = 1, n
206 s( i ) = base ** int( tmp * log( s( i ) ) )
207 30 CONTINUE
208*
209* Compute SCOND = min(S(I)) / max(S(I)).
210*
211 scond = sqrt( smin ) / sqrt( amax )
212 END IF
213*
214 RETURN
215*
216* End of SPOEQUB
217*

◆ sporfs()

subroutine sporfs ( character uplo,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldaf, * ) af,
integer ldaf,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldx, * ) x,
integer ldx,
real, dimension( * ) ferr,
real, dimension( * ) berr,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SPORFS

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

Purpose:
!>
!> SPORFS improves the computed solution to a system of linear
!> equations when the coefficient matrix is symmetric positive definite,
!> and provides error bounds and backward error estimates for the
!> solution.
!>
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 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]A
!>          A is REAL array, dimension (LDA,N)
!>          The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
!>          upper triangular part of A contains the upper triangular part
!>          of the matrix A, and the strictly lower triangular part of A
!>          is not referenced.  If UPLO = 'L', the leading N-by-N lower
!>          triangular part of A contains the lower triangular part of
!>          the matrix A, and the strictly upper triangular part of A is
!>          not referenced.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]AF
!>          AF is REAL array, dimension (LDAF,N)
!>          The triangular factor U or L from the Cholesky factorization
!>          A = U**T*U or A = L*L**T, as computed by SPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>          The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          The right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[in,out]X
!>          X is REAL array, dimension (LDX,NRHS)
!>          On entry, the solution matrix X, as computed by SPOTRS.
!>          On exit, the improved solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[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
!>
Internal Parameters:
!>  ITMAX is the maximum number of steps of iterative refinement.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 181 of file sporfs.f.

183*
184* -- LAPACK computational routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER UPLO
190 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
191* ..
192* .. Array Arguments ..
193 INTEGER IWORK( * )
194 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
195 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
196* ..
197*
198* =====================================================================
199*
200* .. Parameters ..
201 INTEGER ITMAX
202 parameter( itmax = 5 )
203 REAL ZERO
204 parameter( zero = 0.0e+0 )
205 REAL ONE
206 parameter( one = 1.0e+0 )
207 REAL TWO
208 parameter( two = 2.0e+0 )
209 REAL THREE
210 parameter( three = 3.0e+0 )
211* ..
212* .. Local Scalars ..
213 LOGICAL UPPER
214 INTEGER COUNT, I, J, K, KASE, NZ
215 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
216* ..
217* .. Local Arrays ..
218 INTEGER ISAVE( 3 )
219* ..
220* .. External Subroutines ..
221 EXTERNAL saxpy, scopy, slacn2, spotrs, ssymv, xerbla
222* ..
223* .. Intrinsic Functions ..
224 INTRINSIC abs, max
225* ..
226* .. External Functions ..
227 LOGICAL LSAME
228 REAL SLAMCH
229 EXTERNAL lsame, slamch
230* ..
231* .. Executable Statements ..
232*
233* Test the input parameters.
234*
235 info = 0
236 upper = lsame( uplo, 'U' )
237 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
238 info = -1
239 ELSE IF( n.LT.0 ) THEN
240 info = -2
241 ELSE IF( nrhs.LT.0 ) THEN
242 info = -3
243 ELSE IF( lda.LT.max( 1, n ) ) THEN
244 info = -5
245 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
246 info = -7
247 ELSE IF( ldb.LT.max( 1, n ) ) THEN
248 info = -9
249 ELSE IF( ldx.LT.max( 1, n ) ) THEN
250 info = -11
251 END IF
252 IF( info.NE.0 ) THEN
253 CALL xerbla( 'SPORFS', -info )
254 RETURN
255 END IF
256*
257* Quick return if possible
258*
259 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
260 DO 10 j = 1, nrhs
261 ferr( j ) = zero
262 berr( j ) = zero
263 10 CONTINUE
264 RETURN
265 END IF
266*
267* NZ = maximum number of nonzero elements in each row of A, plus 1
268*
269 nz = n + 1
270 eps = slamch( 'Epsilon' )
271 safmin = slamch( 'Safe minimum' )
272 safe1 = nz*safmin
273 safe2 = safe1 / eps
274*
275* Do for each right hand side
276*
277 DO 140 j = 1, nrhs
278*
279 count = 1
280 lstres = three
281 20 CONTINUE
282*
283* Loop until stopping criterion is satisfied.
284*
285* Compute residual R = B - A * X
286*
287 CALL scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
288 CALL ssymv( uplo, n, -one, a, lda, x( 1, j ), 1, one,
289 $ work( n+1 ), 1 )
290*
291* Compute componentwise relative backward error from formula
292*
293* max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
294*
295* where abs(Z) is the componentwise absolute value of the matrix
296* or vector Z. If the i-th component of the denominator is less
297* than SAFE2, then SAFE1 is added to the i-th components of the
298* numerator and denominator before dividing.
299*
300 DO 30 i = 1, n
301 work( i ) = abs( b( i, j ) )
302 30 CONTINUE
303*
304* Compute abs(A)*abs(X) + abs(B).
305*
306 IF( upper ) THEN
307 DO 50 k = 1, n
308 s = zero
309 xk = abs( x( k, j ) )
310 DO 40 i = 1, k - 1
311 work( i ) = work( i ) + abs( a( i, k ) )*xk
312 s = s + abs( a( i, k ) )*abs( x( i, j ) )
313 40 CONTINUE
314 work( k ) = work( k ) + abs( a( k, k ) )*xk + s
315 50 CONTINUE
316 ELSE
317 DO 70 k = 1, n
318 s = zero
319 xk = abs( x( k, j ) )
320 work( k ) = work( k ) + abs( a( k, k ) )*xk
321 DO 60 i = k + 1, n
322 work( i ) = work( i ) + abs( a( i, k ) )*xk
323 s = s + abs( a( i, k ) )*abs( x( i, j ) )
324 60 CONTINUE
325 work( k ) = work( k ) + s
326 70 CONTINUE
327 END IF
328 s = zero
329 DO 80 i = 1, n
330 IF( work( i ).GT.safe2 ) THEN
331 s = max( s, abs( work( n+i ) ) / work( i ) )
332 ELSE
333 s = max( s, ( abs( work( n+i ) )+safe1 ) /
334 $ ( work( i )+safe1 ) )
335 END IF
336 80 CONTINUE
337 berr( j ) = s
338*
339* Test stopping criterion. Continue iterating if
340* 1) The residual BERR(J) is larger than machine epsilon, and
341* 2) BERR(J) decreased by at least a factor of 2 during the
342* last iteration, and
343* 3) At most ITMAX iterations tried.
344*
345 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
346 $ count.LE.itmax ) THEN
347*
348* Update solution and try again.
349*
350 CALL spotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
351 CALL saxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
352 lstres = berr( j )
353 count = count + 1
354 GO TO 20
355 END IF
356*
357* Bound error from formula
358*
359* norm(X - XTRUE) / norm(X) .le. FERR =
360* norm( abs(inv(A))*
361* ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
362*
363* where
364* norm(Z) is the magnitude of the largest component of Z
365* inv(A) is the inverse of A
366* abs(Z) is the componentwise absolute value of the matrix or
367* vector Z
368* NZ is the maximum number of nonzeros in any row of A, plus 1
369* EPS is machine epsilon
370*
371* The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
372* is incremented by SAFE1 if the i-th component of
373* abs(A)*abs(X) + abs(B) is less than SAFE2.
374*
375* Use SLACN2 to estimate the infinity-norm of the matrix
376* inv(A) * diag(W),
377* where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
378*
379 DO 90 i = 1, n
380 IF( work( i ).GT.safe2 ) THEN
381 work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
382 ELSE
383 work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
384 END IF
385 90 CONTINUE
386*
387 kase = 0
388 100 CONTINUE
389 CALL slacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
390 $ kase, isave )
391 IF( kase.NE.0 ) THEN
392 IF( kase.EQ.1 ) THEN
393*
394* Multiply by diag(W)*inv(A**T).
395*
396 CALL spotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
397 DO 110 i = 1, n
398 work( n+i ) = work( i )*work( n+i )
399 110 CONTINUE
400 ELSE IF( kase.EQ.2 ) THEN
401*
402* Multiply by inv(A)*diag(W).
403*
404 DO 120 i = 1, n
405 work( n+i ) = work( i )*work( n+i )
406 120 CONTINUE
407 CALL spotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
408 END IF
409 GO TO 100
410 END IF
411*
412* Normalize error.
413*
414 lstres = zero
415 DO 130 i = 1, n
416 lstres = max( lstres, abs( x( i, j ) ) )
417 130 CONTINUE
418 IF( lstres.NE.zero )
419 $ ferr( j ) = ferr( j ) / lstres
420*
421 140 CONTINUE
422*
423 RETURN
424*
425* End of SPORFS
426*

◆ sporfsx()

subroutine sporfsx ( character uplo,
character equed,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldaf, * ) af,
integer ldaf,
real, dimension( * ) s,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldx, * ) x,
integer ldx,
real rcond,
real, dimension( * ) berr,
integer n_err_bnds,
real, dimension( nrhs, * ) err_bnds_norm,
real, dimension( nrhs, * ) err_bnds_comp,
integer nparams,
real, dimension( * ) params,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SPORFSX

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

Purpose:
!>
!>    SPORFSX improves the computed solution to a system of linear
!>    equations when the coefficient matrix is symmetric positive
!>    definite, and provides error bounds and backward error estimates
!>    for the solution.  In addition to normwise error bound, the code
!>    provides maximum componentwise error bound if possible.  See
!>    comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
!>    error bounds.
!>
!>    The original system of linear equations may have been equilibrated
!>    before calling this routine, as described by arguments EQUED and S
!>    below. In this case, the solution and error bounds returned are
!>    for the original unequilibrated system.
!> 
!>     Some optional parameters are bundled in the PARAMS array.  These
!>     settings determine how refinement is performed, but often the
!>     defaults are acceptable.  If the defaults are acceptable, users
!>     can pass NPARAMS = 0 which prevents the source code from accessing
!>     the PARAMS argument.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>       = 'U':  Upper triangle of A is stored;
!>       = 'L':  Lower triangle of A is stored.
!>
[in]EQUED
!>          EQUED is CHARACTER*1
!>     Specifies the form of equilibration that was done to A
!>     before calling this routine. This is needed to compute
!>     the solution and error bounds correctly.
!>       = 'N':  No equilibration
!>       = 'Y':  Both row and column equilibration, i.e., A has been
!>               replaced by diag(S) * A * diag(S).
!>               The right hand side B has been changed accordingly.
!>
[in]N
!>          N is INTEGER
!>     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]A
!>          A is REAL array, dimension (LDA,N)
!>     The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
!>     upper triangular part of A contains the upper triangular part
!>     of the matrix A, and the strictly lower triangular part of A
!>     is not referenced.  If UPLO = 'L', the leading N-by-N lower
!>     triangular part of A contains the lower triangular part of
!>     the matrix A, and the strictly upper triangular part of A is
!>     not referenced.
!>
[in]LDA
!>          LDA is INTEGER
!>     The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]AF
!>          AF is REAL array, dimension (LDAF,N)
!>     The triangular factor U or L from the Cholesky factorization
!>     A = U**T*U or A = L*L**T, as computed by SPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in,out]S
!>          S is REAL array, dimension (N)
!>     The scale factors for A.  If EQUED = 'Y', A is multiplied on
!>     the left and right by diag(S).  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.  If S is output, each
!>     element of S is a power of the radix. If S is input, each element
!>     of S should be a power of the radix to ensure a reliable solution
!>     and error estimates. Scaling by powers of the radix does not cause
!>     rounding errors unless the result underflows or overflows.
!>     Rounding errors during scaling lead to refining with a matrix that
!>     is not equivalent to the input matrix, producing error estimates
!>     that may not be reliable.
!>
[in]B
!>          B is REAL array, dimension (LDB,NRHS)
!>     The right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>     The leading dimension of the array B.  LDB >= max(1,N).
!>
[in,out]X
!>          X is REAL array, dimension (LDX,NRHS)
!>     On entry, the solution matrix X, as computed by SGETRS.
!>     On exit, the improved solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>     The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is REAL
!>     Reciprocal scaled condition number.  This is an estimate of the
!>     reciprocal Skeel condition number of the matrix A after
!>     equilibration (if done).  If this is less than the machine
!>     precision (in particular, if it is zero), the matrix is singular
!>     to working precision.  Note that the error may still be small even
!>     if this number is very small and the matrix appears ill-
!>     conditioned.
!>
[out]BERR
!>          BERR is REAL array, dimension (NRHS)
!>     Componentwise relative backward error.  This is 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).
!>
[in]N_ERR_BNDS
!>          N_ERR_BNDS is INTEGER
!>     Number of error bounds to return for each right hand side
!>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
!>     ERR_BNDS_COMP below.
!>
[out]ERR_BNDS_NORM
!>          ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     normwise relative error, which is defined as follows:
!>
!>     Normwise relative error in the ith solution vector:
!>             max_j (abs(XTRUE(j,i) - X(j,i)))
!>            ------------------------------
!>                  max_j abs(X(j,i))
!>
!>     The array is indexed by the type of error information as described
!>     below. There currently are up to three pieces of information
!>     returned.
!>
!>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_NORM(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     componentwise relative error, which is defined as follows:
!>
!>     Componentwise relative error in the ith solution vector:
!>                    abs(XTRUE(j,i) - X(j,i))
!>             max_j ----------------------
!>                         abs(X(j,i))
!>
!>     The array is indexed by the right-hand side i (on which the
!>     componentwise relative error depends), and the type of error
!>     information as described below. There currently are up to three
!>     pieces of information returned for each right-hand side. If
!>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
!>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]NPARAMS
!>          NPARAMS is INTEGER
!>     Specifies the number of parameters set in PARAMS.  If <= 0, the
!>     PARAMS array is never referenced and default values are used.
!>
[in,out]PARAMS
!>          PARAMS is REAL array, dimension NPARAMS
!>     Specifies algorithm parameters.  If an entry is < 0.0, then
!>     that entry will be filled with default value used for that
!>     parameter.  Only positions up to NPARAMS are accessed; defaults
!>     are used for higher-numbered parameters.
!>
!>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
!>            refinement or not.
!>         Default: 1.0
!>            = 0.0:  No refinement is performed, and no error bounds are
!>                    computed.
!>            = 1.0:  Use the double-precision refinement algorithm,
!>                    possibly with doubled-single computations if the
!>                    compilation environment does not support DOUBLE
!>                    PRECISION.
!>              (other values are reserved for future use)
!>
!>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
!>            computations allowed for refinement.
!>         Default: 10
!>         Aggressive: Set to 100 to permit convergence using approximate
!>                     factorizations or factorizations other than LU. If
!>                     the factorization uses a technique other than
!>                     Gaussian elimination, the guarantees in
!>                     err_bnds_norm and err_bnds_comp may no longer be
!>                     trustworthy.
!>
!>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
!>            will attempt to find a solution with small componentwise
!>            relative error in the double-precision algorithm.  Positive
!>            is true, 0.0 is false.
!>         Default: 1.0 (attempt componentwise convergence)
!>
[out]WORK
!>          WORK is REAL array, dimension (4*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit. The solution to every right-hand side is
!>         guaranteed.
!>       < 0:  If INFO = -i, the i-th argument had an illegal value
!>       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
!>         has been completed, but the factor U is exactly singular, so
!>         the solution and error bounds could not be computed. RCOND = 0
!>         is returned.
!>       = N+J: The solution corresponding to the Jth right-hand side is
!>         not guaranteed. The solutions corresponding to other right-
!>         hand sides K with K > J may not be guaranteed as well, but
!>         only the first such right-hand side is reported. If a small
!>         componentwise error is not requested (PARAMS(3) = 0.0) then
!>         the Jth right-hand side is the first with a normwise error
!>         bound that is not guaranteed (the smallest J such
!>         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
!>         the Jth right-hand side is the first with either a normwise or
!>         componentwise error bound that is not guaranteed (the smallest
!>         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
!>         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
!>         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
!>         about all of the right-hand sides check ERR_BNDS_NORM or
!>         ERR_BNDS_COMP.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 390 of file sporfsx.f.

394*
395* -- LAPACK computational routine --
396* -- LAPACK is a software package provided by Univ. of Tennessee, --
397* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
398*
399* .. Scalar Arguments ..
400 CHARACTER UPLO, EQUED
401 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
402 $ N_ERR_BNDS
403 REAL RCOND
404* ..
405* .. Array Arguments ..
406 INTEGER IWORK( * )
407 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
408 $ X( LDX, * ), WORK( * )
409 REAL S( * ), PARAMS( * ), BERR( * ),
410 $ ERR_BNDS_NORM( NRHS, * ),
411 $ ERR_BNDS_COMP( NRHS, * )
412* ..
413*
414* ==================================================================
415*
416* .. Parameters ..
417 REAL ZERO, ONE
418 parameter( zero = 0.0e+0, one = 1.0e+0 )
419 REAL ITREF_DEFAULT, ITHRESH_DEFAULT,
420 $ COMPONENTWISE_DEFAULT
421 REAL RTHRESH_DEFAULT, DZTHRESH_DEFAULT
422 parameter( itref_default = 1.0 )
423 parameter( ithresh_default = 10.0 )
424 parameter( componentwise_default = 1.0 )
425 parameter( rthresh_default = 0.5 )
426 parameter( dzthresh_default = 0.25 )
427 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
428 $ LA_LINRX_CWISE_I
429 parameter( la_linrx_itref_i = 1,
430 $ la_linrx_ithresh_i = 2 )
431 parameter( la_linrx_cwise_i = 3 )
432 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
433 $ LA_LINRX_RCOND_I
434 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
435 parameter( la_linrx_rcond_i = 3 )
436* ..
437* .. Local Scalars ..
438 CHARACTER(1) NORM
439 LOGICAL RCEQU
440 INTEGER J, PREC_TYPE, REF_TYPE
441 INTEGER N_NORMS
442 REAL ANORM, RCOND_TMP
443 REAL ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
444 LOGICAL IGNORE_CWISE
445 INTEGER ITHRESH
446 REAL RTHRESH, UNSTABLE_THRESH
447* ..
448* .. External Subroutines ..
449 EXTERNAL xerbla, spocon, sla_porfsx_extended
450* ..
451* .. Intrinsic Functions ..
452 INTRINSIC max, sqrt
453* ..
454* .. External Functions ..
455 EXTERNAL lsame, ilaprec
456 EXTERNAL slamch, slansy, sla_porcond
457 REAL SLAMCH, SLANSY, SLA_PORCOND
458 LOGICAL LSAME
459 INTEGER ILAPREC
460* ..
461* .. Executable Statements ..
462*
463* Check the input parameters.
464*
465 info = 0
466 ref_type = int( itref_default )
467 IF ( nparams .GE. la_linrx_itref_i ) THEN
468 IF ( params( la_linrx_itref_i ) .LT. 0.0 ) THEN
469 params( la_linrx_itref_i ) = itref_default
470 ELSE
471 ref_type = params( la_linrx_itref_i )
472 END IF
473 END IF
474*
475* Set default parameters.
476*
477 illrcond_thresh = real( n ) * slamch( 'Epsilon' )
478 ithresh = int( ithresh_default )
479 rthresh = rthresh_default
480 unstable_thresh = dzthresh_default
481 ignore_cwise = componentwise_default .EQ. 0.0
482*
483 IF ( nparams.GE.la_linrx_ithresh_i ) THEN
484 IF ( params( la_linrx_ithresh_i ).LT.0.0 ) THEN
485 params( la_linrx_ithresh_i ) = ithresh
486 ELSE
487 ithresh = int( params( la_linrx_ithresh_i ) )
488 END IF
489 END IF
490 IF ( nparams.GE.la_linrx_cwise_i ) THEN
491 IF ( params( la_linrx_cwise_i ).LT.0.0 ) THEN
492 IF ( ignore_cwise ) THEN
493 params( la_linrx_cwise_i ) = 0.0
494 ELSE
495 params( la_linrx_cwise_i ) = 1.0
496 END IF
497 ELSE
498 ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0
499 END IF
500 END IF
501 IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
502 n_norms = 0
503 ELSE IF ( ignore_cwise ) THEN
504 n_norms = 1
505 ELSE
506 n_norms = 2
507 END IF
508*
509 rcequ = lsame( equed, 'Y' )
510*
511* Test input parameters.
512*
513 IF (.NOT.lsame(uplo, 'U') .AND. .NOT.lsame(uplo, 'L')) THEN
514 info = -1
515 ELSE IF( .NOT.rcequ .AND. .NOT.lsame( equed, 'N' ) ) THEN
516 info = -2
517 ELSE IF( n.LT.0 ) THEN
518 info = -3
519 ELSE IF( nrhs.LT.0 ) THEN
520 info = -4
521 ELSE IF( lda.LT.max( 1, n ) ) THEN
522 info = -6
523 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
524 info = -8
525 ELSE IF( ldb.LT.max( 1, n ) ) THEN
526 info = -11
527 ELSE IF( ldx.LT.max( 1, n ) ) THEN
528 info = -13
529 END IF
530 IF( info.NE.0 ) THEN
531 CALL xerbla( 'SPORFSX', -info )
532 RETURN
533 END IF
534*
535* Quick return if possible.
536*
537 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
538 rcond = 1.0
539 DO j = 1, nrhs
540 berr( j ) = 0.0
541 IF ( n_err_bnds .GE. 1 ) THEN
542 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
543 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
544 END IF
545 IF ( n_err_bnds .GE. 2 ) THEN
546 err_bnds_norm( j, la_linrx_err_i ) = 0.0
547 err_bnds_comp( j, la_linrx_err_i ) = 0.0
548 END IF
549 IF ( n_err_bnds .GE. 3 ) THEN
550 err_bnds_norm( j, la_linrx_rcond_i ) = 1.0
551 err_bnds_comp( j, la_linrx_rcond_i ) = 1.0
552 END IF
553 END DO
554 RETURN
555 END IF
556*
557* Default to failure.
558*
559 rcond = 0.0
560 DO j = 1, nrhs
561 berr( j ) = 1.0
562 IF ( n_err_bnds .GE. 1 ) THEN
563 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
564 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
565 END IF
566 IF ( n_err_bnds .GE. 2 ) THEN
567 err_bnds_norm( j, la_linrx_err_i ) = 1.0
568 err_bnds_comp( j, la_linrx_err_i ) = 1.0
569 END IF
570 IF ( n_err_bnds .GE. 3 ) THEN
571 err_bnds_norm( j, la_linrx_rcond_i ) = 0.0
572 err_bnds_comp( j, la_linrx_rcond_i ) = 0.0
573 END IF
574 END DO
575*
576* Compute the norm of A and the reciprocal of the condition
577* number of A.
578*
579 norm = 'I'
580 anorm = slansy( norm, uplo, n, a, lda, work )
581 CALL spocon( uplo, n, af, ldaf, anorm, rcond, work,
582 $ iwork, info )
583*
584* Perform refinement on each right-hand side
585*
586 IF ( ref_type .NE. 0 ) THEN
587
588 prec_type = ilaprec( 'D' )
589
590 CALL sla_porfsx_extended( prec_type, uplo, n,
591 $ nrhs, a, lda, af, ldaf, rcequ, s, b,
592 $ ldb, x, ldx, berr, n_norms, err_bnds_norm, err_bnds_comp,
593 $ work( n+1 ), work( 1 ), work( 2*n+1 ), work( 1 ), rcond,
594 $ ithresh, rthresh, unstable_thresh, ignore_cwise,
595 $ info )
596 END IF
597
598 err_lbnd = max( 10.0, sqrt( real( n ) ) ) * slamch( 'Epsilon' )
599 IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
600*
601* Compute scaled normwise condition number cond(A*C).
602*
603 IF ( rcequ ) THEN
604 rcond_tmp = sla_porcond( uplo, n, a, lda, af, ldaf,
605 $ -1, s, info, work, iwork )
606 ELSE
607 rcond_tmp = sla_porcond( uplo, n, a, lda, af, ldaf,
608 $ 0, s, info, work, iwork )
609 END IF
610 DO j = 1, nrhs
611*
612* Cap the error at 1.0.
613*
614 IF ( n_err_bnds .GE. la_linrx_err_i
615 $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0 )
616 $ err_bnds_norm( j, la_linrx_err_i ) = 1.0
617*
618* Threshold the error (see LAWN).
619*
620 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
621 err_bnds_norm( j, la_linrx_err_i ) = 1.0
622 err_bnds_norm( j, la_linrx_trust_i ) = 0.0
623 IF ( info .LE. n ) info = n + j
624 ELSE IF ( err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd )
625 $ THEN
626 err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
627 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
628 END IF
629*
630* Save the condition number.
631*
632 IF (n_err_bnds .GE. la_linrx_rcond_i) THEN
633 err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
634 END IF
635 END DO
636 END IF
637
638 IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
639*
640* Compute componentwise condition number cond(A*diag(Y(:,J))) for
641* each right-hand side using the current solution as an estimate of
642* the true solution. If the componentwise error estimate is too
643* large, then the solution is a lousy estimate of truth and the
644* estimated RCOND may be too optimistic. To avoid misleading users,
645* the inverse condition number is set to 0.0 when the estimated
646* cwise error is at least CWISE_WRONG.
647*
648 cwise_wrong = sqrt( slamch( 'Epsilon' ) )
649 DO j = 1, nrhs
650 IF (err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
651 $ THEN
652 rcond_tmp = sla_porcond( uplo, n, a, lda, af, ldaf, 1,
653 $ x( 1, j ), info, work, iwork )
654 ELSE
655 rcond_tmp = 0.0
656 END IF
657*
658* Cap the error at 1.0.
659*
660 IF ( n_err_bnds .GE. la_linrx_err_i
661 $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0 )
662 $ err_bnds_comp( j, la_linrx_err_i ) = 1.0
663*
664* Threshold the error (see LAWN).
665*
666 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
667 err_bnds_comp( j, la_linrx_err_i ) = 1.0
668 err_bnds_comp( j, la_linrx_trust_i ) = 0.0
669 IF ( params( la_linrx_cwise_i ) .EQ. 1.0
670 $ .AND. info.LT.n + j ) info = n + j
671 ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
672 $ .LT. err_lbnd ) THEN
673 err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
674 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
675 END IF
676*
677* Save the condition number.
678*
679 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
680 err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
681 END IF
682
683 END DO
684 END IF
685*
686 RETURN
687*
688* End of SPORFSX
689*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
ilaprec
integer function ilaprec(prec)
ILAPREC
Definition ilaprec.f:58
sla_porfsx_extended
subroutine sla_porfsx_extended(prec_type, uplo, n, nrhs, a, lda, af, ldaf, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)
SLA_PORFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric or H...
Definition sla_porfsx_extended.f:387
spocon
subroutine spocon(uplo, n, a, lda, anorm, rcond, work, iwork, info)
SPOCON
Definition spocon.f:121
slansy
real function slansy(norm, uplo, n, a, lda, work)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansy.f:122

◆ spotf2()

subroutine spotf2 ( character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
integer info )

SPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm).

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

Purpose:
!>
!> SPOTF2 computes the Cholesky factorization of a real symmetric
!> positive definite matrix A.
!>
!> The factorization has the form
!>    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 lower triangular.
!>
!> This is the unblocked version of the algorithm, calling Level 2 BLAS.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          Specifies whether the upper or lower triangular part of the
!>          symmetric matrix A is stored.
!>          = 'U':  Upper triangular
!>          = 'L':  Lower triangular
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
!>          n by n upper triangular part of A contains the upper
!>          triangular part of the matrix A, and the strictly lower
!>          triangular part of A is not referenced.  If UPLO = 'L', the
!>          leading n by n lower triangular part of A contains the lower
!>          triangular part of the matrix A, and the strictly upper
!>          triangular part of A is not referenced.
!>
!>          On exit, if INFO = 0, the factor U or L from the Cholesky
!>          factorization A = U**T *U  or A = L*L**T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -k, the k-th argument had an illegal value
!>          > 0: if INFO = k, the leading minor of order k is not
!>               positive definite, and the factorization could not be
!>               completed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 108 of file spotf2.f.

109*
110* -- LAPACK computational routine --
111* -- LAPACK is a software package provided by Univ. of Tennessee, --
112* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
113*
114* .. Scalar Arguments ..
115 CHARACTER UPLO
116 INTEGER INFO, LDA, N
117* ..
118* .. Array Arguments ..
119 REAL A( LDA, * )
120* ..
121*
122* =====================================================================
123*
124* .. Parameters ..
125 REAL ONE, ZERO
126 parameter( one = 1.0e+0, zero = 0.0e+0 )
127* ..
128* .. Local Scalars ..
129 LOGICAL UPPER
130 INTEGER J
131 REAL AJJ
132* ..
133* .. External Functions ..
134 LOGICAL LSAME, SISNAN
135 REAL SDOT
136 EXTERNAL lsame, sdot, sisnan
137* ..
138* .. External Subroutines ..
139 EXTERNAL sgemv, sscal, xerbla
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC max, sqrt
143* ..
144* .. Executable Statements ..
145*
146* Test the input parameters.
147*
148 info = 0
149 upper = lsame( uplo, 'U' )
150 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
151 info = -1
152 ELSE IF( n.LT.0 ) THEN
153 info = -2
154 ELSE IF( lda.LT.max( 1, n ) ) THEN
155 info = -4
156 END IF
157 IF( info.NE.0 ) THEN
158 CALL xerbla( 'SPOTF2', -info )
159 RETURN
160 END IF
161*
162* Quick return if possible
163*
164 IF( n.EQ.0 )
165 $ RETURN
166*
167 IF( upper ) THEN
168*
169* Compute the Cholesky factorization A = U**T *U.
170*
171 DO 10 j = 1, n
172*
173* Compute U(J,J) and test for non-positive-definiteness.
174*
175 ajj = a( j, j ) - sdot( j-1, a( 1, j ), 1, a( 1, j ), 1 )
176 IF( ajj.LE.zero.OR.sisnan( ajj ) ) THEN
177 a( j, j ) = ajj
178 GO TO 30
179 END IF
180 ajj = sqrt( ajj )
181 a( j, j ) = ajj
182*
183* Compute elements J+1:N of row J.
184*
185 IF( j.LT.n ) THEN
186 CALL sgemv( 'Transpose', j-1, n-j, -one, a( 1, j+1 ),
187 $ lda, a( 1, j ), 1, one, a( j, j+1 ), lda )
188 CALL sscal( n-j, one / ajj, a( j, j+1 ), lda )
189 END IF
190 10 CONTINUE
191 ELSE
192*
193* Compute the Cholesky factorization A = L*L**T.
194*
195 DO 20 j = 1, n
196*
197* Compute L(J,J) and test for non-positive-definiteness.
198*
199 ajj = a( j, j ) - sdot( j-1, a( j, 1 ), lda, a( j, 1 ),
200 $ lda )
201 IF( ajj.LE.zero.OR.sisnan( ajj ) ) THEN
202 a( j, j ) = ajj
203 GO TO 30
204 END IF
205 ajj = sqrt( ajj )
206 a( j, j ) = ajj
207*
208* Compute elements J+1:N of column J.
209*
210 IF( j.LT.n ) THEN
211 CALL sgemv( 'No transpose', n-j, j-1, -one, a( j+1, 1 ),
212 $ lda, a( j, 1 ), lda, one, a( j+1, j ), 1 )
213 CALL sscal( n-j, one / ajj, a( j+1, j ), 1 )
214 END IF
215 20 CONTINUE
216 END IF
217 GO TO 40
218*
219 30 CONTINUE
220 info = j
221*
222 40 CONTINUE
223 RETURN
224*
225* End of SPOTF2
226*
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:156

◆ spotrf()

subroutine spotrf ( character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
integer info )

SPOTRF

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

Purpose:
!>
!> SPOTRF computes the Cholesky factorization of a real symmetric
!> positive definite matrix A.
!>
!> The factorization has the form
!>    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 lower triangular.
!>
!> This is the block version of the algorithm, calling Level 3 BLAS.
!>
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 order of the matrix A.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
!>          N-by-N upper triangular part of A contains the upper
!>          triangular part of the matrix A, and the strictly lower
!>          triangular part of A is not referenced.  If UPLO = 'L', the
!>          leading N-by-N lower triangular part of A contains the lower
!>          triangular part of the matrix A, and the strictly upper
!>          triangular part of A is not referenced.
!>
!>          On exit, if INFO = 0, the factor U or L from the Cholesky
!>          factorization A = U**T*U or A = L*L**T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= 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 is not
!>                positive definite, and the factorization could not be
!>                completed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 106 of file spotrf.f.

107*
108* -- LAPACK computational routine --
109* -- LAPACK is a software package provided by Univ. of Tennessee, --
110* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
111*
112* .. Scalar Arguments ..
113 CHARACTER UPLO
114 INTEGER INFO, LDA, N
115* ..
116* .. Array Arguments ..
117 REAL A( LDA, * )
118* ..
119*
120* =====================================================================
121*
122* .. Parameters ..
123 REAL ONE
124 parameter( one = 1.0e+0 )
125* ..
126* .. Local Scalars ..
127 LOGICAL UPPER
128 INTEGER J, JB, NB
129* ..
130* .. External Functions ..
131 LOGICAL LSAME
132 INTEGER ILAENV
133 EXTERNAL lsame, ilaenv
134* ..
135* .. External Subroutines ..
136 EXTERNAL sgemm, spotrf2, ssyrk, strsm, xerbla
137* ..
138* .. Intrinsic Functions ..
139 INTRINSIC max, min
140* ..
141* .. Executable Statements ..
142*
143* Test the input parameters.
144*
145 info = 0
146 upper = lsame( uplo, 'U' )
147 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
148 info = -1
149 ELSE IF( n.LT.0 ) THEN
150 info = -2
151 ELSE IF( lda.LT.max( 1, n ) ) THEN
152 info = -4
153 END IF
154 IF( info.NE.0 ) THEN
155 CALL xerbla( 'SPOTRF', -info )
156 RETURN
157 END IF
158*
159* Quick return if possible
160*
161 IF( n.EQ.0 )
162 $ RETURN
163*
164* Determine the block size for this environment.
165*
166 nb = ilaenv( 1, 'SPOTRF', uplo, n, -1, -1, -1 )
167 IF( nb.LE.1 .OR. nb.GE.n ) THEN
168*
169* Use unblocked code.
170*
171 CALL spotrf2( uplo, n, a, lda, info )
172 ELSE
173*
174* Use blocked code.
175*
176 IF( upper ) THEN
177*
178* Compute the Cholesky factorization A = U**T*U.
179*
180 DO 10 j = 1, n, nb
181*
182* Update and factorize the current diagonal block and test
183* for non-positive-definiteness.
184*
185 jb = min( nb, n-j+1 )
186 CALL ssyrk( 'Upper', 'Transpose', jb, j-1, -one,
187 $ a( 1, j ), lda, one, a( j, j ), lda )
188 CALL spotrf2( 'Upper', jb, a( j, j ), lda, info )
189 IF( info.NE.0 )
190 $ GO TO 30
191 IF( j+jb.LE.n ) THEN
192*
193* Compute the current block row.
194*
195 CALL sgemm( 'Transpose', 'No transpose', jb, n-j-jb+1,
196 $ j-1, -one, a( 1, j ), lda, a( 1, j+jb ),
197 $ lda, one, a( j, j+jb ), lda )
198 CALL strsm( 'Left', 'Upper', 'Transpose', 'Non-unit',
199 $ jb, n-j-jb+1, one, a( j, j ), lda,
200 $ a( j, j+jb ), lda )
201 END IF
202 10 CONTINUE
203*
204 ELSE
205*
206* Compute the Cholesky factorization A = L*L**T.
207*
208 DO 20 j = 1, n, nb
209*
210* Update and factorize the current diagonal block and test
211* for non-positive-definiteness.
212*
213 jb = min( nb, n-j+1 )
214 CALL ssyrk( 'Lower', 'No transpose', jb, j-1, -one,
215 $ a( j, 1 ), lda, one, a( j, j ), lda )
216 CALL spotrf2( 'Lower', jb, a( j, j ), lda, info )
217 IF( info.NE.0 )
218 $ GO TO 30
219 IF( j+jb.LE.n ) THEN
220*
221* Compute the current block column.
222*
223 CALL sgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
224 $ j-1, -one, a( j+jb, 1 ), lda, a( j, 1 ),
225 $ lda, one, a( j+jb, j ), lda )
226 CALL strsm( 'Right', 'Lower', 'Transpose', 'Non-unit',
227 $ n-j-jb+1, jb, one, a( j, j ), lda,
228 $ a( j+jb, j ), lda )
229 END IF
230 20 CONTINUE
231 END IF
232 END IF
233 GO TO 40
234*
235 30 CONTINUE
236 info = info + j - 1
237*
238 40 CONTINUE
239 RETURN
240*
241* End of SPOTRF
242*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
spotrf2
recursive subroutine spotrf2(uplo, n, a, lda, info)
SPOTRF2
Definition spotrf2.f:106
ssyrk
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187

◆ spotrf2()

recursive subroutine spotrf2 ( character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
integer info )

SPOTRF2

Purpose:
!>
!> SPOTRF2 computes the Cholesky factorization of a real symmetric
!> positive definite matrix A using the recursive algorithm.
!>
!> The factorization has the form
!>    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 lower triangular.
!>
!> This is the recursive version of the algorithm. It divides
!> the matrix into four submatrices:
!>
!>        [  A11 | A12  ]  where A11 is n1 by n1 and A22 is n2 by n2
!>    A = [ -----|----- ]  with n1 = n/2
!>        [  A21 | A22  ]       n2 = n-n1
!>
!> The subroutine calls itself to factor A11. Update and scale A21
!> or A12, update A22 then call itself to factor A22.
!>
!>
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 order of the matrix A.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
!>          N-by-N upper triangular part of A contains the upper
!>          triangular part of the matrix A, and the strictly lower
!>          triangular part of A is not referenced.  If UPLO = 'L', the
!>          leading N-by-N lower triangular part of A contains the lower
!>          triangular part of the matrix A, and the strictly upper
!>          triangular part of A is not referenced.
!>
!>          On exit, if INFO = 0, the factor U or L from the Cholesky
!>          factorization A = U**T*U or A = L*L**T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= 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 is not
!>                positive definite, and the factorization could not be
!>                completed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 105 of file spotrf2.f.

106*
107* -- LAPACK computational routine --
108* -- LAPACK is a software package provided by Univ. of Tennessee, --
109* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
110*
111* .. Scalar Arguments ..
112 CHARACTER UPLO
113 INTEGER INFO, LDA, N
114* ..
115* .. Array Arguments ..
116 REAL A( LDA, * )
117* ..
118*
119* =====================================================================
120*
121* .. Parameters ..
122 REAL ONE, ZERO
123 parameter( one = 1.0e+0, zero=0.0e+0 )
124* ..
125* .. Local Scalars ..
126 LOGICAL UPPER
127 INTEGER N1, N2, IINFO
128* ..
129* .. External Functions ..
130 LOGICAL LSAME, SISNAN
131 EXTERNAL lsame, sisnan
132* ..
133* .. External Subroutines ..
134 EXTERNAL ssyrk, strsm, xerbla
135* ..
136* .. Intrinsic Functions ..
137 INTRINSIC max, sqrt
138* ..
139* .. Executable Statements ..
140*
141* Test the input parameters
142*
143 info = 0
144 upper = lsame( uplo, 'U' )
145 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
146 info = -1
147 ELSE IF( n.LT.0 ) THEN
148 info = -2
149 ELSE IF( lda.LT.max( 1, n ) ) THEN
150 info = -4
151 END IF
152 IF( info.NE.0 ) THEN
153 CALL xerbla( 'SPOTRF2', -info )
154 RETURN
155 END IF
156*
157* Quick return if possible
158*
159 IF( n.EQ.0 )
160 $ RETURN
161*
162* N=1 case
163*
164 IF( n.EQ.1 ) THEN
165*
166* Test for non-positive-definiteness
167*
168 IF( a( 1, 1 ).LE.zero.OR.sisnan( a( 1, 1 ) ) ) THEN
169 info = 1
170 RETURN
171 END IF
172*
173* Factor
174*
175 a( 1, 1 ) = sqrt( a( 1, 1 ) )
176*
177* Use recursive code
178*
179 ELSE
180 n1 = n/2
181 n2 = n-n1
182*
183* Factor A11
184*
185 CALL spotrf2( uplo, n1, a( 1, 1 ), lda, iinfo )
186 IF ( iinfo.NE.0 ) THEN
187 info = iinfo
188 RETURN
189 END IF
190*
191* Compute the Cholesky factorization A = U**T*U
192*
193 IF( upper ) THEN
194*
195* Update and scale A12
196*
197 CALL strsm( 'L', 'U', 'T', 'N', n1, n2, one,
198 $ a( 1, 1 ), lda, a( 1, n1+1 ), lda )
199*
200* Update and factor A22
201*
202 CALL ssyrk( uplo, 'T', n2, n1, -one, a( 1, n1+1 ), lda,
203 $ one, a( n1+1, n1+1 ), lda )
204 CALL spotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
205 IF ( iinfo.NE.0 ) THEN
206 info = iinfo + n1
207 RETURN
208 END IF
209*
210* Compute the Cholesky factorization A = L*L**T
211*
212 ELSE
213*
214* Update and scale A21
215*
216 CALL strsm( 'R', 'L', 'T', 'N', n2, n1, one,
217 $ a( 1, 1 ), lda, a( n1+1, 1 ), lda )
218*
219* Update and factor A22
220*
221 CALL ssyrk( uplo, 'N', n2, n1, -one, a( n1+1, 1 ), lda,
222 $ one, a( n1+1, n1+1 ), lda )
223 CALL spotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
224 IF ( iinfo.NE.0 ) THEN
225 info = iinfo + n1
226 RETURN
227 END IF
228 END IF
229 END IF
230 RETURN
231*
232* End of SPOTRF2
233*

◆ spotri()

subroutine spotri ( character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
integer info )

SPOTRI

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

Purpose:
!>
!> SPOTRI computes the inverse of a real symmetric positive definite
!> matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
!> computed by SPOTRF.
!>
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 order of the matrix A.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N)
!>          On entry, the triangular factor U or L from the Cholesky
!>          factorization A = U**T*U or A = L*L**T, as computed by
!>          SPOTRF.
!>          On exit, the upper or lower triangle of the (symmetric)
!>          inverse of A, overwriting the input factor U or L.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= 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 (i,i) element of the factor U or L is
!>                zero, and the inverse could not be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 94 of file spotri.f.

95*
96* -- LAPACK computational routine --
97* -- LAPACK is a software package provided by Univ. of Tennessee, --
98* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
99*
100* .. Scalar Arguments ..
101 CHARACTER UPLO
102 INTEGER INFO, LDA, N
103* ..
104* .. Array Arguments ..
105 REAL A( LDA, * )
106* ..
107*
108* =====================================================================
109*
110* .. External Functions ..
111 LOGICAL LSAME
112 EXTERNAL lsame
113* ..
114* .. External Subroutines ..
115 EXTERNAL slauum, strtri, xerbla
116* ..
117* .. Intrinsic Functions ..
118 INTRINSIC max
119* ..
120* .. Executable Statements ..
121*
122* Test the input parameters.
123*
124 info = 0
125 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
126 info = -1
127 ELSE IF( n.LT.0 ) THEN
128 info = -2
129 ELSE IF( lda.LT.max( 1, n ) ) THEN
130 info = -4
131 END IF
132 IF( info.NE.0 ) THEN
133 CALL xerbla( 'SPOTRI', -info )
134 RETURN
135 END IF
136*
137* Quick return if possible
138*
139 IF( n.EQ.0 )
140 $ RETURN
141*
142* Invert the triangular Cholesky factor U or L.
143*
144 CALL strtri( uplo, 'Non-unit', n, a, lda, info )
145 IF( info.GT.0 )
146 $ RETURN
147*
148* Form inv(U) * inv(U)**T or inv(L)**T * inv(L).
149*
150 CALL slauum( uplo, n, a, lda, info )
151*
152 RETURN
153*
154* End of SPOTRI
155*
slauum
subroutine slauum(uplo, n, a, lda, info)
SLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked...
Definition slauum.f:102
strtri
subroutine strtri(uplo, diag, n, a, lda, info)
STRTRI
Definition strtri.f:109

◆ spotrs()

subroutine spotrs ( character uplo,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
integer info )

SPOTRS

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

Purpose:
!>
!> SPOTRS solves a system of linear equations A*X = B with a symmetric
!> positive definite matrix A using the Cholesky factorization
!> A = U**T*U or A = L*L**T computed by SPOTRF.
!>
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 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]A
!>          A is REAL array, dimension (LDA,N)
!>          The triangular factor U or L from the Cholesky factorization
!>          A = U**T*U or A = L*L**T, as computed by SPOTRF.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          On entry, the right hand side matrix B.
!>          On exit, the 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
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 109 of file spotrs.f.

110*
111* -- LAPACK computational routine --
112* -- LAPACK is a software package provided by Univ. of Tennessee, --
113* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
114*
115* .. Scalar Arguments ..
116 CHARACTER UPLO
117 INTEGER INFO, LDA, LDB, N, NRHS
118* ..
119* .. Array Arguments ..
120 REAL A( LDA, * ), B( LDB, * )
121* ..
122*
123* =====================================================================
124*
125* .. Parameters ..
126 REAL ONE
127 parameter( one = 1.0e+0 )
128* ..
129* .. Local Scalars ..
130 LOGICAL UPPER
131* ..
132* .. External Functions ..
133 LOGICAL LSAME
134 EXTERNAL lsame
135* ..
136* .. External Subroutines ..
137 EXTERNAL strsm, xerbla
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC max
141* ..
142* .. Executable Statements ..
143*
144* Test the input parameters.
145*
146 info = 0
147 upper = lsame( uplo, 'U' )
148 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
149 info = -1
150 ELSE IF( n.LT.0 ) THEN
151 info = -2
152 ELSE IF( nrhs.LT.0 ) THEN
153 info = -3
154 ELSE IF( lda.LT.max( 1, n ) ) THEN
155 info = -5
156 ELSE IF( ldb.LT.max( 1, n ) ) THEN
157 info = -7
158 END IF
159 IF( info.NE.0 ) THEN
160 CALL xerbla( 'SPOTRS', -info )
161 RETURN
162 END IF
163*
164* Quick return if possible
165*
166 IF( n.EQ.0 .OR. nrhs.EQ.0 )
167 $ RETURN
168*
169 IF( upper ) THEN
170*
171* Solve A*X = B where A = U**T *U.
172*
173* Solve U**T *X = B, overwriting B with X.
174*
175 CALL strsm( 'Left', 'Upper', 'Transpose', 'Non-unit', n, nrhs,
176 $ one, a, lda, b, ldb )
177*
178* Solve U*X = B, overwriting B with X.
179*
180 CALL strsm( 'Left', 'Upper', 'No transpose', 'Non-unit', n,
181 $ nrhs, one, a, lda, b, ldb )
182 ELSE
183*
184* Solve A*X = B where A = L*L**T.
185*
186* Solve L*X = B, overwriting B with X.
187*
188 CALL strsm( 'Left', 'Lower', 'No transpose', 'Non-unit', n,
189 $ nrhs, one, a, lda, b, ldb )
190*
191* Solve L**T *X = B, overwriting B with X.
192*
193 CALL strsm( 'Left', 'Lower', 'Transpose', 'Non-unit', n, nrhs,
194 $ one, a, lda, b, ldb )
195 END IF
196*
197 RETURN
198*
199* End of SPOTRS
200*