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

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. More...
 
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. More...
 
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. More...
 
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. More...
 
subroutine zpocon (UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
 ZPOCON More...
 
subroutine zpoequ (N, A, LDA, S, SCOND, AMAX, INFO)
 ZPOEQU More...
 
subroutine zpoequb (N, A, LDA, S, SCOND, AMAX, INFO)
 ZPOEQUB More...
 
subroutine zporfs (UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
 ZPORFS More...
 
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 More...
 
subroutine zpotf2 (UPLO, N, A, LDA, INFO)
 ZPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm). More...
 
subroutine zpotrf (UPLO, N, A, LDA, INFO)
 ZPOTRF More...
 
recursive subroutine zpotrf2 (UPLO, N, A, LDA, INFO)
 ZPOTRF2 More...
 
subroutine zpotri (UPLO, N, A, LDA, INFO)
 ZPOTRI More...
 
subroutine zpotrs (UPLO, N, NRHS, A, LDA, B, LDB, INFO)
 ZPOTRS More...
 

Detailed Description

This is the group of complex16 computational functions for PO matrices

Function Documentation

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.
[in]WORK
          WORK is COMPLEX*16 array, dimension (2*N).
     Workspace.
[in]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 134 of file zla_porcond_c.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
[in]WORK
          WORK is COMPLEX*16 array, dimension (2*N).
     Workspace.
[in]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 127 of file zla_porcond_x.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 resonsible for setting the second fields of
 ERR_BNDS_NORM and ERR_BNDS_COMP.
Parameters
[in]PREC_TYPE
          PREC_TYPE is INTEGER
     Specifies the intermediate precision to be used in refinement.
     The value is defined by ILAPREC(P) where P is a CHARACTER and
     P    = 'S':  Single
          = 'D':  Double
          = 'I':  Indigenous
          = 'X', 'E':  Extra
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
     The number of linear equations, i.e., the order of the
     matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
     The number of right-hand-sides, i.e., the number of columns of the
     matrix B.
[in]A
          A is 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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * slamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated normwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * slamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*A, where S scales each row by a power of the
              radix so all absolute row sums of Z are approximately 1.

     This subroutine is only responsible for setting the second field
     above.
     See Lapack Working Note 165 for further details and extra
     cautions.
[in,out]ERR_BNDS_COMP
          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension
                    (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     componentwise relative error, which is defined as follows:

     Componentwise relative error in the ith solution vector:
                    abs(XTRUE(j,i) - X(j,i))
             max_j ----------------------
                         abs(X(j,i))

     The array is indexed by the right-hand side i (on which the
     componentwise relative error depends), and the type of error
     information as described below. There currently are up to three
     pieces of information returned for each right-hand side. If
     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
     the first (:,N_ERR_BNDS) entries are returned.

     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_COMP(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * slamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated componentwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * slamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*(A*diag(x)), where x is the solution for the
              current right-hand side and S scales each row of
              A*diag(x) by a power of the radix so all absolute row
              sums of Z are approximately 1.

     This subroutine is only responsible for setting the second field
     above.
     See Lapack Working Note 165 for further details and extra
     cautions.
[in]RES
          RES is 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 definte as the relative
     change in each component being less than DZ_UB. The default value
     is 0.25, requiring the first bit to be stable. See LAWN 165 for
     more details.
[in]IGNORE_CWISE
          IGNORE_CWISE is LOGICAL
     If .TRUE. then ignore componentwise convergence. Default value
     is .FALSE..
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit.
       < 0:  if INFO = -i, the ith argument to ZPOTRS had an illegal
             value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 392 of file zla_porfsx_extended.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 "max absolute element" norm is used. If this is
 much less than 1, the stability of the LU factorization of the
 (equilibrated) matrix A could be poor. This also means that the
 solution X, estimated condition numbers, and error bounds could be
 unreliable.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]NCOLS
          NCOLS is INTEGER
     The number of columns of the matrix A. NCOLS >= 0.
[in]A
          A is 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]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 109 of file zla_porpvgrw.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 123 of file zpocon.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 115 of file zpoequ.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine 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
 symmetric positive definite matrix A and reduce its condition number
 (with respect to the two-norm).  S contains the scale factors,
 S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
 elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
 choice of S puts the condition number of B within a factor N of the
 smallest possible condition number over all possible diagonal
 scalings.
Parameters
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The N-by-N symmetric positive definite matrix whose scaling
          factors are to be computed.  Only the diagonal elements of A
          are referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]S
          S is DOUBLE PRECISION array, dimension (N)
          If INFO = 0, S contains the scale factors for A.
[out]SCOND
          SCOND is DOUBLE PRECISION
          If INFO = 0, S contains the ratio of the smallest S(i) to
          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
          large nor too small, it is not worth scaling by S.
[out]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 115 of file zpoequb.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 185 of file zporfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 symmetric positive
    definite, and provides error bounds and backward error estimates
    for the solution.  In addition to normwise error bound, the code
    provides maximum componentwise error bound if possible.  See
    comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
    error bounds.

    The original system of linear equations may have been equilibrated
    before calling this routine, as described by arguments EQUED and S
    below. In this case, the solution and error bounds returned are
    for the original unequilibrated system.
     Some optional parameters are bundled in the PARAMS array.  These
     settings determine how refinement is performed, but often the
     defaults are acceptable.  If the defaults are acceptable, users
     can pass NPARAMS = 0 which prevents the source code from accessing
     the PARAMS argument.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]EQUED
          EQUED is CHARACTER*1
     Specifies the form of equilibration that was done to A
     before calling this routine. This is needed to compute
     the solution and error bounds correctly.
       = 'N':  No equilibration
       = 'Y':  Both row and column equilibration, i.e., A has been
               replaced by diag(S) * A * diag(S).
               The right hand side B has been changed accordingly.
[in]N
          N is INTEGER
     The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
     The number of right hand sides, i.e., the number of columns
     of the matrices B and X.  NRHS >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
     The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
     upper triangular part of A contains the upper triangular part
     of the matrix A, and the strictly lower triangular part of A
     is not referenced.  If UPLO = 'L', the leading N-by-N lower
     triangular part of A contains the lower triangular part of
     the matrix A, and the strictly upper triangular part of A is
     not referenced.
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is 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 DPOTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= max(1,N).
[in,out]S
          S is DOUBLE PRECISION array, dimension (N)
     The row scale factors for A.  If EQUED = 'Y', A is multiplied on
     the left and right by diag(S).  S is an input argument if FACT =
     'F'; otherwise, S is an output argument.  If FACT = 'F' and EQUED
     = 'Y', each element of S must be positive.  If S is output, each
     element of S is a power of the radix. If S is input, each element
     of S should be a power of the radix to ensure a reliable solution
     and error estimates. Scaling by powers of the radix does not cause
     rounding errors unless the result underflows or overflows.
     Rounding errors during scaling lead to refining with a matrix that
     is not equivalent to the input matrix, producing error estimates
     that may not be reliable.
[in]B
          B is 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 DGETRS.
     On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
     The leading dimension of the array X.  LDX >= max(1,N).
[out]RCOND
          RCOND is DOUBLE PRECISION
     Reciprocal scaled condition number.  This is an estimate of the
     reciprocal Skeel condition number of the matrix A after
     equilibration (if done).  If this is less than the machine
     precision (in particular, if it is zero), the matrix is singular
     to working precision.  Note that the error may still be small even
     if this number is very small and the matrix appears ill-
     conditioned.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
     Componentwise relative backward error.  This is the
     componentwise relative backward error of each solution vector X(j)
     (i.e., the smallest relative change in any element of A or B that
     makes X(j) an exact solution).
[in]N_ERR_BNDS
          N_ERR_BNDS is INTEGER
     Number of error bounds to return for each right hand side
     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
     ERR_BNDS_COMP below.
[out]ERR_BNDS_NORM
          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     normwise relative error, which is defined as follows:

     Normwise relative error in the ith solution vector:
             max_j (abs(XTRUE(j,i) - X(j,i)))
            ------------------------------
                  max_j abs(X(j,i))

     The array is indexed by the type of error information as described
     below. There currently are up to three pieces of information
     returned.

     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_NORM(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated normwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*A, where S scales each row by a power of the
              radix so all absolute row sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[out]ERR_BNDS_COMP
          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     componentwise relative error, which is defined as follows:

     Componentwise relative error in the ith solution vector:
                    abs(XTRUE(j,i) - X(j,i))
             max_j ----------------------
                         abs(X(j,i))

     The array is indexed by the right-hand side i (on which the
     componentwise relative error depends), and the type of error
     information as described below. There currently are up to three
     pieces of information returned for each right-hand side. If
     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
     the first (:,N_ERR_BNDS) entries are returned.

     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_COMP(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated componentwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*(A*diag(x)), where x is the solution for the
              current right-hand side and S scales each row of
              A*diag(x) by a power of the radix so all absolute row
              sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[in]NPARAMS
          NPARAMS is INTEGER
     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
     PARAMS array is never referenced and default values are used.
[in,out]PARAMS
          PARAMS is DOUBLE PRECISION array, dimension NPARAMS
     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
     that entry will be filled with default value used for that
     parameter.  Only positions up to NPARAMS are accessed; defaults
     are used for higher-numbered parameters.

       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
            refinement or not.
         Default: 1.0D+0
            = 0.0 : No refinement is performed, and no error bounds are
                    computed.
            = 1.0 : Use the double-precision refinement algorithm,
                    possibly with doubled-single computations if the
                    compilation environment does not support DOUBLE
                    PRECISION.
              (other values are reserved for future use)

       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
            computations allowed for refinement.
         Default: 10
         Aggressive: Set to 100 to permit convergence using approximate
                     factorizations or factorizations other than LU. If
                     the factorization uses a technique other than
                     Gaussian elimination, the guarantees in
                     err_bnds_norm and err_bnds_comp may no longer be
                     trustworthy.

       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
            will attempt to find a solution with small componentwise
            relative error in the double-precision algorithm.  Positive
            is true, 0.0 is false.
         Default: 1.0 (attempt componentwise convergence)
[out]WORK
          WORK is 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.
Date
April 2012

Definition at line 395 of file zporfsx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 111 of file zpotf2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 109 of file zpotrf.f.

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

Here is the call graph for this function:

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 real symmetric
 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 symmetric matrix A.  If UPLO = 'U', the leading
          N-by-N upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading N-by-N lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.

          On exit, if INFO = 0, the factor U or L from the Cholesky
          factorization A = U**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.
Date
November 2015

Definition at line 108 of file zpotrf2.f.

108 *
109 * -- LAPACK computational routine (version 3.6.0) --
110 * -- LAPACK is a software package provided by Univ. of Tennessee, --
111 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
112 * November 2015
113 *
114 * .. Scalar Arguments ..
115  CHARACTER uplo
116  INTEGER info, lda, n
117 * ..
118 * .. Array Arguments ..
119  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 n1, n2, iinfo
133  DOUBLE PRECISION ajj
134 * ..
135 * .. External Functions ..
136  LOGICAL lsame, disnan
137  EXTERNAL lsame, disnan
138 * ..
139 * .. External Subroutines ..
140  EXTERNAL zherk, ztrsm, xerbla
141 * ..
142 * .. Intrinsic Functions ..
143  INTRINSIC max, dble, sqrt
144 * ..
145 * .. Executable Statements ..
146 *
147 * Test the input parameters
148 *
149  info = 0
150  upper = lsame( uplo, 'U' )
151  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
152  info = -1
153  ELSE IF( n.LT.0 ) THEN
154  info = -2
155  ELSE IF( lda.LT.max( 1, n ) ) THEN
156  info = -4
157  END IF
158  IF( info.NE.0 ) THEN
159  CALL xerbla( 'ZPOTRF2', -info )
160  RETURN
161  END IF
162 *
163 * Quick return if possible
164 *
165  IF( n.EQ.0 )
166  $ RETURN
167 *
168 * N=1 case
169 *
170  IF( n.EQ.1 ) THEN
171 *
172 * Test for non-positive-definiteness
173 *
174  ajj = dble( a( 1, 1 ) )
175  IF( ajj.LE.zero.OR.disnan( ajj ) ) THEN
176  info = 1
177  RETURN
178  END IF
179 *
180 * Factor
181 *
182  a( 1, 1 ) = sqrt( ajj )
183 *
184 * Use recursive code
185 *
186  ELSE
187  n1 = n/2
188  n2 = n-n1
189 *
190 * Factor A11
191 *
192  CALL zpotrf2( uplo, n1, a( 1, 1 ), lda, iinfo )
193  IF ( iinfo.NE.0 ) THEN
194  info = iinfo
195  RETURN
196  END IF
197 *
198 * Compute the Cholesky factorization A = U**H*U
199 *
200  IF( upper ) THEN
201 *
202 * Update and scale A12
203 *
204  CALL ztrsm( 'L', 'U', 'C', 'N', n1, n2, cone,
205  $ a( 1, 1 ), lda, a( 1, n1+1 ), lda )
206 *
207 * Update and factor A22
208 *
209  CALL zherk( uplo, 'C', n2, n1, -one, a( 1, n1+1 ), lda,
210  $ one, a( n1+1, n1+1 ), lda )
211  CALL zpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
212  IF ( iinfo.NE.0 ) THEN
213  info = iinfo + n1
214  RETURN
215  END IF
216 *
217 * Compute the Cholesky factorization A = L*L**H
218 *
219  ELSE
220 *
221 * Update and scale A21
222 *
223  CALL ztrsm( 'R', 'L', 'C', 'N', n2, n1, cone,
224  $ a( 1, 1 ), lda, a( n1+1, 1 ), lda )
225 *
226 * Update and factor A22
227 *
228  CALL zherk( uplo, 'N', n2, n1, -one, a( n1+1, 1 ), lda,
229  $ one, a( n1+1, n1+1 ), lda )
230  CALL zpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo )
231  IF ( iinfo.NE.0 ) THEN
232  info = iinfo + n1
233  RETURN
234  END IF
235  END IF
236  END IF
237  RETURN
238 *
239 * End of ZPOTRF2
240 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
recursive subroutine zpotrf2(UPLO, N, A, LDA, INFO)
ZPOTRF2
Definition: zpotrf2.f:108
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:175

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 97 of file zpotri.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 112 of file zpotrs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: