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

Functions

double precision function dla_porcond (UPLO, N, A, LDA, AF, LDAF, CMODE, C, INFO, WORK, IWORK)
 DLA_PORCOND estimates the Skeel condition number for a symmetric positive-definite matrix. More...
 
subroutine dla_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)
 DLA_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...
 
double precision function dla_porpvgrw (UPLO, NCOLS, A, LDA, AF, LDAF, WORK)
 DLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian positive-definite matrix. More...
 
subroutine dpocon (UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
 DPOCON More...
 
subroutine dpoequ (N, A, LDA, S, SCOND, AMAX, INFO)
 DPOEQU More...
 
subroutine dpoequb (N, A, LDA, S, SCOND, AMAX, INFO)
 DPOEQUB More...
 
subroutine dporfs (UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
 DPORFS More...
 
subroutine dporfsx (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)
 DPORFSX More...
 
subroutine dpotf2 (UPLO, N, A, LDA, INFO)
 DPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm). More...
 
subroutine dpotrf (UPLO, N, A, LDA, INFO)
 DPOTRF More...
 
recursive subroutine dpotrf2 (UPLO, N, A, LDA, INFO)
 DPOTRF2 More...
 
subroutine dpotri (UPLO, N, A, LDA, INFO)
 DPOTRI More...
 
subroutine dpotrs (UPLO, N, NRHS, A, LDA, B, LDB, INFO)
 DPOTRS More...
 

Detailed Description

This is the group of double computational functions for PO matrices

Function Documentation

double precision function dla_porcond ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
integer  CMODE,
double precision, dimension( * )  C,
integer  INFO,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK 
)

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

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

Purpose:
    DLA_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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPOTRF.
[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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 144 of file dla_porcond.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dla_porfsx_extended ( integer  PREC_TYPE,
character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
logical  COLEQU,
double precision, dimension( * )  C,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldy, * )  Y,
integer  LDY,
double precision, dimension( * )  BERR_OUT,
integer  N_NORMS,
double precision, dimension( nrhs, * )  ERR_BNDS_NORM,
double precision, dimension( nrhs, * )  ERR_BNDS_COMP,
double precision, dimension( * )  RES,
double precision, dimension(*)  AYB,
double precision, dimension( * )  DY,
double precision, dimension( * )  Y_TAIL,
double precision  RCOND,
integer  ITHRESH,
double precision  RTHRESH,
double precision  DZ_UB,
logical  IGNORE_CWISE,
integer  INFO 
)

DLA_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 DLA_PORFSX_EXTENDED + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLA_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 DPORFSX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPOTRF.
[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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension
                    (LDY,NRHS)
     On entry, the solution matrix X, as computed by DPOTRS.
     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 DOUBLE PRECISION 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 DLA_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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
     Workspace to hold the intermediate residual.
[in]AYB
          AYB is DOUBLE PRECISION array, dimension (N)
     Workspace. This can be the same workspace passed for Y_TAIL.
[in]DY
          DY is DOUBLE PRECISION array, dimension (N)
     Workspace to hold the intermediate solution.
[in]Y_TAIL
          Y_TAIL is DOUBLE PRECISION array, dimension (N)
     Workspace to hold the trailing bits of the intermediate solution.
[in]RCOND
          RCOND is DOUBLE PRECISION
     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 DOUBLE PRECISION
     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 DOUBLE PRECISION
     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 DPOTRS had an illegal
             value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 392 of file dla_porfsx_extended.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function dla_porpvgrw ( character*1  UPLO,
integer  NCOLS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
double precision, dimension( * )  WORK 
)

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

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

Purpose:
 DLA_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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPOTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= max(1,N).
[in]WORK
          WORK is DOUBLE PRECISION array, dimension (2*N)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 108 of file dla_porpvgrw.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dpocon ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision  ANORM,
double precision  RCOND,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DPOCON

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

Purpose:
 DPOCON 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 DPOTRF.

 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 DOUBLE PRECISION 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 DPOTRF.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]ANORM
          ANORM is DOUBLE PRECISION
          The 1-norm (or infinity-norm) of the symmetric matrix A.
[out]RCOND
          RCOND is DOUBLE PRECISION
          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 DOUBLE PRECISION 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 dpocon.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  DOUBLE PRECISION anorm, rcond
133 * ..
134 * .. Array Arguments ..
135  INTEGER iwork( * )
136  DOUBLE PRECISION a( lda, * ), work( * )
137 * ..
138 *
139 * =====================================================================
140 *
141 * .. Parameters ..
142  DOUBLE PRECISION one, zero
143  parameter( one = 1.0d+0, zero = 0.0d+0 )
144 * ..
145 * .. Local Scalars ..
146  LOGICAL upper
147  CHARACTER normin
148  INTEGER ix, kase
149  DOUBLE PRECISION ainvnm, scale, scalel, scaleu, smlnum
150 * ..
151 * .. Local Arrays ..
152  INTEGER isave( 3 )
153 * ..
154 * .. External Functions ..
155  LOGICAL lsame
156  INTEGER idamax
157  DOUBLE PRECISION dlamch
158  EXTERNAL lsame, idamax, dlamch
159 * ..
160 * .. External Subroutines ..
161  EXTERNAL dlacn2, dlatrs, drscl, 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( 'DPOCON', -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 = dlamch( 'Safe minimum' )
197 *
198 * Estimate the 1-norm of inv(A).
199 *
200  kase = 0
201  normin = 'N'
202  10 CONTINUE
203  CALL dlacn2( 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 dlatrs( '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 dlatrs( '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 dlatrs( '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 dlatrs( '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 = idamax( n, work, 1 )
236  IF( scale.LT.abs( work( ix ) )*smlnum .OR. scale.EQ.zero )
237  $ GO TO 20
238  CALL drscl( 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 DPOCON
252 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
DLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
Definition: dlatrs.f:240
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:138
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
subroutine drscl(N, SA, SX, INCX)
DRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: drscl.f:86

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dpoequ ( integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  S,
double precision  SCOND,
double precision  AMAX,
integer  INFO 
)

DPOEQU

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

Purpose:
 DPOEQU 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
          If INFO = 0, S contains the scale factors for A.
[out]SCOND
          SCOND is DOUBLE PRECISION
          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 DOUBLE PRECISION
          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 dpoequ.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  DOUBLE PRECISION amax, scond
123 * ..
124 * .. Array Arguments ..
125  DOUBLE PRECISION a( lda, * ), s( * )
126 * ..
127 *
128 * =====================================================================
129 *
130 * .. Parameters ..
131  DOUBLE PRECISION zero, one
132  parameter( zero = 0.0d+0, one = 1.0d+0 )
133 * ..
134 * .. Local Scalars ..
135  INTEGER i
136  DOUBLE PRECISION 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( 'DPOEQU', -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 DPOEQU
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 dpoequb ( integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  S,
double precision  SCOND,
double precision  AMAX,
integer  INFO 
)

DPOEQUB

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

Purpose:
 DPOEQU 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
          If INFO = 0, S contains the scale factors for A.
[out]SCOND
          SCOND is DOUBLE PRECISION
          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 DOUBLE PRECISION
          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 dpoequb.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  DOUBLE PRECISION amax, scond
123 * ..
124 * .. Array Arguments ..
125  DOUBLE PRECISION a( lda, * ), s( * )
126 * ..
127 *
128 * =====================================================================
129 *
130 * .. Parameters ..
131  DOUBLE PRECISION zero, one
132  parameter( zero = 0.0d+0, one = 1.0d+0 )
133 * ..
134 * .. Local Scalars ..
135  INTEGER i
136  DOUBLE PRECISION smin, base, tmp
137 * ..
138 * .. External Functions ..
139  DOUBLE PRECISION dlamch
140  EXTERNAL dlamch
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( 'DPOEQUB', -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 = dlamch( 'B' )
174  tmp = -0.5d+0 / 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 DPOEQUB
214 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dporfs ( character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DPORFS

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

Purpose:
 DPORFS 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPOTRF.
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= max(1,N).
[in]B
          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by DPOTRS.
          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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dporfs.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  DOUBLE PRECISION 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  DOUBLE PRECISION zero
207  parameter( zero = 0.0d+0 )
208  DOUBLE PRECISION one
209  parameter( one = 1.0d+0 )
210  DOUBLE PRECISION two
211  parameter( two = 2.0d+0 )
212  DOUBLE PRECISION three
213  parameter( three = 3.0d+0 )
214 * ..
215 * .. Local Scalars ..
216  LOGICAL upper
217  INTEGER count, i, j, k, kase, nz
218  DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
219 * ..
220 * .. Local Arrays ..
221  INTEGER isave( 3 )
222 * ..
223 * .. External Subroutines ..
224  EXTERNAL daxpy, dcopy, dlacn2, dpotrs, dsymv, xerbla
225 * ..
226 * .. Intrinsic Functions ..
227  INTRINSIC abs, max
228 * ..
229 * .. External Functions ..
230  LOGICAL lsame
231  DOUBLE PRECISION dlamch
232  EXTERNAL lsame, dlamch
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( 'DPORFS', -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 = dlamch( 'Epsilon' )
274  safmin = dlamch( '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 dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
291  CALL dsymv( 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 dpotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
354  CALL daxpy( 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 DLACN2 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 dlacn2( 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 dpotrs( 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 dpotrs( 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 DPORFS
429 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dsymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DSYMV
Definition: dsymv.f:154
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:138
subroutine dpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
DPOTRS
Definition: dpotrs.f:112
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dporfsx ( character  UPLO,
character  EQUED,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
double precision, dimension( * )  S,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision  RCOND,
double precision, dimension( * )  BERR,
integer  N_ERR_BNDS,
double precision, dimension( nrhs, * )  ERR_BNDS_NORM,
double precision, dimension( nrhs, * )  ERR_BNDS_COMP,
integer  NPARAMS,
double precision, dimension( * )  PARAMS,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DPORFSX

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

Purpose:
    DPORFSX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPOTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= max(1,N).
[in,out]S
          S is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
     On entry, the solution matrix X, as computed by DGETRS.
     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 DOUBLE PRECISION
     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 DOUBLE PRECISION 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 DOUBLE PRECISION 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) * dlamch('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) * dlamch('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) * dlamch('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 DOUBLE PRECISION 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) * dlamch('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) * dlamch('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) * dlamch('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 DOUBLE PRECISION 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.0D+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 DOUBLE PRECISION 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 dporfsx.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  DOUBLE PRECISION rcond
407 * ..
408 * .. Array Arguments ..
409  INTEGER iwork( * )
410  DOUBLE PRECISION a( lda, * ), af( ldaf, * ), b( ldb, * ),
411  $ x( ldx, * ), work( * )
412  DOUBLE PRECISION s( * ), params( * ), berr( * ),
413  $ err_bnds_norm( nrhs, * ),
414  $ err_bnds_comp( nrhs, * )
415 * ..
416 *
417 * ==================================================================
418 *
419 * .. Parameters ..
420  DOUBLE PRECISION zero, one
421  parameter( zero = 0.0d+0, one = 1.0d+0 )
422  DOUBLE PRECISION itref_default, ithresh_default
423  DOUBLE PRECISION componentwise_default, rthresh_default
424  DOUBLE PRECISION dzthresh_default
425  parameter( itref_default = 1.0d+0 )
426  parameter( ithresh_default = 10.0d+0 )
427  parameter( componentwise_default = 1.0d+0 )
428  parameter( rthresh_default = 0.5d+0 )
429  parameter( dzthresh_default = 0.25d+0 )
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  DOUBLE PRECISION anorm, rcond_tmp
446  DOUBLE PRECISION illrcond_thresh, err_lbnd, cwise_wrong
447  LOGICAL ignore_cwise
448  INTEGER ithresh
449  DOUBLE PRECISION 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 dlamch, dlansy, dla_porcond
460  DOUBLE PRECISION dlamch, dlansy, dla_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.0d+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 = dble( n ) * dlamch( 'Epsilon' )
482  ithresh = int( ithresh_default )
483  rthresh = rthresh_default
484  unstable_thresh = dzthresh_default
485  ignore_cwise = componentwise_default .EQ. 0.0d+0
486 *
487  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
488  IF ( params( la_linrx_ithresh_i ).LT.0.0d+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.0d+0 ) THEN
496  IF ( ignore_cwise ) THEN
497  params( la_linrx_cwise_i ) = 0.0d+0
498  ELSE
499  params( la_linrx_cwise_i ) = 1.0d+0
500  END IF
501  ELSE
502  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+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( 'DPORFSX', -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.0d+0
543  DO j = 1, nrhs
544  berr( j ) = 0.0d+0
545  IF ( n_err_bnds .GE. 1 ) THEN
546  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
547  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
548  END IF
549  IF ( n_err_bnds .GE. 2 ) THEN
550  err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
551  err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
552  END IF
553  IF ( n_err_bnds .GE. 3 ) THEN
554  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
555  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
556  END IF
557  END DO
558  RETURN
559  END IF
560 *
561 * Default to failure.
562 *
563  rcond = 0.0d+0
564  DO j = 1, nrhs
565  berr( j ) = 1.0d+0
566  IF ( n_err_bnds .GE. 1 ) THEN
567  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
568  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
569  END IF
570  IF ( n_err_bnds .GE. 2 ) THEN
571  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
572  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
573  END IF
574  IF ( n_err_bnds .GE. 3 ) THEN
575  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
576  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+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 = dlansy( norm, uplo, n, a, lda, work )
585  CALL dpocon( 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( 'E' )
593 
594  CALL dla_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.0d+0, sqrt( dble( n ) ) ) * dlamch( '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 = dla_porcond( uplo, n, a, lda, af, ldaf,
609  $ -1, s, info, work, iwork )
610  ELSE
611  rcond_tmp = dla_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.0d+0 )
620  $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+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.0d+0
626  err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+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.0d+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( dlamch( 'Epsilon' ) )
653  DO j = 1, nrhs
654  IF (err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
655  $ THEN
656  rcond_tmp = dla_porcond( uplo, n, a, lda, af, ldaf, 1,
657  $ x( 1, j ), info, work, iwork )
658  ELSE
659  rcond_tmp = 0.0d+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.0d+0 )
666  $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+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.0d+0
672  err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
673  IF ( params( la_linrx_cwise_i ) .EQ. 1.0d+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.0d+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 DPORFSX
693 *
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dla_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)
DLA_PORFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric or H...
double precision function dla_porcond(UPLO, N, A, LDA, AF, LDAF, CMODE, C, INFO, WORK, IWORK)
DLA_PORCOND estimates the Skeel condition number for a symmetric positive-definite matrix...
Definition: dla_porcond.f:144
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:60
subroutine dpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
DPOCON
Definition: dpocon.f:123
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY 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: dlansy.f:124

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dpotf2 ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer  INFO 
)

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

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

Purpose:
 DPOTF2 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 DOUBLE PRECISION 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 dpotf2.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  DOUBLE PRECISION a( lda, * )
123 * ..
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128  DOUBLE PRECISION one, zero
129  parameter( one = 1.0d+0, zero = 0.0d+0 )
130 * ..
131 * .. Local Scalars ..
132  LOGICAL upper
133  INTEGER j
134  DOUBLE PRECISION ajj
135 * ..
136 * .. External Functions ..
137  LOGICAL lsame, disnan
138  DOUBLE PRECISION ddot
139  EXTERNAL lsame, ddot, disnan
140 * ..
141 * .. External Subroutines ..
142  EXTERNAL dgemv, dscal, 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( 'DPOTF2', -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 ) - ddot( j-1, a( 1, j ), 1, a( 1, j ), 1 )
179  IF( ajj.LE.zero.OR.disnan( 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 dgemv( 'Transpose', j-1, n-j, -one, a( 1, j+1 ),
190  $ lda, a( 1, j ), 1, one, a( j, j+1 ), lda )
191  CALL dscal( 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 ) - ddot( j-1, a( j, 1 ), lda, a( j, 1 ),
203  $ lda )
204  IF( ajj.LE.zero.OR.disnan( 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 dgemv( '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 dscal( 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 DPOTF2
229 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dpotrf ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer  INFO 
)

DPOTRF

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

Purpose:
 DPOTRF 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 DOUBLE PRECISION 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 dpotrf.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  DOUBLE PRECISION a( lda, * )
121 * ..
122 *
123 * =====================================================================
124 *
125 * .. Parameters ..
126  DOUBLE PRECISION one
127  parameter( one = 1.0d+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 dgemm, dpotrf2, dsyrk, dtrsm, 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( 'DPOTRF', -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, 'DPOTRF', uplo, n, -1, -1, -1 )
170  IF( nb.LE.1 .OR. nb.GE.n ) THEN
171 *
172 * Use unblocked code.
173 *
174  CALL dpotrf2( 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 dsyrk( 'Upper', 'Transpose', jb, j-1, -one,
190  $ a( 1, j ), lda, one, a( j, j ), lda )
191  CALL dpotrf2( '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 dgemm( '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 dtrsm( '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 dsyrk( 'Lower', 'No transpose', jb, j-1, -one,
218  $ a( j, 1 ), lda, one, a( j, j ), lda )
219  CALL dpotrf2( '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 dgemm( '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 dtrsm( '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 DPOTRF
245 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
recursive subroutine dpotrf2(UPLO, N, A, LDA, INFO)
DPOTRF2
Definition: dpotrf2.f:108
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
Definition: dsyrk.f:171
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:183
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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 dpotrf2 ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer  INFO 
)

DPOTRF2

Purpose:
 DPOTRF2 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 calls 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 DOUBLE PRECISION 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 dpotrf2.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  DOUBLE PRECISION a( lda, * )
120 * ..
121 *
122 * =====================================================================
123 *
124 * .. Parameters ..
125  DOUBLE PRECISION one, zero
126  parameter( one = 1.0d+0, zero = 0.0d+0 )
127 * ..
128 * .. Local Scalars ..
129  LOGICAL upper
130  INTEGER n1, n2, iinfo
131 * ..
132 * .. External Functions ..
133  LOGICAL lsame, disnan
134  EXTERNAL lsame, disnan
135 * ..
136 * .. External Subroutines ..
137  EXTERNAL dsyrk, dtrsm, 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( 'DPOTRF2', -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.disnan( 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 dpotrf2( 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 dtrsm( '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 dsyrk( uplo, 'T', n2, n1, -one, a( 1, n1+1 ), lda,
206  $ one, a( n1+1, n1+1 ), lda )
207  CALL dpotrf2( 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 dtrsm( '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 dsyrk( uplo, 'N', n2, n1, -one, a( n1+1, 1 ), lda,
225  $ one, a( n1+1, n1+1 ), lda )
226  CALL dpotrf2( 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 DPOTRF2
236 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
recursive subroutine dpotrf2(UPLO, N, A, LDA, INFO)
DPOTRF2
Definition: dpotrf2.f:108
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
Definition: dsyrk.f:171
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.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:

subroutine dpotri ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer  INFO 
)

DPOTRI

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

Purpose:
 DPOTRI 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 DPOTRF.
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 DOUBLE PRECISION 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
          DPOTRF.
          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 dpotri.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  DOUBLE PRECISION a( lda, * )
109 * ..
110 *
111 * =====================================================================
112 *
113 * .. External Functions ..
114  LOGICAL lsame
115  EXTERNAL lsame
116 * ..
117 * .. External Subroutines ..
118  EXTERNAL dlauum, dtrtri, 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( 'DPOTRI', -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 dtrtri( 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 dlauum( uplo, n, a, lda, info )
154 *
155  RETURN
156 *
157 * End of DPOTRI
158 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlauum(UPLO, N, A, LDA, INFO)
DLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked...
Definition: dlauum.f:104
subroutine dtrtri(UPLO, DIAG, N, A, LDA, INFO)
DTRTRI
Definition: dtrtri.f:111

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dpotrs ( character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

DPOTRS

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

Purpose:
 DPOTRS 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 DPOTRF.
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 DOUBLE PRECISION 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 DPOTRF.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in,out]B
          B is DOUBLE PRECISION 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 dpotrs.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  DOUBLE PRECISION a( lda, * ), b( ldb, * )
124 * ..
125 *
126 * =====================================================================
127 *
128 * .. Parameters ..
129  DOUBLE PRECISION one
130  parameter( one = 1.0d+0 )
131 * ..
132 * .. Local Scalars ..
133  LOGICAL upper
134 * ..
135 * .. External Functions ..
136  LOGICAL lsame
137  EXTERNAL lsame
138 * ..
139 * .. External Subroutines ..
140  EXTERNAL dtrsm, 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( 'DPOTRS', -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 dtrsm( '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 dtrsm( '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 dtrsm( '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 dtrsm( 'Left', 'Lower', 'Transpose', 'Non-unit', n, nrhs,
197  $ one, a, lda, b, ldb )
198  END IF
199 *
200  RETURN
201 *
202 * End of DPOTRS
203 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.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: