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

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.
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.
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.
subroutine dpocon (uplo, n, a, lda, anorm, rcond, work, iwork, info)
 DPOCON
subroutine dpoequ (n, a, lda, s, scond, amax, info)
 DPOEQU
subroutine dpoequb (n, a, lda, s, scond, amax, info)
 DPOEQUB
subroutine dporfs (uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx, ferr, berr, work, iwork, info)
 DPORFS
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
subroutine dpotf2 (uplo, n, a, lda, info)
 DPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm).
subroutine dpotrf (uplo, n, a, lda, info)
 DPOTRF
recursive subroutine dpotrf2 (uplo, n, a, lda, info)
 DPOTRF2
subroutine dpotri (uplo, n, a, lda, info)
 DPOTRI
subroutine dpotrs (uplo, n, nrhs, a, lda, b, ldb, info)
 DPOTRS

Detailed Description

This is the group of double computational functions for PO matrices

Function Documentation

◆ dla_porcond()

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.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (3*N).
!>     Workspace.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 139 of file dla_porcond.f.

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

◆ dla_porfsx_extended()

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 responsible for setting the second fields of
!> ERR_BNDS_NORM and ERR_BNDS_COMP.
!>
Parameters
[in]PREC_TYPE
!>          PREC_TYPE is INTEGER
!>     Specifies the intermediate precision to be used in refinement.
!>     The value is defined by ILAPREC(P) where P is a CHARACTER and P
!>          = 'S':  Single
!>          = 'D':  Double
!>          = 'I':  Indigenous
!>          = 'X' or 'E':  Extra
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>       = 'U':  Upper triangle of A is stored;
!>       = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>     The number of linear equations, i.e., the order of the
!>     matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>     The number of right-hand-sides, i.e., the number of columns of the
!>     matrix B.
!>
[in]A
!>          A is 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  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     This subroutine is only responsible for setting the second field
!>     above.
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in,out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is 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 < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     This subroutine is only responsible for setting the second field
!>     above.
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]RES
!>          RES is 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 define as the relative
!>     change in each component being less than DZ_UB. The default value
!>     is 0.25, requiring the first bit to be stable. See LAWN 165 for
!>     more details.
!>
[in]IGNORE_CWISE
!>          IGNORE_CWISE is LOGICAL
!>     If .TRUE. then ignore componentwise convergence. Default value
!>     is .FALSE..
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>       < 0:  if INFO = -i, the ith argument to DPOTRS had an illegal
!>             value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 380 of file dla_porfsx_extended.f.

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

◆ dla_porpvgrw()

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  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).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (2*N)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 104 of file dla_porpvgrw.f.

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

◆ dpocon()

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.

Definition at line 119 of file dpocon.f.

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

◆ dpoequ()

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.

Definition at line 111 of file dpoequ.f.

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

◆ dpoequb()

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:
!>
!> DPOEQUB computes row and column scalings intended to equilibrate a
!> symmetric positive definite matrix A and reduce its condition number
!> (with respect to the two-norm).  S contains the scale factors,
!> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
!> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
!> choice of S puts the condition number of B within a factor N of the
!> smallest possible condition number over all possible diagonal
!> scalings.
!>
!> This routine differs from DPOEQU by restricting the scaling factors
!> to a power of the radix.  Barring over- and underflow, scaling by
!> these factors introduces no additional rounding errors.  However, the
!> scaled diagonal entries are no longer approximately 1 but lie
!> between sqrt(radix) and 1/sqrt(radix).
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]A
!>          A is 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.

Definition at line 117 of file dpoequb.f.

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

◆ dporfs()

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.

Definition at line 181 of file dporfs.f.

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

◆ dporfsx()

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 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  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * 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 . 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 < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * 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 . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]NPARAMS
!>          NPARAMS is INTEGER
!>     Specifies the number of parameters set in PARAMS.  If <= 0, the
!>     PARAMS array is never referenced and default values are used.
!>
[in,out]PARAMS
!>          PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
!>     Specifies algorithm parameters.  If an entry is < 0.0, then
!>     that entry will be filled with default value used for that
!>     parameter.  Only positions up to NPARAMS are accessed; defaults
!>     are used for higher-numbered parameters.
!>
!>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
!>            refinement or not.
!>         Default: 1.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.

Definition at line 390 of file dporfsx.f.

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

◆ dpotf2()

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.

Definition at line 108 of file dpotf2.f.

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

◆ dpotrf()

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.

Definition at line 106 of file dpotrf.f.

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

◆ dpotrf2()

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.

Definition at line 105 of file dpotrf2.f.

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

◆ dpotri()

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.

Definition at line 94 of file dpotri.f.

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

◆ dpotrs()

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.

Definition at line 109 of file dpotrs.f.

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