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

Functions

real function cla_porcond_c (uplo, n, a, lda, af, ldaf, c, capply, info, work, rwork)
 CLA_PORCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for Hermitian positive-definite matrices.
real function cla_porcond_x (uplo, n, a, lda, af, ldaf, x, info, work, rwork)
 CLA_PORCOND_X computes the infinity norm condition number of op(A)*diag(x) for Hermitian positive-definite matrices.
subroutine cla_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)
 CLA_PORFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric or Hermitian positive-definite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.
real function cla_porpvgrw (uplo, ncols, a, lda, af, ldaf, work)
 CLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian positive-definite matrix.
subroutine cpocon (uplo, n, a, lda, anorm, rcond, work, rwork, info)
 CPOCON
subroutine cpoequ (n, a, lda, s, scond, amax, info)
 CPOEQU
subroutine cpoequb (n, a, lda, s, scond, amax, info)
 CPOEQUB
subroutine cporfs (uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx, ferr, berr, work, rwork, info)
 CPORFS
subroutine cporfsx (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, rwork, info)
 CPORFSX
subroutine cpotf2 (uplo, n, a, lda, info)
 CPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm).
subroutine cpotrf (uplo, n, a, lda, info)
 CPOTRF
recursive subroutine cpotrf2 (uplo, n, a, lda, info)
 CPOTRF2
subroutine cpotri (uplo, n, a, lda, info)
 CPOTRI
subroutine cpotrs (uplo, n, nrhs, a, lda, b, ldb, info)
 CPOTRS

Detailed Description

This is the group of complex computational functions for PO matrices

Function Documentation

◆ cla_porcond_c()

real function cla_porcond_c ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldaf, * ) af,
integer ldaf,
real, dimension( * ) c,
logical capply,
integer info,
complex, dimension( * ) work,
real, dimension( * ) rwork )

CLA_PORCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for Hermitian positive-definite matrices.

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

Purpose:
!>
!>    CLA_PORCOND_C Computes the infinity norm condition number of
!>    op(A) * inv(diag(C)) where C is a REAL vector
!>
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 COMPLEX 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 COMPLEX array, dimension (LDAF,N)
!>     The triangular factor U or L from the Cholesky factorization
!>     A = U**H*U or A = L*L**H, as computed by CPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]C
!>          C is REAL array, dimension (N)
!>     The vector C in the formula op(A) * inv(diag(C)).
!>
[in]CAPPLY
!>          CAPPLY is LOGICAL
!>     If .TRUE. then access the vector C in the formula above.
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>     i > 0:  The ith argument is invalid.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N).
!>     Workspace.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 128 of file cla_porcond_c.f.

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

◆ cla_porcond_x()

real function cla_porcond_x ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldaf, * ) af,
integer ldaf,
complex, dimension( * ) x,
integer info,
complex, dimension( * ) work,
real, dimension( * ) rwork )

CLA_PORCOND_X computes the infinity norm condition number of op(A)*diag(x) for Hermitian positive-definite matrices.

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

Purpose:
!>
!>    CLA_PORCOND_X Computes the infinity norm condition number of
!>    op(A) * diag(X) where X is a COMPLEX vector.
!>
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 COMPLEX 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 COMPLEX array, dimension (LDAF,N)
!>     The triangular factor U or L from the Cholesky factorization
!>     A = U**H*U or A = L*L**H, as computed by CPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]X
!>          X is COMPLEX array, dimension (N)
!>     The vector X in the formula op(A) * diag(X).
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>     i > 0:  The ith argument is invalid.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N).
!>     Workspace.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 121 of file cla_porcond_x.f.

123*
124* -- LAPACK computational routine --
125* -- LAPACK is a software package provided by Univ. of Tennessee, --
126* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
127*
128* .. Scalar Arguments ..
129 CHARACTER UPLO
130 INTEGER N, LDA, LDAF, INFO
131* ..
132* .. Array Arguments ..
133 COMPLEX A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
134 REAL RWORK( * )
135* ..
136*
137* =====================================================================
138*
139* .. Local Scalars ..
140 INTEGER KASE, I, J
141 REAL AINVNM, ANORM, TMP
142 LOGICAL UP, UPPER
143 COMPLEX ZDUM
144* ..
145* .. Local Arrays ..
146 INTEGER ISAVE( 3 )
147* ..
148* .. External Functions ..
149 LOGICAL LSAME
150 EXTERNAL lsame
151* ..
152* .. External Subroutines ..
153 EXTERNAL clacn2, cpotrs, xerbla
154* ..
155* .. Intrinsic Functions ..
156 INTRINSIC abs, max, real, aimag
157* ..
158* .. Statement Functions ..
159 REAL CABS1
160* ..
161* .. Statement Function Definitions ..
162 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
163* ..
164* .. Executable Statements ..
165*
166 cla_porcond_x = 0.0e+0
167*
168 info = 0
169 upper = lsame( uplo, 'U' )
170 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
171 info = -1
172 ELSE IF ( n.LT.0 ) THEN
173 info = -2
174 ELSE IF( lda.LT.max( 1, n ) ) THEN
175 info = -4
176 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
177 info = -6
178 END IF
179 IF( info.NE.0 ) THEN
180 CALL xerbla( 'CLA_PORCOND_X', -info )
181 RETURN
182 END IF
183 up = .false.
184 IF ( lsame( uplo, 'U' ) ) up = .true.
185*
186* Compute norm of op(A)*op2(C).
187*
188 anorm = 0.0
189 IF ( up ) THEN
190 DO i = 1, n
191 tmp = 0.0e+0
192 DO j = 1, i
193 tmp = tmp + cabs1( a( j, i ) * x( j ) )
194 END DO
195 DO j = i+1, n
196 tmp = tmp + cabs1( a( i, j ) * x( j ) )
197 END DO
198 rwork( i ) = tmp
199 anorm = max( anorm, tmp )
200 END DO
201 ELSE
202 DO i = 1, n
203 tmp = 0.0e+0
204 DO j = 1, i
205 tmp = tmp + cabs1( a( i, j ) * x( j ) )
206 END DO
207 DO j = i+1, n
208 tmp = tmp + cabs1( a( j, i ) * x( j ) )
209 END DO
210 rwork( i ) = tmp
211 anorm = max( anorm, tmp )
212 END DO
213 END IF
214*
215* Quick return if possible.
216*
217 IF( n.EQ.0 ) THEN
218 cla_porcond_x = 1.0e+0
219 RETURN
220 ELSE IF( anorm .EQ. 0.0e+0 ) THEN
221 RETURN
222 END IF
223*
224* Estimate the norm of inv(op(A)).
225*
226 ainvnm = 0.0e+0
227*
228 kase = 0
229 10 CONTINUE
230 CALL clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
231 IF( kase.NE.0 ) THEN
232 IF( kase.EQ.2 ) THEN
233*
234* Multiply by R.
235*
236 DO i = 1, n
237 work( i ) = work( i ) * rwork( i )
238 END DO
239*
240 IF ( up ) THEN
241 CALL cpotrs( 'U', n, 1, af, ldaf,
242 $ work, n, info )
243 ELSE
244 CALL cpotrs( 'L', n, 1, af, ldaf,
245 $ work, n, info )
246 ENDIF
247*
248* Multiply by inv(X).
249*
250 DO i = 1, n
251 work( i ) = work( i ) / x( i )
252 END DO
253 ELSE
254*
255* Multiply by inv(X**H).
256*
257 DO i = 1, n
258 work( i ) = work( i ) / x( i )
259 END DO
260*
261 IF ( up ) THEN
262 CALL cpotrs( 'U', n, 1, af, ldaf,
263 $ work, n, info )
264 ELSE
265 CALL cpotrs( 'L', n, 1, af, ldaf,
266 $ work, n, info )
267 END IF
268*
269* Multiply by R.
270*
271 DO i = 1, n
272 work( i ) = work( i ) * rwork( i )
273 END DO
274 END IF
275 GO TO 10
276 END IF
277*
278* Compute the estimate of the reciprocal condition number.
279*
280 IF( ainvnm .NE. 0.0e+0 )
281 $ cla_porcond_x = 1.0e+0 / ainvnm
282*
283 RETURN
284*
285* End of CLA_PORCOND_X
286*
cla_porcond_x
real function cla_porcond_x(uplo, n, a, lda, af, ldaf, x, info, work, rwork)
CLA_PORCOND_X computes the infinity norm condition number of op(A)*diag(x) for Hermitian positive-def...
Definition cla_porcond_x.f:123

◆ cla_porfsx_extended()

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

CLA_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 CLA_PORFSX_EXTENDED + dependencies [TGZ] [ZIP] [TXT]

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

Definition at line 380 of file cla_porfsx_extended.f.

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

◆ cla_porpvgrw()

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

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

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

Purpose:
!>
!>
!> CLA_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 COMPLEX 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 COMPLEX 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 CPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[out]WORK
!>          WORK is REAL array, dimension (2*N)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 104 of file cla_porpvgrw.f.

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

◆ cpocon()

subroutine cpocon ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real anorm,
real rcond,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

CPOCON

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

Purpose:
!>
!> CPOCON estimates the reciprocal of the condition number (in the
!> 1-norm) of a complex Hermitian positive definite matrix using the
!> Cholesky factorization A = U**H*U or A = L*L**H computed by CPOTRF.
!>
!> 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 COMPLEX array, dimension (LDA,N)
!>          The triangular factor U or L from the Cholesky factorization
!>          A = U**H*U or A = L*L**H, as computed by CPOTRF.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]ANORM
!>          ANORM is REAL
!>          The 1-norm (or infinity-norm) of the Hermitian matrix A.
!>
[out]RCOND
!>          RCOND is REAL
!>          The reciprocal of the condition number of the matrix A,
!>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
!>          estimate of the 1-norm of inv(A) computed in this routine.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL 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 cpocon.f.

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

◆ cpoequ()

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

CPOEQU

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

Purpose:
!>
!> CPOEQU computes row and column scalings intended to equilibrate a
!> Hermitian 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 COMPLEX array, dimension (LDA,N)
!>          The N-by-N Hermitian positive definite matrix whose scaling
!>          factors are to be computed.  Only the diagonal elements of A
!>          are referenced.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]S
!>          S is REAL array, dimension (N)
!>          If INFO = 0, S contains the scale factors for A.
!>
[out]SCOND
!>          SCOND is REAL
!>          If INFO = 0, S contains the ratio of the smallest S(i) to
!>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
!>          large nor too small, it is not worth scaling by S.
!>
[out]AMAX
!>          AMAX is REAL
!>          Absolute value of largest matrix element.  If AMAX is very
!>          close to overflow or very close to underflow, the matrix
!>          should be scaled.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 112 of file cpoequ.f.

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

◆ cpoequb()

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

CPOEQUB

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

Purpose:
!>
!> CPOEQUB computes row and column scalings intended to equilibrate a
!> Hermitian 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 CPOEQU 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 COMPLEX array, dimension (LDA,N)
!>          The N-by-N Hermitian positive definite matrix whose scaling
!>          factors are to be computed.  Only the diagonal elements of A
!>          are referenced.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]S
!>          S is REAL array, dimension (N)
!>          If INFO = 0, S contains the scale factors for A.
!>
[out]SCOND
!>          SCOND is REAL
!>          If INFO = 0, S contains the ratio of the smallest S(i) to
!>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
!>          large nor too small, it is not worth scaling by S.
!>
[out]AMAX
!>          AMAX is REAL
!>          Absolute value of largest matrix element.  If AMAX is very
!>          close to overflow or very close to underflow, the matrix
!>          should be scaled.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 118 of file cpoequb.f.

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

◆ cporfs()

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

CPORFS

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

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

◆ cporfsx()

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

CPORFSX

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

Purpose:
!>
!>    CPORFSX improves the computed solution to a system of linear
!>    equations when the coefficient matrix is Hermitian 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 COMPLEX array, dimension (LDA,N)
!>     The Hermitian 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 COMPLEX array, dimension (LDAF,N)
!>     The triangular factor U or L from the Cholesky factorization
!>     A = U**H*U or A = L*L**H, as computed by CPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in,out]S
!>          S is REAL array, dimension (N)
!>     The scale factors for A.  If EQUED = 'Y', A is multiplied on
!>     the left and right by diag(S).  S is an input argument if FACT =
!>     'F'; otherwise, S is an output argument.  If FACT = 'F' and EQUED
!>     = 'Y', each element of S must be positive.  If S is output, each
!>     element of S is a power of the radix. If S is input, each element
!>     of S should be a power of the radix to ensure a reliable solution
!>     and error estimates. Scaling by powers of the radix does not cause
!>     rounding errors unless the result underflows or overflows.
!>     Rounding errors during scaling lead to refining with a matrix that
!>     is not equivalent to the input matrix, producing error estimates
!>     that may not be reliable.
!>
[in]B
!>          B is COMPLEX 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 COMPLEX array, dimension (LDX,NRHS)
!>     On entry, the solution matrix X, as computed by SGETRS.
!>     On exit, the improved solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>     The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is REAL
!>     Reciprocal scaled condition number.  This is an estimate of the
!>     reciprocal Skeel condition number of the matrix A after
!>     equilibration (if done).  If this is less than the machine
!>     precision (in particular, if it is zero), the matrix is singular
!>     to working precision.  Note that the error may still be small even
!>     if this number is very small and the matrix appears ill-
!>     conditioned.
!>
[out]BERR
!>          BERR is REAL array, dimension (NRHS)
!>     Componentwise relative backward error.  This is the
!>     componentwise relative backward error of each solution vector X(j)
!>     (i.e., the smallest relative change in any element of A or B that
!>     makes X(j) an exact solution).
!>
[in]N_ERR_BNDS
!>          N_ERR_BNDS is INTEGER
!>     Number of error bounds to return for each right hand side
!>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
!>     ERR_BNDS_COMP below.
!>
[out]ERR_BNDS_NORM
!>          ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     normwise relative error, which is defined as follows:
!>
!>     Normwise relative error in the ith solution vector:
!>             max_j (abs(XTRUE(j,i) - X(j,i)))
!>            ------------------------------
!>                  max_j abs(X(j,i))
!>
!>     The array is indexed by the type of error information as described
!>     below. There currently are up to three pieces of information
!>     returned.
!>
!>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_NORM(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     componentwise relative error, which is defined as follows:
!>
!>     Componentwise relative error in the ith solution vector:
!>                    abs(XTRUE(j,i) - X(j,i))
!>             max_j ----------------------
!>                         abs(X(j,i))
!>
!>     The array is indexed by the right-hand side i (on which the
!>     componentwise relative error depends), and the type of error
!>     information as described below. There currently are up to three
!>     pieces of information returned for each right-hand side. If
!>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
!>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]NPARAMS
!>          NPARAMS is INTEGER
!>     Specifies the number of parameters set in PARAMS.  If <= 0, the
!>     PARAMS array is never referenced and default values are used.
!>
[in,out]PARAMS
!>          PARAMS is REAL array, dimension NPARAMS
!>     Specifies algorithm parameters.  If an entry is < 0.0, then
!>     that entry will be filled with default value used for that
!>     parameter.  Only positions up to NPARAMS are accessed; defaults
!>     are used for higher-numbered parameters.
!>
!>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
!>            refinement or not.
!>         Default: 1.0
!>            = 0.0:  No refinement is performed, and no error bounds are
!>                    computed.
!>            = 1.0:  Use the double-precision refinement algorithm,
!>                    possibly with doubled-single computations if the
!>                    compilation environment does not support DOUBLE
!>                    PRECISION.
!>              (other values are reserved for future use)
!>
!>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
!>            computations allowed for refinement.
!>         Default: 10
!>         Aggressive: Set to 100 to permit convergence using approximate
!>                     factorizations or factorizations other than LU. If
!>                     the factorization uses a technique other than
!>                     Gaussian elimination, the guarantees in
!>                     err_bnds_norm and err_bnds_comp may no longer be
!>                     trustworthy.
!>
!>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
!>            will attempt to find a solution with small componentwise
!>            relative error in the double-precision algorithm.  Positive
!>            is true, 0.0 is false.
!>         Default: 1.0 (attempt componentwise convergence)
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL array, dimension (2*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 389 of file cporfsx.f.

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

◆ cpotf2()

subroutine cpotf2 ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
integer info )

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

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

Purpose:
!>
!> CPOTF2 computes the Cholesky factorization of a complex Hermitian
!> positive definite matrix A.
!>
!> The factorization has the form
!>    A = U**H * U ,  if UPLO = 'U', or
!>    A = L  * L**H,  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
!>          Hermitian 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 COMPLEX array, dimension (LDA,N)
!>          On entry, the Hermitian 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**H *U  or A = L*L**H.
!>
[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 cpotf2.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 COMPLEX A( LDA, * )
120* ..
121*
122* =====================================================================
123*
124* .. Parameters ..
125 REAL ONE, ZERO
126 parameter( one = 1.0e+0, zero = 0.0e+0 )
127 COMPLEX CONE
128 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
129* ..
130* .. Local Scalars ..
131 LOGICAL UPPER
132 INTEGER J
133 REAL AJJ
134* ..
135* .. External Functions ..
136 LOGICAL LSAME, SISNAN
137 COMPLEX CDOTC
138 EXTERNAL lsame, cdotc, sisnan
139* ..
140* .. External Subroutines ..
141 EXTERNAL cgemv, clacgv, csscal, xerbla
142* ..
143* .. Intrinsic Functions ..
144 INTRINSIC max, real, sqrt
145* ..
146* .. Executable Statements ..
147*
148* Test the input parameters.
149*
150 info = 0
151 upper = lsame( uplo, 'U' )
152 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
153 info = -1
154 ELSE IF( n.LT.0 ) THEN
155 info = -2
156 ELSE IF( lda.LT.max( 1, n ) ) THEN
157 info = -4
158 END IF
159 IF( info.NE.0 ) THEN
160 CALL xerbla( 'CPOTF2', -info )
161 RETURN
162 END IF
163*
164* Quick return if possible
165*
166 IF( n.EQ.0 )
167 $ RETURN
168*
169 IF( upper ) THEN
170*
171* Compute the Cholesky factorization A = U**H *U.
172*
173 DO 10 j = 1, n
174*
175* Compute U(J,J) and test for non-positive-definiteness.
176*
177 ajj = real( real( a( j, j ) ) - cdotc( j-1, a( 1, j ), 1,
178 $ a( 1, j ), 1 ) )
179 IF( ajj.LE.zero.OR.sisnan( ajj ) ) THEN
180 a( j, j ) = ajj
181 GO TO 30
182 END IF
183 ajj = sqrt( ajj )
184 a( j, j ) = ajj
185*
186* Compute elements J+1:N of row J.
187*
188 IF( j.LT.n ) THEN
189 CALL clacgv( j-1, a( 1, j ), 1 )
190 CALL cgemv( 'Transpose', j-1, n-j, -cone, a( 1, j+1 ),
191 $ lda, a( 1, j ), 1, cone, a( j, j+1 ), lda )
192 CALL clacgv( j-1, a( 1, j ), 1 )
193 CALL csscal( n-j, one / ajj, a( j, j+1 ), lda )
194 END IF
195 10 CONTINUE
196 ELSE
197*
198* Compute the Cholesky factorization A = L*L**H.
199*
200 DO 20 j = 1, n
201*
202* Compute L(J,J) and test for non-positive-definiteness.
203*
204 ajj = real( real( a( j, j ) ) - cdotc( j-1, a( j, 1 ), lda,
205 $ a( j, 1 ), lda ) )
206 IF( ajj.LE.zero.OR.sisnan( ajj ) ) THEN
207 a( j, j ) = ajj
208 GO TO 30
209 END IF
210 ajj = sqrt( ajj )
211 a( j, j ) = ajj
212*
213* Compute elements J+1:N of column J.
214*
215 IF( j.LT.n ) THEN
216 CALL clacgv( j-1, a( j, 1 ), lda )
217 CALL cgemv( 'No transpose', n-j, j-1, -cone, a( j+1, 1 ),
218 $ lda, a( j, 1 ), lda, cone, a( j+1, j ), 1 )
219 CALL clacgv( j-1, a( j, 1 ), lda )
220 CALL csscal( n-j, one / ajj, a( j+1, j ), 1 )
221 END IF
222 20 CONTINUE
223 END IF
224 GO TO 40
225*
226 30 CONTINUE
227 info = j
228*
229 40 CONTINUE
230 RETURN
231*
232* End of CPOTF2
233*
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
cdotc
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158

◆ cpotrf()

subroutine cpotrf ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
integer info )

CPOTRF

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

Purpose:
!>
!> CPOTRF computes the Cholesky factorization of a complex Hermitian
!> positive definite matrix A.
!>
!> The factorization has the form
!>    A = U**H * U,  if UPLO = 'U', or
!>    A = L  * L**H,  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 COMPLEX array, dimension (LDA,N)
!>          On entry, the Hermitian 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**H*U or A = L*L**H.
!>
[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 cpotrf.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 COMPLEX A( LDA, * )
118* ..
119*
120* =====================================================================
121*
122* .. Parameters ..
123 REAL ONE
124 COMPLEX CONE
125 parameter( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+0 ) )
126* ..
127* .. Local Scalars ..
128 LOGICAL UPPER
129 INTEGER J, JB, NB
130* ..
131* .. External Functions ..
132 LOGICAL LSAME
133 INTEGER ILAENV
134 EXTERNAL lsame, ilaenv
135* ..
136* .. External Subroutines ..
137 EXTERNAL cgemm, cherk, cpotrf2, ctrsm, xerbla
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC max, min
141* ..
142* .. Executable Statements ..
143*
144* Test the input parameters.
145*
146 info = 0
147 upper = lsame( uplo, 'U' )
148 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
149 info = -1
150 ELSE IF( n.LT.0 ) THEN
151 info = -2
152 ELSE IF( lda.LT.max( 1, n ) ) THEN
153 info = -4
154 END IF
155 IF( info.NE.0 ) THEN
156 CALL xerbla( 'CPOTRF', -info )
157 RETURN
158 END IF
159*
160* Quick return if possible
161*
162 IF( n.EQ.0 )
163 $ RETURN
164*
165* Determine the block size for this environment.
166*
167 nb = ilaenv( 1, 'CPOTRF', uplo, n, -1, -1, -1 )
168 IF( nb.LE.1 .OR. nb.GE.n ) THEN
169*
170* Use unblocked code.
171*
172 CALL cpotrf2( uplo, n, a, lda, info )
173 ELSE
174*
175* Use blocked code.
176*
177 IF( upper ) THEN
178*
179* Compute the Cholesky factorization A = U**H *U.
180*
181 DO 10 j = 1, n, nb
182*
183* Update and factorize the current diagonal block and test
184* for non-positive-definiteness.
185*
186 jb = min( nb, n-j+1 )
187 CALL cherk( 'Upper', 'Conjugate transpose', jb, j-1,
188 $ -one, a( 1, j ), lda, one, a( j, j ), lda )
189 CALL cpotrf2( 'Upper', jb, a( j, j ), lda, info )
190 IF( info.NE.0 )
191 $ GO TO 30
192 IF( j+jb.LE.n ) THEN
193*
194* Compute the current block row.
195*
196 CALL cgemm( 'Conjugate transpose', 'No transpose', jb,
197 $ n-j-jb+1, j-1, -cone, a( 1, j ), lda,
198 $ a( 1, j+jb ), lda, cone, a( j, j+jb ),
199 $ lda )
200 CALL ctrsm( 'Left', 'Upper', 'Conjugate transpose',
201 $ 'Non-unit', jb, n-j-jb+1, cone, a( j, j ),
202 $ lda, a( j, j+jb ), lda )
203 END IF
204 10 CONTINUE
205*
206 ELSE
207*
208* Compute the Cholesky factorization A = L*L**H.
209*
210 DO 20 j = 1, n, nb
211*
212* Update and factorize the current diagonal block and test
213* for non-positive-definiteness.
214*
215 jb = min( nb, n-j+1 )
216 CALL cherk( 'Lower', 'No transpose', jb, j-1, -one,
217 $ a( j, 1 ), lda, one, a( j, j ), lda )
218 CALL cpotrf2( 'Lower', jb, a( j, j ), lda, info )
219 IF( info.NE.0 )
220 $ GO TO 30
221 IF( j+jb.LE.n ) THEN
222*
223* Compute the current block column.
224*
225 CALL cgemm( 'No transpose', 'Conjugate transpose',
226 $ n-j-jb+1, jb, j-1, -cone, a( j+jb, 1 ),
227 $ lda, a( j, 1 ), lda, cone, a( j+jb, j ),
228 $ lda )
229 CALL ctrsm( 'Right', 'Lower', 'Conjugate transpose',
230 $ 'Non-unit', n-j-jb+1, jb, cone, a( j, j ),
231 $ lda, a( j+jb, j ), lda )
232 END IF
233 20 CONTINUE
234 END IF
235 END IF
236 GO TO 40
237*
238 30 CONTINUE
239 info = info + j - 1
240*
241 40 CONTINUE
242 RETURN
243*
244* End of CPOTRF
245*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
cpotrf2
recursive subroutine cpotrf2(uplo, n, a, lda, info)
CPOTRF2
Definition cpotrf2.f:106
ctrsm
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180
cherk
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
Definition cherk.f:173
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187

◆ cpotrf2()

recursive subroutine cpotrf2 ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
integer info )

CPOTRF2

Purpose:
!>
!> CPOTRF2 computes the Cholesky factorization of a Hermitian
!> positive definite matrix A using the recursive algorithm.
!>
!> The factorization has the form
!>    A = U**H * U,  if UPLO = 'U', or
!>    A = L  * L**H,  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 COMPLEX array, dimension (LDA,N)
!>          On entry, the Hermitian 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**H*U or A = L*L**H.
!>
[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 cpotrf2.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 COMPLEX A( LDA, * )
117* ..
118*
119* =====================================================================
120*
121* .. Parameters ..
122 REAL ONE, ZERO
123 parameter( one = 1.0e+0, zero = 0.0e+0 )
124 COMPLEX CONE
125 parameter( cone = (1.0e+0, 0.0e+0) )
126* ..
127* .. Local Scalars ..
128 LOGICAL UPPER
129 INTEGER N1, N2, IINFO
130 REAL AJJ
131* ..
132* .. External Functions ..
133 LOGICAL LSAME, SISNAN
134 EXTERNAL lsame, sisnan
135* ..
136* .. External Subroutines ..
137 EXTERNAL cherk, ctrsm, xerbla
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC max, real, sqrt
141* ..
142* .. Executable Statements ..
143*
144* Test the input parameters
145*
146 info = 0
147 upper = lsame( uplo, 'U' )
148 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
149 info = -1
150 ELSE IF( n.LT.0 ) THEN
151 info = -2
152 ELSE IF( lda.LT.max( 1, n ) ) THEN
153 info = -4
154 END IF
155 IF( info.NE.0 ) THEN
156 CALL xerbla( 'CPOTRF2', -info )
157 RETURN
158 END IF
159*
160* Quick return if possible
161*
162 IF( n.EQ.0 )
163 $ RETURN
164*
165* N=1 case
166*
167 IF( n.EQ.1 ) THEN
168*
169* Test for non-positive-definiteness
170*
171 ajj = real( a( 1, 1 ) )
172 IF( ajj.LE.zero.OR.sisnan( ajj ) ) THEN
173 info = 1
174 RETURN
175 END IF
176*
177* Factor
178*
179 a( 1, 1 ) = sqrt( ajj )
180*
181* Use recursive code
182*
183 ELSE
184 n1 = n/2
185 n2 = n-n1
186*
187* Factor A11
188*
189 CALL cpotrf2( uplo, n1, a( 1, 1 ), lda, iinfo )
190 IF ( iinfo.NE.0 ) THEN
191 info = iinfo
192 RETURN
193 END IF
194*
195* Compute the Cholesky factorization A = U**H*U
196*
197 IF( upper ) THEN
198*
199* Update and scale A12
200*
201 CALL ctrsm( 'L', 'U', 'C', 'N', n1, n2, cone,
202 $ a( 1, 1 ), lda, a( 1, n1+1 ), lda )
203*
204* Update and factor A22
205*
206 CALL cherk( uplo, 'C', n2, n1, -one, a( 1, n1+1 ), lda,
207 $ one, a( n1+1, n1+1 ), lda )
208*
209 CALL cpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
210*
211 IF ( iinfo.NE.0 ) THEN
212 info = iinfo + n1
213 RETURN
214 END IF
215*
216* Compute the Cholesky factorization A = L*L**H
217*
218 ELSE
219*
220* Update and scale A21
221*
222 CALL ctrsm( 'R', 'L', 'C', 'N', n2, n1, cone,
223 $ a( 1, 1 ), lda, a( n1+1, 1 ), lda )
224*
225* Update and factor A22
226*
227 CALL cherk( uplo, 'N', n2, n1, -one, a( n1+1, 1 ), lda,
228 $ one, a( n1+1, n1+1 ), lda )
229*
230 CALL cpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
231*
232 IF ( iinfo.NE.0 ) THEN
233 info = iinfo + n1
234 RETURN
235 END IF
236*
237 END IF
238 END IF
239 RETURN
240*
241* End of CPOTRF2
242*

◆ cpotri()

subroutine cpotri ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
integer info )

CPOTRI

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

Purpose:
!>
!> CPOTRI computes the inverse of a complex Hermitian positive definite
!> matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
!> computed by CPOTRF.
!>
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 COMPLEX array, dimension (LDA,N)
!>          On entry, the triangular factor U or L from the Cholesky
!>          factorization A = U**H*U or A = L*L**H, as computed by
!>          CPOTRF.
!>          On exit, the upper or lower triangle of the (Hermitian)
!>          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 cpotri.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 COMPLEX A( LDA, * )
106* ..
107*
108* =====================================================================
109*
110* .. External Functions ..
111 LOGICAL LSAME
112 EXTERNAL lsame
113* ..
114* .. External Subroutines ..
115 EXTERNAL clauum, ctrtri, 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( 'CPOTRI', -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 ctrtri( uplo, 'Non-unit', n, a, lda, info )
145 IF( info.GT.0 )
146 $ RETURN
147*
148* Form inv(U) * inv(U)**H or inv(L)**H * inv(L).
149*
150 CALL clauum( uplo, n, a, lda, info )
151*
152 RETURN
153*
154* End of CPOTRI
155*
clauum
subroutine clauum(uplo, n, a, lda, info)
CLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked...
Definition clauum.f:102
ctrtri
subroutine ctrtri(uplo, diag, n, a, lda, info)
CTRTRI
Definition ctrtri.f:109

◆ cpotrs()

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

CPOTRS

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

Purpose:
!>
!> CPOTRS solves a system of linear equations A*X = B with a Hermitian
!> positive definite matrix A using the Cholesky factorization
!> A = U**H*U or A = L*L**H computed by CPOTRF.
!>
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 COMPLEX array, dimension (LDA,N)
!>          The triangular factor U or L from the Cholesky factorization
!>          A = U**H*U or A = L*L**H, as computed by CPOTRF.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is COMPLEX 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 cpotrs.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 COMPLEX A( LDA, * ), B( LDB, * )
121* ..
122*
123* =====================================================================
124*
125* .. Parameters ..
126 COMPLEX ONE
127 parameter( one = ( 1.0e+0, 0.0e+0 ) )
128* ..
129* .. Local Scalars ..
130 LOGICAL UPPER
131* ..
132* .. External Functions ..
133 LOGICAL LSAME
134 EXTERNAL lsame
135* ..
136* .. External Subroutines ..
137 EXTERNAL ctrsm, 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( 'CPOTRS', -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**H *U.
172*
173* Solve U**H *X = B, overwriting B with X.
174*
175 CALL ctrsm( 'Left', 'Upper', 'Conjugate transpose', 'Non-unit',
176 $ n, nrhs, one, a, lda, b, ldb )
177*
178* Solve U*X = B, overwriting B with X.
179*
180 CALL ctrsm( '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**H.
185*
186* Solve L*X = B, overwriting B with X.
187*
188 CALL ctrsm( 'Left', 'Lower', 'No transpose', 'Non-unit', n,
189 $ nrhs, one, a, lda, b, ldb )
190*
191* Solve L**H *X = B, overwriting B with X.
192*
193 CALL ctrsm( 'Left', 'Lower', 'Conjugate transpose', 'Non-unit',
194 $ n, nrhs, one, a, lda, b, ldb )
195 END IF
196*
197 RETURN
198*
199* End of CPOTRS
200*