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

Functions

double precision function zla_porcond_c (uplo, n, a, lda, af, ldaf, c, capply, info, work, rwork)
 ZLA_PORCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for Hermitian positive-definite matrices.
double precision function zla_porcond_x (uplo, n, a, lda, af, ldaf, x, info, work, rwork)
 ZLA_PORCOND_X computes the infinity norm condition number of op(A)*diag(x) for Hermitian positive-definite matrices.
subroutine zla_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)
 ZLA_PORFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric or Hermitian positive-definite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.
double precision function zla_porpvgrw (uplo, ncols, a, lda, af, ldaf, work)
 ZLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian positive-definite matrix.
subroutine zpocon (uplo, n, a, lda, anorm, rcond, work, rwork, info)
 ZPOCON
subroutine zpoequ (n, a, lda, s, scond, amax, info)
 ZPOEQU
subroutine zpoequb (n, a, lda, s, scond, amax, info)
 ZPOEQUB
subroutine zporfs (uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx, ferr, berr, work, rwork, info)
 ZPORFS
subroutine zporfsx (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)
 ZPORFSX
subroutine zpotf2 (uplo, n, a, lda, info)
 ZPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm).
subroutine zpotrf (uplo, n, a, lda, info)
 ZPOTRF
recursive subroutine zpotrf2 (uplo, n, a, lda, info)
 ZPOTRF2
subroutine zpotri (uplo, n, a, lda, info)
 ZPOTRI
subroutine zpotrs (uplo, n, nrhs, a, lda, b, ldb, info)
 ZPOTRS

Detailed Description

This is the group of complex16 computational functions for PO matrices

Function Documentation

◆ zla_porcond_c()

double precision function zla_porcond_c ( character uplo,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldaf, * ) af,
integer ldaf,
double precision, dimension( * ) c,
logical capply,
integer info,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork )

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

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

Purpose:
!>
!>    ZLA_PORCOND_C Computes the infinity norm condition number of
!>    op(A) * inv(diag(C)) where C is a DOUBLE PRECISION 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*16 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*16 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 ZPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]C
!>          C is DOUBLE PRECISION 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*16 array, dimension (2*N).
!>     Workspace.
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 129 of file zla_porcond_c.f.

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

◆ zla_porcond_x()

double precision function zla_porcond_x ( character uplo,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldaf, * ) af,
integer ldaf,
complex*16, dimension( * ) x,
integer info,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork )

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

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

Purpose:
!>
!>    ZLA_PORCOND_X Computes the infinity norm condition number of
!>    op(A) * diag(X) where X is a COMPLEX*16 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*16 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*16 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 ZPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]X
!>          X is COMPLEX*16 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*16 array, dimension (2*N).
!>     Workspace.
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 122 of file zla_porcond_x.f.

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

◆ zla_porfsx_extended()

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

ZLA_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 ZLA_PORFSX_EXTENDED + dependencies [TGZ] [ZIP] [TXT]

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

Definition at line 380 of file zla_porfsx_extended.f.

387*
388* -- LAPACK computational routine --
389* -- LAPACK is a software package provided by Univ. of Tennessee, --
390* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
391*
392* .. Scalar Arguments ..
393 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
394 $ N_NORMS, ITHRESH
395 CHARACTER UPLO
396 LOGICAL COLEQU, IGNORE_CWISE
397 DOUBLE PRECISION RTHRESH, DZ_UB
398* ..
399* .. Array Arguments ..
400 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
401 $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
402 DOUBLE PRECISION C( * ), AYB( * ), RCOND, BERR_OUT( * ),
403 $ ERR_BNDS_NORM( NRHS, * ),
404 $ ERR_BNDS_COMP( NRHS, * )
405* ..
406*
407* =====================================================================
408*
409* .. Local Scalars ..
410 INTEGER UPLO2, CNT, I, J, X_STATE, Z_STATE,
411 $ Y_PREC_STATE
412 DOUBLE PRECISION 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*16 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 zaxpy, zcopy, zpotrs, zhemv, blas_zhemv_x,
452 $ blas_zhemv2_x, zla_heamv, zla_wwaddw,
453 $ zla_lin_berr, dlamch
454 DOUBLE PRECISION DLAMCH
455* ..
456* .. Intrinsic Functions ..
457 INTRINSIC abs, dble, dimag, max, min
458* ..
459* .. Statement Functions ..
460 DOUBLE PRECISION CABS1
461* ..
462* .. Statement Function Definitions ..
463 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
464* ..
465* .. Executable Statements ..
466*
467 IF (info.NE.0) RETURN
468 eps = dlamch( 'Epsilon' )
469 hugeval = dlamch( 'Overflow' )
470* Force HUGEVAL to Inf
471 hugeval = hugeval * hugeval
472* Using HUGEVAL may lead to spurious underflows.
473 incr_thresh = dble(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.0d+0
486 END DO
487 END IF
488
489 dxrat = 0.0d+0
490 dxratmax = 0.0d+0
491 dzrat = 0.0d+0
492 dzratmax = 0.0d+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 zcopy( n, b( 1, j ), 1, res, 1 )
510 IF (y_prec_state .EQ. base_residual) THEN
511 CALL zhemv(uplo, n, dcmplx(-1.0d+0), a, lda, y(1,j), 1,
512 $ dcmplx(1.0d+0), res, 1)
513 ELSE IF (y_prec_state .EQ. extra_residual) THEN
514 CALL blas_zhemv_x(uplo2, n, dcmplx(-1.0d+0), a, lda,
515 $ y( 1, j ), 1, dcmplx(1.0d+0), res, 1, prec_type)
516 ELSE
517 CALL blas_zhemv2_x(uplo2, n, dcmplx(-1.0d+0), a, lda,
518 $ y(1, j), y_tail, 1, dcmplx(1.0d+0), res, 1,
519 $ prec_type)
520 END IF
521
522! XXX: RES is no longer needed.
523 CALL zcopy( n, res, 1, dy, 1 )
524 CALL zpotrs( uplo, n, 1, af, ldaf, dy, n, info)
525*
526* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
527*
528 normx = 0.0d+0
529 normy = 0.0d+0
530 normdx = 0.0d+0
531 dz_z = 0.0d+0
532 ymin = hugeval
533
534 DO i = 1, n
535 yk = cabs1(y(i, j))
536 dyk = cabs1(dy(i))
537
538 IF (yk .NE. 0.0d+0) THEN
539 dz_z = max( dz_z, dyk / yk )
540 ELSE IF (dyk .NE. 0.0d+0) THEN
541 dz_z = hugeval
542 END IF
543
544 ymin = min( ymin, yk )
545
546 normy = max( normy, yk )
547
548 IF ( colequ ) THEN
549 normx = max(normx, yk * c(i))
550 normdx = max(normdx, dyk * c(i))
551 ELSE
552 normx = normy
553 normdx = max(normdx, dyk)
554 END IF
555 END DO
556
557 IF (normx .NE. 0.0d+0) THEN
558 dx_x = normdx / normx
559 ELSE IF (normdx .EQ. 0.0d+0) THEN
560 dx_x = 0.0d+0
561 ELSE
562 dx_x = hugeval
563 END IF
564
565 dxrat = normdx / prevnormdx
566 dzrat = dz_z / prev_dz_z
567*
568* Check termination criteria.
569*
570 IF (ymin*rcond .LT. incr_thresh*normy
571 $ .AND. y_prec_state .LT. extra_y)
572 $ incr_prec = .true.
573
574 IF (x_state .EQ. noprog_state .AND. dxrat .LE. rthresh)
575 $ x_state = working_state
576 IF (x_state .EQ. working_state) THEN
577 IF (dx_x .LE. eps) THEN
578 x_state = conv_state
579 ELSE IF (dxrat .GT. rthresh) THEN
580 IF (y_prec_state .NE. extra_y) THEN
581 incr_prec = .true.
582 ELSE
583 x_state = noprog_state
584 END IF
585 ELSE
586 IF (dxrat .GT. dxratmax) dxratmax = dxrat
587 END IF
588 IF (x_state .GT. working_state) final_dx_x = dx_x
589 END IF
590
591 IF (z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub)
592 $ z_state = working_state
593 IF (z_state .EQ. noprog_state .AND. dzrat .LE. rthresh)
594 $ z_state = working_state
595 IF (z_state .EQ. working_state) THEN
596 IF (dz_z .LE. eps) THEN
597 z_state = conv_state
598 ELSE IF (dz_z .GT. dz_ub) THEN
599 z_state = unstable_state
600 dzratmax = 0.0d+0
601 final_dz_z = hugeval
602 ELSE IF (dzrat .GT. rthresh) THEN
603 IF (y_prec_state .NE. extra_y) THEN
604 incr_prec = .true.
605 ELSE
606 z_state = noprog_state
607 END IF
608 ELSE
609 IF (dzrat .GT. dzratmax) dzratmax = dzrat
610 END IF
611 IF (z_state .GT. working_state) final_dz_z = dz_z
612 END IF
613
614 IF ( x_state.NE.working_state.AND.
615 $ (ignore_cwise.OR.z_state.NE.working_state) )
616 $ GOTO 666
617
618 IF (incr_prec) THEN
619 incr_prec = .false.
620 y_prec_state = y_prec_state + 1
621 DO i = 1, n
622 y_tail( i ) = 0.0d+0
623 END DO
624 END IF
625
626 prevnormdx = normdx
627 prev_dz_z = dz_z
628*
629* Update soluton.
630*
631 IF (y_prec_state .LT. extra_y) THEN
632 CALL zaxpy( n, dcmplx(1.0d+0), dy, 1, y(1,j), 1 )
633 ELSE
634 CALL zla_wwaddw(n, y(1,j), y_tail, dy)
635 END IF
636
637 END DO
638* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
639 666 CONTINUE
640*
641* Set final_* when cnt hits ithresh.
642*
643 IF (x_state .EQ. working_state) final_dx_x = dx_x
644 IF (z_state .EQ. working_state) final_dz_z = dz_z
645*
646* Compute error bounds.
647*
648 IF (n_norms .GE. 1) THEN
649 err_bnds_norm( j, la_linrx_err_i ) =
650 $ final_dx_x / (1 - dxratmax)
651 END IF
652 IF (n_norms .GE. 2) THEN
653 err_bnds_comp( j, la_linrx_err_i ) =
654 $ final_dz_z / (1 - dzratmax)
655 END IF
656*
657* Compute componentwise relative backward error from formula
658* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
659* where abs(Z) is the componentwise absolute value of the matrix
660* or vector Z.
661*
662* Compute residual RES = B_s - op(A_s) * Y,
663* op(A) = A, A**T, or A**H depending on TRANS (and type).
664*
665 CALL zcopy( n, b( 1, j ), 1, res, 1 )
666 CALL zhemv(uplo, n, dcmplx(-1.0d+0), a, lda, y(1,j), 1,
667 $ dcmplx(1.0d+0), res, 1)
668
669 DO i = 1, n
670 ayb( i ) = cabs1( b( i, j ) )
671 END DO
672*
673* Compute abs(op(A_s))*abs(Y) + abs(B_s).
674*
675 CALL zla_heamv (uplo2, n, 1.0d+0,
676 $ a, lda, y(1, j), 1, 1.0d+0, ayb, 1)
677
678 CALL zla_lin_berr (n, n, 1, res, ayb, berr_out(j))
679*
680* End of loop for each RHS.
681*
682 END DO
683*
684 RETURN
685*
686* End of ZLA_PORFSX_EXTENDED
687*
ilauplo
integer function ilauplo(uplo)
ILAUPLO
Definition ilauplo.f:58
zla_heamv
subroutine zla_heamv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
ZLA_HEAMV computes a matrix-vector product using a Hermitian indefinite matrix to calculate error bou...
Definition zla_heamv.f:178
zla_lin_berr
subroutine zla_lin_berr(n, nz, nrhs, res, ayb, berr)
ZLA_LIN_BERR computes a component-wise relative backward error.
Definition zla_lin_berr.f:101
zla_wwaddw
subroutine zla_wwaddw(n, x, y, w)
ZLA_WWADDW adds a vector into a doubled-single vector.
Definition zla_wwaddw.f:81
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zhemv
subroutine zhemv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
ZHEMV
Definition zhemv.f:154
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
min
#define min(a, b)
Definition macros.h:20

◆ zla_porpvgrw()

double precision function zla_porpvgrw ( character*1 uplo,
integer ncols,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldaf, * ) af,
integer ldaf,
double precision, dimension( * ) work )

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

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

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

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

◆ zpocon()

subroutine zpocon ( character uplo,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
double precision anorm,
double precision rcond,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZPOCON

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

Purpose:
!>
!> ZPOCON 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 ZPOTRF.
!>
!> 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*16 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 ZPOTRF.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]ANORM
!>          ANORM is DOUBLE PRECISION
!>          The 1-norm (or infinity-norm) of the Hermitian matrix A.
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>          The reciprocal of the condition number of the matrix A,
!>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
!>          estimate of the 1-norm of inv(A) computed in this routine.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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 zpocon.f.

121*
122* -- LAPACK computational routine --
123* -- LAPACK is a software package provided by Univ. of Tennessee, --
124* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
125*
126* .. Scalar Arguments ..
127 CHARACTER UPLO
128 INTEGER INFO, LDA, N
129 DOUBLE PRECISION ANORM, RCOND
130* ..
131* .. Array Arguments ..
132 DOUBLE PRECISION RWORK( * )
133 COMPLEX*16 A( LDA, * ), WORK( * )
134* ..
135*
136* =====================================================================
137*
138* .. Parameters ..
139 DOUBLE PRECISION ONE, ZERO
140 parameter( one = 1.0d+0, zero = 0.0d+0 )
141* ..
142* .. Local Scalars ..
143 LOGICAL UPPER
144 CHARACTER NORMIN
145 INTEGER IX, KASE
146 DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
147 COMPLEX*16 ZDUM
148* ..
149* .. Local Arrays ..
150 INTEGER ISAVE( 3 )
151* ..
152* .. External Functions ..
153 LOGICAL LSAME
154 INTEGER IZAMAX
155 DOUBLE PRECISION DLAMCH
156 EXTERNAL lsame, izamax, dlamch
157* ..
158* .. External Subroutines ..
159 EXTERNAL xerbla, zdrscl, zlacn2, zlatrs
160* ..
161* .. Intrinsic Functions ..
162 INTRINSIC abs, dble, dimag, max
163* ..
164* .. Statement Functions ..
165 DOUBLE PRECISION CABS1
166* ..
167* .. Statement Function definitions ..
168 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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( 'ZPOCON', -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 = dlamch( 'Safe minimum' )
201*
202* Estimate the 1-norm of inv(A).
203*
204 kase = 0
205 normin = 'N'
206 10 CONTINUE
207 CALL zlacn2( 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 zlatrs( '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 zlatrs( '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 zlatrs( '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 zlatrs( '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 = izamax( n, work, 1 )
240 IF( scale.LT.cabs1( work( ix ) )*smlnum .OR. scale.EQ.zero )
241 $ GO TO 20
242 CALL zdrscl( 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 ZPOCON
256*
izamax
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
zdrscl
subroutine zdrscl(n, sa, sx, incx)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition zdrscl.f:84
zlatrs
subroutine zlatrs(uplo, trans, diag, normin, n, a, lda, x, scale, cnorm, info)
ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition zlatrs.f:239

◆ zpoequ()

subroutine zpoequ ( integer n,
complex*16, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) s,
double precision scond,
double precision amax,
integer info )

ZPOEQU

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

Purpose:
!>
!> ZPOEQU 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*16 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 DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, S contains the scale factors for A.
!>
[out]SCOND
!>          SCOND is DOUBLE PRECISION
!>          If INFO = 0, S contains the ratio of the smallest S(i) to
!>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
!>          large nor too small, it is not worth scaling by S.
!>
[out]AMAX
!>          AMAX is DOUBLE PRECISION
!>          Absolute value of largest matrix element.  If AMAX is very
!>          close to overflow or very close to underflow, the matrix
!>          should be scaled.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 112 of file zpoequ.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 DOUBLE PRECISION AMAX, SCOND
121* ..
122* .. Array Arguments ..
123 DOUBLE PRECISION S( * )
124 COMPLEX*16 A( LDA, * )
125* ..
126*
127* =====================================================================
128*
129* .. Parameters ..
130 DOUBLE PRECISION ZERO, ONE
131 parameter( zero = 0.0d+0, one = 1.0d+0 )
132* ..
133* .. Local Scalars ..
134 INTEGER I
135 DOUBLE PRECISION SMIN
136* ..
137* .. External Subroutines ..
138 EXTERNAL xerbla
139* ..
140* .. Intrinsic Functions ..
141 INTRINSIC dble, max, min, 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( 'ZPOEQU', -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 ) = dble( a( 1, 1 ) )
169 smin = s( 1 )
170 amax = s( 1 )
171 DO 10 i = 2, n
172 s( i ) = dble( 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 ZPOEQU
203*

◆ zpoequb()

subroutine zpoequb ( integer n,
complex*16, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) s,
double precision scond,
double precision amax,
integer info )

ZPOEQUB

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

Purpose:
!>
!> ZPOEQUB 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 ZPOEQU 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*16 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 DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, S contains the scale factors for A.
!>
[out]SCOND
!>          SCOND is DOUBLE PRECISION
!>          If INFO = 0, S contains the ratio of the smallest S(i) to
!>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
!>          large nor too small, it is not worth scaling by S.
!>
[out]AMAX
!>          AMAX is DOUBLE PRECISION
!>          Absolute value of largest matrix element.  If AMAX is very
!>          close to overflow or very close to underflow, the matrix
!>          should be scaled.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 118 of file zpoequb.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 DOUBLE PRECISION AMAX, SCOND
127* ..
128* .. Array Arguments ..
129 COMPLEX*16 A( LDA, * )
130 DOUBLE PRECISION S( * )
131* ..
132*
133* =====================================================================
134*
135* .. Parameters ..
136 DOUBLE PRECISION ZERO, ONE
137 parameter( zero = 0.0d+0, one = 1.0d+0 )
138* ..
139* .. Local Scalars ..
140 INTEGER I
141 DOUBLE PRECISION SMIN, BASE, TMP
142* ..
143* .. External Functions ..
144 DOUBLE PRECISION DLAMCH
145 EXTERNAL dlamch
146* ..
147* .. External Subroutines ..
148 EXTERNAL xerbla
149* ..
150* .. Intrinsic Functions ..
151 INTRINSIC max, min, sqrt, log, int, real, dimag
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( 'ZPOEQUB', -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 = dlamch( 'B' )
179 tmp = -0.5d+0 / log( base )
180*
181* Find the minimum and maximum diagonal elements.
182*
183 s( 1 ) = dble( a( 1, 1 ) )
184 smin = s( 1 )
185 amax = s( 1 )
186 DO 10 i = 2, n
187 s( i ) = dble( 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 ZPOEQUB
219*

◆ zporfs()

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

ZPORFS

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

Purpose:
!>
!> ZPORFS 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*16 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*16 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 ZPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>          The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]B
!>          B is COMPLEX*16 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*16 array, dimension (LDX,NRHS)
!>          On entry, the solution matrix X, as computed by ZPOTRS.
!>          On exit, the improved solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]FERR
!>          FERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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 zporfs.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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
194 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
195 $ WORK( * ), X( LDX, * )
196* ..
197*
198* ====================================================================
199*
200* .. Parameters ..
201 INTEGER ITMAX
202 parameter( itmax = 5 )
203 DOUBLE PRECISION ZERO
204 parameter( zero = 0.0d+0 )
205 COMPLEX*16 ONE
206 parameter( one = ( 1.0d+0, 0.0d+0 ) )
207 DOUBLE PRECISION TWO
208 parameter( two = 2.0d+0 )
209 DOUBLE PRECISION THREE
210 parameter( three = 3.0d+0 )
211* ..
212* .. Local Scalars ..
213 LOGICAL UPPER
214 INTEGER COUNT, I, J, K, KASE, NZ
215 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
216 COMPLEX*16 ZDUM
217* ..
218* .. Local Arrays ..
219 INTEGER ISAVE( 3 )
220* ..
221* .. External Subroutines ..
222 EXTERNAL xerbla, zaxpy, zcopy, zhemv, zlacn2, zpotrs
223* ..
224* .. Intrinsic Functions ..
225 INTRINSIC abs, dble, dimag, max
226* ..
227* .. External Functions ..
228 LOGICAL LSAME
229 DOUBLE PRECISION DLAMCH
230 EXTERNAL lsame, dlamch
231* ..
232* .. Statement Functions ..
233 DOUBLE PRECISION CABS1
234* ..
235* .. Statement Function definitions ..
236 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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( 'ZPORFS', -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 = dlamch( 'Epsilon' )
278 safmin = dlamch( '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 zcopy( n, b( 1, j ), 1, work, 1 )
295 CALL zhemv( 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( dble( 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( dble( 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 zpotrs( uplo, n, 1, af, ldaf, work, n, info )
357 CALL zaxpy( 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 ZLACN2 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 zlacn2( 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 zpotrs( 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 zpotrs( 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 ZPORFS
432*

◆ zporfsx()

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

ZPORFSX

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

Purpose:
!>
!>    ZPORFSX 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*16 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*16 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 ZPOTRF.
!>
[in]LDAF
!>          LDAF is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in,out]S
!>          S is DOUBLE PRECISION array, dimension (N)
!>     The scale factors for A.  If EQUED = 'Y', A is multiplied on
!>     the left and right by diag(S).  S is an input argument if FACT =
!>     'F'; otherwise, S is an output argument.  If FACT = 'F' and EQUED
!>     = 'Y', each element of S must be positive.  If S is output, each
!>     element of S is a power of the radix. If S is input, each element
!>     of S should be a power of the radix to ensure a reliable solution
!>     and error estimates. Scaling by powers of the radix does not cause
!>     rounding errors unless the result underflows or overflows.
!>     Rounding errors during scaling lead to refining with a matrix that
!>     is not equivalent to the input matrix, producing error estimates
!>     that may not be reliable.
!>
[in]B
!>          B is COMPLEX*16 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*16 array, dimension (LDX,NRHS)
!>     On entry, the solution matrix X, as computed by ZGETRS.
!>     On exit, the improved solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>     The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>     Reciprocal scaled condition number.  This is an estimate of the
!>     reciprocal Skeel condition number of the matrix A after
!>     equilibration (if done).  If this is less than the machine
!>     precision (in particular, if it is zero), the matrix is singular
!>     to working precision.  Note that the error may still be small even
!>     if this number is very small and the matrix appears ill-
!>     conditioned.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>     Componentwise relative backward error.  This is the
!>     componentwise relative backward error of each solution vector X(j)
!>     (i.e., the smallest relative change in any element of A or B that
!>     makes X(j) an exact solution).
!>
[in]N_ERR_BNDS
!>          N_ERR_BNDS is INTEGER
!>     Number of error bounds to return for each right hand side
!>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
!>     ERR_BNDS_COMP below.
!>
[out]ERR_BNDS_NORM
!>          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     normwise relative error, which is defined as follows:
!>
!>     Normwise relative error in the ith solution vector:
!>             max_j (abs(XTRUE(j,i) - X(j,i)))
!>            ------------------------------
!>                  max_j abs(X(j,i))
!>
!>     The array is indexed by the type of error information as described
!>     below. There currently are up to three pieces of information
!>     returned.
!>
!>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_NORM(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * dlamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * dlamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     componentwise relative error, which is defined as follows:
!>
!>     Componentwise relative error in the ith solution vector:
!>                    abs(XTRUE(j,i) - X(j,i))
!>             max_j ----------------------
!>                         abs(X(j,i))
!>
!>     The array is indexed by the right-hand side i (on which the
!>     componentwise relative error depends), and the type of error
!>     information as described below. There currently are up to three
!>     pieces of information returned for each right-hand side. If
!>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
!>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * dlamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * dlamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]NPARAMS
!>          NPARAMS is INTEGER
!>     Specifies the number of parameters set in PARAMS.  If <= 0, the
!>     PARAMS array is never referenced and default values are used.
!>
[in,out]PARAMS
!>          PARAMS is DOUBLE PRECISION array, dimension NPARAMS
!>     Specifies algorithm parameters.  If an entry is < 0.0, then
!>     that entry will be filled with default value used for that
!>     parameter.  Only positions up to NPARAMS are accessed; defaults
!>     are used for higher-numbered parameters.
!>
!>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
!>            refinement or not.
!>         Default: 1.0D+0
!>            = 0.0:  No refinement is performed, and no error bounds are
!>                    computed.
!>            = 1.0:  Use the double-precision refinement algorithm,
!>                    possibly with doubled-single computations if the
!>                    compilation environment does not support DOUBLE
!>                    PRECISION.
!>              (other values are reserved for future use)
!>
!>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
!>            computations allowed for refinement.
!>         Default: 10
!>         Aggressive: Set to 100 to permit convergence using approximate
!>                     factorizations or factorizations other than LU. If
!>                     the factorization uses a technique other than
!>                     Gaussian elimination, the guarantees in
!>                     err_bnds_norm and err_bnds_comp may no longer be
!>                     trustworthy.
!>
!>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
!>            will attempt to find a solution with small componentwise
!>            relative error in the double-precision algorithm.  Positive
!>            is true, 0.0 is false.
!>         Default: 1.0 (attempt componentwise convergence)
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION 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 zporfsx.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 DOUBLE PRECISION RCOND
403* ..
404* .. Array Arguments ..
405 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
406 $ X( LDX, * ), WORK( * )
407 DOUBLE PRECISION RWORK( * ), S( * ), PARAMS(*), BERR( * ),
408 $ ERR_BNDS_NORM( NRHS, * ),
409 $ ERR_BNDS_COMP( NRHS, * )
410* ..
411*
412* ==================================================================
413*
414* .. Parameters ..
415 DOUBLE PRECISION ZERO, ONE
416 parameter( zero = 0.0d+0, one = 1.0d+0 )
417 DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT
418 DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
419 DOUBLE PRECISION DZTHRESH_DEFAULT
420 parameter( itref_default = 1.0d+0 )
421 parameter( ithresh_default = 10.0d+0 )
422 parameter( componentwise_default = 1.0d+0 )
423 parameter( rthresh_default = 0.5d+0 )
424 parameter( dzthresh_default = 0.25d+0 )
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 DOUBLE PRECISION ANORM, RCOND_TMP
441 DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
442 LOGICAL IGNORE_CWISE
443 INTEGER ITHRESH
444 DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH
445* ..
446* .. External Subroutines ..
447 EXTERNAL xerbla, zpocon, zla_porfsx_extended
448* ..
449* .. Intrinsic Functions ..
450 INTRINSIC max, sqrt, transfer
451* ..
452* .. External Functions ..
453 EXTERNAL lsame, ilaprec
454 EXTERNAL dlamch, zlanhe, zla_porcond_x, zla_porcond_c
455 DOUBLE PRECISION DLAMCH, ZLANHE, ZLA_PORCOND_X, ZLA_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.0d+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 = dble( n ) * dlamch( 'Epsilon' )
476 ithresh = int( ithresh_default )
477 rthresh = rthresh_default
478 unstable_thresh = dzthresh_default
479 ignore_cwise = componentwise_default .EQ. 0.0d+0
480*
481 IF ( nparams.GE.la_linrx_ithresh_i ) THEN
482 IF ( params(la_linrx_ithresh_i ).LT.0.0d+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.0d+0 ) THEN
490 IF ( ignore_cwise ) THEN
491 params( la_linrx_cwise_i ) = 0.0d+0
492 ELSE
493 params( la_linrx_cwise_i ) = 1.0d+0
494 END IF
495 ELSE
496 ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+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( 'ZPORFSX', -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.0d+0
537 DO j = 1, nrhs
538 berr( j ) = 0.0d+0
539 IF ( n_err_bnds .GE. 1 ) THEN
540 err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
541 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
542 END IF
543 IF ( n_err_bnds .GE. 2 ) THEN
544 err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
545 err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
546 END IF
547 IF ( n_err_bnds .GE. 3 ) THEN
548 err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
549 err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
550 END IF
551 END DO
552 RETURN
553 END IF
554*
555* Default to failure.
556*
557 rcond = 0.0d+0
558 DO j = 1, nrhs
559 berr( j ) = 1.0d+0
560 IF ( n_err_bnds .GE. 1 ) THEN
561 err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
562 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
563 END IF
564 IF ( n_err_bnds .GE. 2 ) THEN
565 err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
566 err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
567 END IF
568 IF ( n_err_bnds .GE. 3 ) THEN
569 err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
570 err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+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 = zlanhe( norm, uplo, n, a, lda, rwork )
579 CALL zpocon( 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( 'E' )
587
588 CALL zla_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.0d+0, sqrt( dble( n ) ) ) * dlamch( '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 = zla_porcond_c( uplo, n, a, lda, af, ldaf,
604 $ s, .true., info, work, rwork )
605 ELSE
606 rcond_tmp = zla_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.0d+0 )
615 $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+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.0d+0
621 err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+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.0d+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( dlamch( 'Epsilon' ) )
649 DO j = 1, nrhs
650 IF (err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
651 $ THEN
652 rcond_tmp = zla_porcond_x( uplo, n, a, lda, af, ldaf,
653 $ x(1,j), info, work, rwork )
654 ELSE
655 rcond_tmp = 0.0d+0
656 END IF
657*
658* Cap the error at 1.0.
659*
660 IF ( n_err_bnds .GE. la_linrx_err_i
661 $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
662 $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
663*
664* Threshold the error (see LAWN).
665*
666 IF (rcond_tmp .LT. illrcond_thresh) THEN
667 err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
668 err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
669 IF ( params( la_linrx_cwise_i ) .EQ. 1.0d+0
670 $ .AND. info.LT.n + j ) info = n + j
671 ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
672 $ .LT. err_lbnd ) THEN
673 err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
674 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
675 END IF
676*
677* Save the condition number.
678*
679 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
680 err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
681 END IF
682
683 END DO
684 END IF
685*
686 RETURN
687*
688* End of ZPORFSX
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
zlanhe
double precision function zlanhe(norm, uplo, n, a, lda, work)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlanhe.f:124
zla_porfsx_extended
subroutine zla_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)
ZLA_PORFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric or H...
Definition zla_porfsx_extended.f:387
zpocon
subroutine zpocon(uplo, n, a, lda, anorm, rcond, work, rwork, info)
ZPOCON
Definition zpocon.f:121

◆ zpotf2()

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

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

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

Purpose:
!>
!> ZPOTF2 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*16 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 zpotf2.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*16 A( LDA, * )
120* ..
121*
122* =====================================================================
123*
124* .. Parameters ..
125 DOUBLE PRECISION ONE, ZERO
126 parameter( one = 1.0d+0, zero = 0.0d+0 )
127 COMPLEX*16 CONE
128 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
129* ..
130* .. Local Scalars ..
131 LOGICAL UPPER
132 INTEGER J
133 DOUBLE PRECISION AJJ
134* ..
135* .. External Functions ..
136 LOGICAL LSAME, DISNAN
137 COMPLEX*16 ZDOTC
138 EXTERNAL lsame, zdotc, disnan
139* ..
140* .. External Subroutines ..
141 EXTERNAL xerbla, zdscal, zgemv, zlacgv
142* ..
143* .. Intrinsic Functions ..
144 INTRINSIC dble, max, 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( 'ZPOTF2', -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 = dble( a( j, j ) ) - dble( zdotc( j-1, a( 1, j ), 1,
178 $ a( 1, j ), 1 ) )
179 IF( ajj.LE.zero.OR.disnan( ajj ) ) THEN
180 a( j, j ) = ajj
181 GO TO 30
182 END IF
183 ajj = sqrt( ajj )
184 a( j, j ) = ajj
185*
186* Compute elements J+1:N of row J.
187*
188 IF( j.LT.n ) THEN
189 CALL zlacgv( j-1, a( 1, j ), 1 )
190 CALL zgemv( 'Transpose', j-1, n-j, -cone, a( 1, j+1 ),
191 $ lda, a( 1, j ), 1, cone, a( j, j+1 ), lda )
192 CALL zlacgv( j-1, a( 1, j ), 1 )
193 CALL zdscal( 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 = dble( a( j, j ) ) - dble( zdotc( j-1, a( j, 1 ), lda,
205 $ a( j, 1 ), lda ) )
206 IF( ajj.LE.zero.OR.disnan( 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 zlacgv( j-1, a( j, 1 ), lda )
217 CALL zgemv( '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 zlacgv( j-1, a( j, 1 ), lda )
220 CALL zdscal( 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 ZPOTF2
233*
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:59
zlacgv
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:74
zdotc
complex *16 function zdotc(n, zx, incx, zy, incy)
ZDOTC
Definition zdotc.f:83
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:158

◆ zpotrf()

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

ZPOTRF

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

Purpose:
!>
!> ZPOTRF 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*16 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 zpotrf.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*16 A( LDA, * )
118* ..
119*
120* =====================================================================
121*
122* .. Parameters ..
123 DOUBLE PRECISION ONE
124 COMPLEX*16 CONE
125 parameter( one = 1.0d+0, cone = ( 1.0d+0, 0.0d+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 xerbla, zgemm, zherk, zpotrf2, ztrsm
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( 'ZPOTRF', -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, 'ZPOTRF', uplo, n, -1, -1, -1 )
168 IF( nb.LE.1 .OR. nb.GE.n ) THEN
169*
170* Use unblocked code.
171*
172 CALL zpotrf2( 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 zherk( 'Upper', 'Conjugate transpose', jb, j-1,
188 $ -one, a( 1, j ), lda, one, a( j, j ), lda )
189 CALL zpotrf2( '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 zgemm( '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 ztrsm( '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 zherk( 'Lower', 'No transpose', jb, j-1, -one,
217 $ a( j, 1 ), lda, one, a( j, j ), lda )
218 CALL zpotrf2( '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 zgemm( '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 ztrsm( '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 ZPOTRF
245*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
zpotrf2
recursive subroutine zpotrf2(uplo, n, a, lda, info)
ZPOTRF2
Definition zpotrf2.f:106
zherk
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
Definition zherk.f:173
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:187
ztrsm
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180

◆ zpotrf2()

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

ZPOTRF2

Purpose:
!>
!> ZPOTRF2 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 call itself to factor A22.
!>
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in,out]A
!>          A is COMPLEX*16 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 zpotrf2.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*16 A( LDA, * )
117* ..
118*
119* =====================================================================
120*
121* .. Parameters ..
122 DOUBLE PRECISION ONE, ZERO
123 parameter( one = 1.0d+0, zero = 0.0d+0 )
124 COMPLEX*16 CONE
125 parameter( cone = (1.0d+0, 0.0d+0) )
126* ..
127* .. Local Scalars ..
128 LOGICAL UPPER
129 INTEGER N1, N2, IINFO
130 DOUBLE PRECISION AJJ
131* ..
132* .. External Functions ..
133 LOGICAL LSAME, DISNAN
134 EXTERNAL lsame, disnan
135* ..
136* .. External Subroutines ..
137 EXTERNAL zherk, ztrsm, xerbla
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC max, dble, 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( 'ZPOTRF2', -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 = dble( a( 1, 1 ) )
172 IF( ajj.LE.zero.OR.disnan( 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 zpotrf2( 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 ztrsm( '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 zherk( uplo, 'C', n2, n1, -one, a( 1, n1+1 ), lda,
207 $ one, a( n1+1, n1+1 ), lda )
208 CALL zpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
209 IF ( iinfo.NE.0 ) THEN
210 info = iinfo + n1
211 RETURN
212 END IF
213*
214* Compute the Cholesky factorization A = L*L**H
215*
216 ELSE
217*
218* Update and scale A21
219*
220 CALL ztrsm( 'R', 'L', 'C', 'N', n2, n1, cone,
221 $ a( 1, 1 ), lda, a( n1+1, 1 ), lda )
222*
223* Update and factor A22
224*
225 CALL zherk( uplo, 'N', n2, n1, -one, a( n1+1, 1 ), lda,
226 $ one, a( n1+1, n1+1 ), lda )
227 CALL zpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
228 IF ( iinfo.NE.0 ) THEN
229 info = iinfo + n1
230 RETURN
231 END IF
232 END IF
233 END IF
234 RETURN
235*
236* End of ZPOTRF2
237*

◆ zpotri()

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

ZPOTRI

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

Purpose:
!>
!> ZPOTRI 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 ZPOTRF.
!>
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*16 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
!>          ZPOTRF.
!>          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 zpotri.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*16 A( LDA, * )
106* ..
107*
108* =====================================================================
109*
110* .. External Functions ..
111 LOGICAL LSAME
112 EXTERNAL lsame
113* ..
114* .. External Subroutines ..
115 EXTERNAL xerbla, zlauum, ztrtri
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( 'ZPOTRI', -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 ztrtri( 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 zlauum( uplo, n, a, lda, info )
151*
152 RETURN
153*
154* End of ZPOTRI
155*
zlauum
subroutine zlauum(uplo, n, a, lda, info)
ZLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked...
Definition zlauum.f:102
ztrtri
subroutine ztrtri(uplo, diag, n, a, lda, info)
ZTRTRI
Definition ztrtri.f:109

◆ zpotrs()

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

ZPOTRS

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

Purpose:
!>
!> ZPOTRS 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 ZPOTRF.
!>
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*16 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 ZPOTRF.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is COMPLEX*16 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 zpotrs.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*16 A( LDA, * ), B( LDB, * )
121* ..
122*
123* =====================================================================
124*
125* .. Parameters ..
126 COMPLEX*16 ONE
127 parameter( one = ( 1.0d+0, 0.0d+0 ) )
128* ..
129* .. Local Scalars ..
130 LOGICAL UPPER
131* ..
132* .. External Functions ..
133 LOGICAL LSAME
134 EXTERNAL lsame
135* ..
136* .. External Subroutines ..
137 EXTERNAL xerbla, ztrsm
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( 'ZPOTRS', -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 ztrsm( '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 ztrsm( '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 ztrsm( '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 ztrsm( 'Left', 'Lower', 'Conjugate transpose', 'Non-unit',
194 $ n, nrhs, one, a, lda, b, ldb )
195 END IF
196*
197 RETURN
198*
199* End of ZPOTRS
200*