LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for real:

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. More...
 
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. More...
 
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. More...
 
subroutine spocon (UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
 SPOCON More...
 
subroutine spoequ (N, A, LDA, S, SCOND, AMAX, INFO)
 SPOEQU More...
 
subroutine spoequb (N, A, LDA, S, SCOND, AMAX, INFO)
 SPOEQUB More...
 
subroutine sporfs (UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
 SPORFS More...
 
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 More...
 
subroutine spotf2 (UPLO, N, A, LDA, INFO)
 SPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm). More...
 
subroutine spotrf (UPLO, N, A, LDA, INFO)
 SPOTRF More...
 
recursive subroutine spotrf2 (UPLO, N, A, LDA, INFO)
 SPOTRF2 More...
 
subroutine spotri (UPLO, N, A, LDA, INFO)
 SPOTRI More...
 
subroutine spotrs (UPLO, N, NRHS, A, LDA, B, LDB, INFO)
 SPOTRS More...
 

Detailed Description

This is the group of real computational functions for PO matrices

Function Documentation

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.
[in]WORK
          WORK is REAL array, dimension (3*N).
     Workspace.
[in]IWORK
          IWORK is INTEGER array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 142 of file sla_porcond.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 resonsible 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', '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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" 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 "guaranteed". 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 .LT. 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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" 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 "guaranteed". 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 definte 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.
Date
September 2012

Definition at line 389 of file sla_porfsx_extended.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 "max absolute element" 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).
[in]WORK
          WORK is REAL array, dimension (2*N)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 106 of file sla_porpvgrw.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 123 of file spocon.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 114 of file spoequ.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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:
 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.
Date
November 2011

Definition at line 114 of file spoequb.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 185 of file sporfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 row 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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" 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 "guaranteed". 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 .LT. 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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" 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 "guaranteed". 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 .LE. 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 .LT. 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.
Date
April 2012

Definition at line 396 of file sporfsx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
September 2012

Definition at line 111 of file spotf2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2015

Definition at line 109 of file spotrf.f.

109 *
110 * -- LAPACK computational routine (version 3.6.0) --
111 * -- LAPACK is a software package provided by Univ. of Tennessee, --
112 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
113 * November 2015
114 *
115 * .. Scalar Arguments ..
116  CHARACTER uplo
117  INTEGER info, lda, n
118 * ..
119 * .. Array Arguments ..
120  REAL a( lda, * )
121 * ..
122 *
123 * =====================================================================
124 *
125 * .. Parameters ..
126  REAL one
127  parameter( one = 1.0e+0 )
128 * ..
129 * .. Local Scalars ..
130  LOGICAL upper
131  INTEGER j, jb, nb
132 * ..
133 * .. External Functions ..
134  LOGICAL lsame
135  INTEGER ilaenv
136  EXTERNAL lsame, ilaenv
137 * ..
138 * .. External Subroutines ..
139  EXTERNAL sgemm, spotrf2, ssyrk, strsm, xerbla
140 * ..
141 * .. Intrinsic Functions ..
142  INTRINSIC max, min
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( 'SPOTRF', -info )
159  RETURN
160  END IF
161 *
162 * Quick return if possible
163 *
164  IF( n.EQ.0 )
165  $ RETURN
166 *
167 * Determine the block size for this environment.
168 *
169  nb = ilaenv( 1, 'SPOTRF', uplo, n, -1, -1, -1 )
170  IF( nb.LE.1 .OR. nb.GE.n ) THEN
171 *
172 * Use unblocked code.
173 *
174  CALL spotrf2( uplo, n, a, lda, info )
175  ELSE
176 *
177 * Use blocked code.
178 *
179  IF( upper ) THEN
180 *
181 * Compute the Cholesky factorization A = U**T*U.
182 *
183  DO 10 j = 1, n, nb
184 *
185 * Update and factorize the current diagonal block and test
186 * for non-positive-definiteness.
187 *
188  jb = min( nb, n-j+1 )
189  CALL ssyrk( 'Upper', 'Transpose', jb, j-1, -one,
190  $ a( 1, j ), lda, one, a( j, j ), lda )
191  CALL spotrf2( 'Upper', jb, a( j, j ), lda, info )
192  IF( info.NE.0 )
193  $ GO TO 30
194  IF( j+jb.LE.n ) THEN
195 *
196 * Compute the current block row.
197 *
198  CALL sgemm( 'Transpose', 'No transpose', jb, n-j-jb+1,
199  $ j-1, -one, a( 1, j ), lda, a( 1, j+jb ),
200  $ lda, one, a( j, j+jb ), lda )
201  CALL strsm( 'Left', 'Upper', 'Transpose', 'Non-unit',
202  $ jb, n-j-jb+1, one, a( j, j ), lda,
203  $ a( j, j+jb ), lda )
204  END IF
205  10 CONTINUE
206 *
207  ELSE
208 *
209 * Compute the Cholesky factorization A = L*L**T.
210 *
211  DO 20 j = 1, n, nb
212 *
213 * Update and factorize the current diagonal block and test
214 * for non-positive-definiteness.
215 *
216  jb = min( nb, n-j+1 )
217  CALL ssyrk( 'Lower', 'No transpose', jb, j-1, -one,
218  $ a( j, 1 ), lda, one, a( j, j ), lda )
219  CALL spotrf2( 'Lower', jb, a( j, j ), lda, info )
220  IF( info.NE.0 )
221  $ GO TO 30
222  IF( j+jb.LE.n ) THEN
223 *
224 * Compute the current block column.
225 *
226  CALL sgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
227  $ j-1, -one, a( j+jb, 1 ), lda, a( j, 1 ),
228  $ lda, one, a( j+jb, j ), lda )
229  CALL strsm( 'Right', 'Lower', 'Transpose', 'Non-unit',
230  $ n-j-jb+1, jb, one, a( j, j ), lda,
231  $ a( j+jb, j ), lda )
232  END IF
233  20 CONTINUE
234  END IF
235  END IF
236  GO TO 40
237 *
238  30 CONTINUE
239  info = info + j - 1
240 *
241  40 CONTINUE
242  RETURN
243 *
244 * End of SPOTRF
245 *
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:171
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183
recursive subroutine spotrf2(UPLO, N, A, LDA, INFO)
SPOTRF2
Definition: spotrf2.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2015

Definition at line 108 of file spotrf2.f.

108 *
109 * -- LAPACK computational routine (version 3.6.0) --
110 * -- LAPACK is a software package provided by Univ. of Tennessee, --
111 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
112 * November 2015
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 n1, n2, iinfo
131 * ..
132 * .. External Functions ..
133  LOGICAL lsame, sisnan
134  EXTERNAL lsame, sisnan
135 * ..
136 * .. External Subroutines ..
137  EXTERNAL sgemm, ssyrk, xerbla
138 * ..
139 * .. Intrinsic Functions ..
140  INTRINSIC max, sqrt
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( lda.LT.max( 1, n ) ) THEN
153  info = -4
154  END IF
155  IF( info.NE.0 ) THEN
156  CALL xerbla( 'SPOTRF2', -info )
157  RETURN
158  END IF
159 *
160 * Quick return if possible
161 *
162  IF( n.EQ.0 )
163  $ RETURN
164 *
165 * N=1 case
166 *
167  IF( n.EQ.1 ) THEN
168 *
169 * Test for non-positive-definiteness
170 *
171  IF( a( 1, 1 ).LE.zero.OR.sisnan( a( 1, 1 ) ) ) THEN
172  info = 1
173  RETURN
174  END IF
175 *
176 * Factor
177 *
178  a( 1, 1 ) = sqrt( a( 1, 1 ) )
179 *
180 * Use recursive code
181 *
182  ELSE
183  n1 = n/2
184  n2 = n-n1
185 *
186 * Factor A11
187 *
188  CALL spotrf2( uplo, n1, a( 1, 1 ), lda, iinfo )
189  IF ( iinfo.NE.0 ) THEN
190  info = iinfo
191  RETURN
192  END IF
193 *
194 * Compute the Cholesky factorization A = U**T*U
195 *
196  IF( upper ) THEN
197 *
198 * Update and scale A12
199 *
200  CALL strsm( 'L', 'U', 'T', 'N', n1, n2, one,
201  $ a( 1, 1 ), lda, a( 1, n1+1 ), lda )
202 *
203 * Update and factor A22
204 *
205  CALL ssyrk( uplo, 'T', n2, n1, -one, a( 1, n1+1 ), lda,
206  $ one, a( n1+1, n1+1 ), lda )
207  CALL spotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
208  IF ( iinfo.NE.0 ) THEN
209  info = iinfo + n1
210  RETURN
211  END IF
212 *
213 * Compute the Cholesky factorization A = L*L**T
214 *
215  ELSE
216 *
217 * Update and scale A21
218 *
219  CALL strsm( 'R', 'L', 'T', 'N', n2, n1, one,
220  $ a( 1, 1 ), lda, a( n1+1, 1 ), lda )
221 *
222 * Update and factor A22
223 *
224  CALL ssyrk( uplo, 'N', n2, n1, -one, a( n1+1, 1 ), lda,
225  $ one, a( n1+1, n1+1 ), lda )
226  CALL spotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
227  IF ( iinfo.NE.0 ) THEN
228  info = iinfo + n1
229  RETURN
230  END IF
231  END IF
232  END IF
233  RETURN
234 *
235 * End of SPOTRF2
236 *
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:171
logical function sisnan(SIN)
SISNAN tests input for NaN.
Definition: sisnan.f:61
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183
recursive subroutine spotrf2(UPLO, N, A, LDA, INFO)
SPOTRF2
Definition: spotrf2.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 97 of file spotri.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 112 of file spotrs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: