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

Functions

subroutine zcposv (UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK, SWORK, RWORK, ITER, INFO)
  ZCPOSV computes the solution to system of linear equations A * X = B for PO matrices More...
 
subroutine zposv (UPLO, N, NRHS, A, LDA, B, LDB, INFO)
  ZPOSV computes the solution to system of linear equations A * X = B for PO matrices More...
 
subroutine zposvx (FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO)
  ZPOSVX computes the solution to system of linear equations A * X = B for PO matrices More...
 
subroutine zposvxx (FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, S, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO)
  ZPOSVXX computes the solution to system of linear equations A * X = B for PO matrices More...
 

Detailed Description

This is the group of complex16 solve driver functions for PO matrices

Function Documentation

subroutine zcposv ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
complex*16, dimension( n, * )  WORK,
complex, dimension( * )  SWORK,
double precision, dimension( * )  RWORK,
integer  ITER,
integer  INFO 
)

ZCPOSV computes the solution to system of linear equations A * X = B for PO matrices

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

Purpose:
 ZCPOSV computes the solution to a complex system of linear equations
    A * X = B,
 where A is an N-by-N Hermitian positive definite matrix and X and B
 are N-by-NRHS matrices.

 ZCPOSV first attempts to factorize the matrix in COMPLEX and use this
 factorization within an iterative refinement procedure to produce a
 solution with COMPLEX*16 normwise backward error quality (see below).
 If the approach fails the method switches to a COMPLEX*16
 factorization and solve.

 The iterative refinement is not going to be a winning strategy if
 the ratio COMPLEX performance over COMPLEX*16 performance is too
 small. A reasonable strategy should take the number of right-hand
 sides and the size of the matrix into account. This might be done
 with a call to ILAENV in the future. Up to now, we always try
 iterative refinement.

 The iterative refinement process is stopped if
     ITER > ITERMAX
 or for all the RHS we have:
     RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
 where
     o ITER is the number of the current iteration in the iterative
       refinement process
     o RNRM is the infinity-norm of the residual
     o XNRM is the infinity-norm of the solution
     o ANRM is the infinity-operator-norm of the matrix A
     o EPS is the machine epsilon returned by DLAMCH('Epsilon')
 The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
 respectively.
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]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 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.

          Note that the imaginary parts of the diagonal
          elements need not be set and are assumed to be zero.

          On exit, if iterative refinement has been successfully used
          (INFO.EQ.0 and ITER.GE.0, see description below), then A is
          unchanged, if double precision factorization has been used
          (INFO.EQ.0 and ITER.LT.0, see description below), then the
          array A contains 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).
[in]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          The N-by-NRHS right hand side matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
          If INFO = 0, the N-by-NRHS solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]WORK
          WORK is COMPLEX*16 array, dimension (N*NRHS)
          This array is used to hold the residual vectors.
[out]SWORK
          SWORK is COMPLEX array, dimension (N*(N+NRHS))
          This array is used to use the single precision matrix and the
          right-hand sides or solutions in single precision.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[out]ITER
          ITER is INTEGER
          < 0: iterative refinement has failed, COMPLEX*16
               factorization has been performed
               -1 : the routine fell back to full precision for
                    implementation- or machine-specific reasons
               -2 : narrowing the precision induced an overflow,
                    the routine fell back to full precision
               -3 : failure of CPOTRF
               -31: stop the iterative refinement after the 30th
                    iterations
          > 0: iterative refinement has been sucessfully used.
               Returns the number of iterations
[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 of
                (COMPLEX*16) A is not positive definite, so the
                factorization could not be completed, and the solution
                has not been computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 211 of file zcposv.f.

211 *
212 * -- LAPACK driver routine (version 3.4.0) --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215 * November 2011
216 *
217 * .. Scalar Arguments ..
218  CHARACTER uplo
219  INTEGER info, iter, lda, ldb, ldx, n, nrhs
220 * ..
221 * .. Array Arguments ..
222  DOUBLE PRECISION rwork( * )
223  COMPLEX swork( * )
224  COMPLEX*16 a( lda, * ), b( ldb, * ), work( n, * ),
225  $ x( ldx, * )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  LOGICAL doitref
232  parameter( doitref = .true. )
233 *
234  INTEGER itermax
235  parameter( itermax = 30 )
236 *
237  DOUBLE PRECISION bwdmax
238  parameter( bwdmax = 1.0e+00 )
239 *
240  COMPLEX*16 negone, one
241  parameter( negone = ( -1.0d+00, 0.0d+00 ),
242  $ one = ( 1.0d+00, 0.0d+00 ) )
243 *
244 * .. Local Scalars ..
245  INTEGER i, iiter, ptsa, ptsx
246  DOUBLE PRECISION anrm, cte, eps, rnrm, xnrm
247  COMPLEX*16 zdum
248 *
249 * .. External Subroutines ..
250  EXTERNAL zaxpy, zhemm, zlacpy, zlat2c, zlag2c, clag2z,
251  $ cpotrf, cpotrs, xerbla
252 * ..
253 * .. External Functions ..
254  INTEGER izamax
255  DOUBLE PRECISION dlamch, zlanhe
256  LOGICAL lsame
257  EXTERNAL izamax, dlamch, zlanhe, lsame
258 * ..
259 * .. Intrinsic Functions ..
260  INTRINSIC abs, dble, max, sqrt
261 * .. Statement Functions ..
262  DOUBLE PRECISION cabs1
263 * ..
264 * .. Statement Function definitions ..
265  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
266 * ..
267 * .. Executable Statements ..
268 *
269  info = 0
270  iter = 0
271 *
272 * Test the input parameters.
273 *
274  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
275  info = -1
276  ELSE IF( n.LT.0 ) THEN
277  info = -2
278  ELSE IF( nrhs.LT.0 ) THEN
279  info = -3
280  ELSE IF( lda.LT.max( 1, n ) ) THEN
281  info = -5
282  ELSE IF( ldb.LT.max( 1, n ) ) THEN
283  info = -7
284  ELSE IF( ldx.LT.max( 1, n ) ) THEN
285  info = -9
286  END IF
287  IF( info.NE.0 ) THEN
288  CALL xerbla( 'ZCPOSV', -info )
289  RETURN
290  END IF
291 *
292 * Quick return if (N.EQ.0).
293 *
294  IF( n.EQ.0 )
295  $ RETURN
296 *
297 * Skip single precision iterative refinement if a priori slower
298 * than double precision factorization.
299 *
300  IF( .NOT.doitref ) THEN
301  iter = -1
302  GO TO 40
303  END IF
304 *
305 * Compute some constants.
306 *
307  anrm = zlanhe( 'I', uplo, n, a, lda, rwork )
308  eps = dlamch( 'Epsilon' )
309  cte = anrm*eps*sqrt( dble( n ) )*bwdmax
310 *
311 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
312 *
313  ptsa = 1
314  ptsx = ptsa + n*n
315 *
316 * Convert B from double precision to single precision and store the
317 * result in SX.
318 *
319  CALL zlag2c( n, nrhs, b, ldb, swork( ptsx ), n, info )
320 *
321  IF( info.NE.0 ) THEN
322  iter = -2
323  GO TO 40
324  END IF
325 *
326 * Convert A from double precision to single precision and store the
327 * result in SA.
328 *
329  CALL zlat2c( uplo, n, a, lda, swork( ptsa ), n, info )
330 *
331  IF( info.NE.0 ) THEN
332  iter = -2
333  GO TO 40
334  END IF
335 *
336 * Compute the Cholesky factorization of SA.
337 *
338  CALL cpotrf( uplo, n, swork( ptsa ), n, info )
339 *
340  IF( info.NE.0 ) THEN
341  iter = -3
342  GO TO 40
343  END IF
344 *
345 * Solve the system SA*SX = SB.
346 *
347  CALL cpotrs( uplo, n, nrhs, swork( ptsa ), n, swork( ptsx ), n,
348  $ info )
349 *
350 * Convert SX back to COMPLEX*16
351 *
352  CALL clag2z( n, nrhs, swork( ptsx ), n, x, ldx, info )
353 *
354 * Compute R = B - AX (R is WORK).
355 *
356  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
357 *
358  CALL zhemm( 'Left', uplo, n, nrhs, negone, a, lda, x, ldx, one,
359  $ work, n )
360 *
361 * Check whether the NRHS normwise backward errors satisfy the
362 * stopping criterion. If yes, set ITER=0 and return.
363 *
364  DO i = 1, nrhs
365  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
366  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
367  IF( rnrm.GT.xnrm*cte )
368  $ GO TO 10
369  END DO
370 *
371 * If we are here, the NRHS normwise backward errors satisfy the
372 * stopping criterion. We are good to exit.
373 *
374  iter = 0
375  RETURN
376 *
377  10 CONTINUE
378 *
379  DO 30 iiter = 1, itermax
380 *
381 * Convert R (in WORK) from double precision to single precision
382 * and store the result in SX.
383 *
384  CALL zlag2c( n, nrhs, work, n, swork( ptsx ), n, info )
385 *
386  IF( info.NE.0 ) THEN
387  iter = -2
388  GO TO 40
389  END IF
390 *
391 * Solve the system SA*SX = SR.
392 *
393  CALL cpotrs( uplo, n, nrhs, swork( ptsa ), n, swork( ptsx ), n,
394  $ info )
395 *
396 * Convert SX back to double precision and update the current
397 * iterate.
398 *
399  CALL clag2z( n, nrhs, swork( ptsx ), n, work, n, info )
400 *
401  DO i = 1, nrhs
402  CALL zaxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
403  END DO
404 *
405 * Compute R = B - AX (R is WORK).
406 *
407  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
408 *
409  CALL zhemm( 'L', uplo, n, nrhs, negone, a, lda, x, ldx, one,
410  $ work, n )
411 *
412 * Check whether the NRHS normwise backward errors satisfy the
413 * stopping criterion. If yes, set ITER=IITER>0 and return.
414 *
415  DO i = 1, nrhs
416  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
417  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
418  IF( rnrm.GT.xnrm*cte )
419  $ GO TO 20
420  END DO
421 *
422 * If we are here, the NRHS normwise backward errors satisfy the
423 * stopping criterion, we are good to exit.
424 *
425  iter = iiter
426 *
427  RETURN
428 *
429  20 CONTINUE
430 *
431  30 CONTINUE
432 *
433 * If we are at this place of the code, this is because we have
434 * performed ITER=ITERMAX iterations and never satisified the
435 * stopping criterion, set up the ITER flag accordingly and follow
436 * up on double precision routine.
437 *
438  iter = -itermax - 1
439 *
440  40 CONTINUE
441 *
442 * Single-precision iterative refinement failed to converge to a
443 * satisfactory solution, so we resort to double precision.
444 *
445  CALL zpotrf( uplo, n, a, lda, info )
446 *
447  IF( info.NE.0 )
448  $ RETURN
449 *
450  CALL zlacpy( 'All', n, nrhs, b, ldb, x, ldx )
451  CALL zpotrs( uplo, n, nrhs, a, lda, x, ldx, info )
452 *
453  RETURN
454 *
455 * End of ZCPOSV.
456 *
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
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 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
subroutine zhemm(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZHEMM
Definition: zhemm.f:193
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 cpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
CPOTRS
Definition: cpotrs.f:112
subroutine cpotrf(UPLO, N, A, LDA, INFO)
CPOTRF
Definition: cpotrf.f:109
subroutine zlat2c(UPLO, N, A, LDA, SA, LDSA, INFO)
ZLAT2C converts a double complex triangular matrix to a complex triangular matrix.
Definition: zlat2c.f:113
subroutine clag2z(M, N, SA, LDSA, A, LDA, INFO)
CLAG2Z converts a complex single precision matrix to a complex double precision matrix.
Definition: clag2z.f:105
subroutine zlag2c(M, N, A, LDA, SA, LDSA, INFO)
ZLAG2C converts a complex double precision matrix to a complex single precision matrix.
Definition: zlag2c.f:109

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zposv ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

ZPOSV computes the solution to system of linear equations A * X = B for PO matrices

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

Purpose:
 ZPOSV computes the solution to a complex system of linear equations
    A * X = B,
 where A is an N-by-N Hermitian positive definite matrix and X and B
 are N-by-NRHS matrices.

 The Cholesky decomposition is used to factor A as
    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 a lower triangular
 matrix.  The factored form of A is then used to solve the system of
 equations A * X = B.
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]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 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).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          On entry, the N-by-NRHS right hand side matrix B.
          On exit, if INFO = 0, the N-by-NRHS 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
          > 0:  if INFO = i, the leading minor of order i of A is not
                positive definite, so the factorization could not be
                completed, and the solution has not been computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 132 of file zposv.f.

132 *
133 * -- LAPACK driver routine (version 3.4.0) --
134 * -- LAPACK is a software package provided by Univ. of Tennessee, --
135 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136 * November 2011
137 *
138 * .. Scalar Arguments ..
139  CHARACTER uplo
140  INTEGER info, lda, ldb, n, nrhs
141 * ..
142 * .. Array Arguments ..
143  COMPLEX*16 a( lda, * ), b( ldb, * )
144 * ..
145 *
146 * =====================================================================
147 *
148 * .. External Functions ..
149  LOGICAL lsame
150  EXTERNAL lsame
151 * ..
152 * .. External Subroutines ..
153  EXTERNAL xerbla, zpotrf, zpotrs
154 * ..
155 * .. Intrinsic Functions ..
156  INTRINSIC max
157 * ..
158 * .. Executable Statements ..
159 *
160 * Test the input parameters.
161 *
162  info = 0
163  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
164  info = -1
165  ELSE IF( n.LT.0 ) THEN
166  info = -2
167  ELSE IF( nrhs.LT.0 ) THEN
168  info = -3
169  ELSE IF( lda.LT.max( 1, n ) ) THEN
170  info = -5
171  ELSE IF( ldb.LT.max( 1, n ) ) THEN
172  info = -7
173  END IF
174  IF( info.NE.0 ) THEN
175  CALL xerbla( 'ZPOSV ', -info )
176  RETURN
177  END IF
178 *
179 * Compute the Cholesky factorization A = U**H *U or A = L*L**H.
180 *
181  CALL zpotrf( uplo, n, a, lda, info )
182  IF( info.EQ.0 ) THEN
183 *
184 * Solve the system A*X = B, overwriting B with X.
185 *
186  CALL zpotrs( uplo, n, nrhs, a, lda, b, ldb, info )
187 *
188  END IF
189  RETURN
190 *
191 * End of ZPOSV
192 *
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102
subroutine zpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
ZPOTRS
Definition: zpotrs.f:112
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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 zposvx ( character  FACT,
character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
character  EQUED,
double precision, dimension( * )  S,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
double precision  RCOND,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZPOSVX computes the solution to system of linear equations A * X = B for PO matrices

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

Purpose:
 ZPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
 compute the solution to a complex system of linear equations
    A * X = B,
 where A is an N-by-N Hermitian positive definite matrix and X and B
 are N-by-NRHS matrices.

 Error bounds on the solution and a condition estimate are also
 provided.
Description:
 The following steps are performed:

 1. If FACT = 'E', real scaling factors are computed to equilibrate
    the system:
       diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
    Whether or not the system will be equilibrated depends on the
    scaling of the matrix A, but if equilibration is used, A is
    overwritten by diag(S)*A*diag(S) and B by diag(S)*B.

 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
    factor the matrix A (after equilibration if FACT = 'E') as
       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 a lower triangular
    matrix.

 3. If the leading i-by-i principal minor is not positive definite,
    then the routine returns with INFO = i. Otherwise, the factored
    form of A is used to estimate the condition number of the matrix
    A.  If the reciprocal of the condition number is less than machine
    precision, INFO = N+1 is returned as a warning, but the routine
    still goes on to solve for X and compute error bounds as
    described below.

 4. The system of equations is solved for X using the factored form
    of A.

 5. Iterative refinement is applied to improve the computed solution
    matrix and calculate error bounds and backward error estimates
    for it.

 6. If equilibration was used, the matrix X is premultiplied by
    diag(S) so that it solves the original system before
    equilibration.
Parameters
[in]FACT
          FACT is CHARACTER*1
          Specifies whether or not the factored form of the matrix A is
          supplied on entry, and if not, whether the matrix A should be
          equilibrated before it is factored.
          = 'F':  On entry, AF contains the factored form of A.
                  If EQUED = 'Y', the matrix A has been equilibrated
                  with scaling factors given by S.  A and AF will not
                  be modified.
          = 'N':  The matrix A will be copied to AF and factored.
          = 'E':  The matrix A will be equilibrated if necessary, then
                  copied to AF and factored.
[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 matrices B and X.  NRHS >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the Hermitian matrix A, except if FACT = 'F' and
          EQUED = 'Y', then A must contain the equilibrated matrix
          diag(S)*A*diag(S).  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.  A is not modified if
          FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.

          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
          diag(S)*A*diag(S).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in,out]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
          If FACT = 'F', then AF is an input argument and on entry
          contains the triangular factor U or L from the Cholesky
          factorization A = U**H *U or A = L*L**H, in the same storage
          format as A.  If EQUED .ne. 'N', then AF is the factored form
          of the equilibrated matrix diag(S)*A*diag(S).

          If FACT = 'N', then AF is an output argument and on exit
          returns the triangular factor U or L from the Cholesky
          factorization A = U**H *U or A = L*L**H of the original
          matrix A.

          If FACT = 'E', then AF is an output argument and on exit
          returns the triangular factor U or L from the Cholesky
          factorization A = U**H *U or A = L*L**H of the equilibrated
          matrix A (see the description of A for the form of the
          equilibrated matrix).
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= max(1,N).
[in,out]EQUED
          EQUED is CHARACTER*1
          Specifies the form of equilibration that was done.
          = 'N':  No equilibration (always true if FACT = 'N').
          = 'Y':  Equilibration was done, i.e., A has been replaced by
                  diag(S) * A * diag(S).
          EQUED is an input argument if FACT = 'F'; otherwise, it is an
          output argument.
[in,out]S
          S is DOUBLE PRECISION array, dimension (N)
          The scale factors for A; not accessed if EQUED = 'N'.  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.
[in,out]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          On entry, the N-by-NRHS righthand side matrix B.
          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
          B is overwritten by diag(S) * B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
          the original system of equations.  Note that if EQUED = 'Y',
          A and B are modified on exit, and the solution to the
          equilibrated system is inv(diag(S))*X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]RCOND
          RCOND is DOUBLE PRECISION
          The estimate of the reciprocal condition number of the matrix
          A after equilibration (if done).  If RCOND is less than the
          machine precision (in particular, if RCOND = 0), the matrix
          is singular to working precision.  This condition is
          indicated by a return code of INFO > 0.
[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
          > 0: if INFO = i, and i is
                <= N:  the leading minor of order i of A is
                       not positive definite, so the factorization
                       could not be completed, and the solution has not
                       been computed. RCOND = 0 is returned.
                = N+1: U is nonsingular, but RCOND is less than machine
                       precision, meaning that the matrix is singular
                       to working precision.  Nevertheless, the
                       solution and error bounds are computed because
                       there are a number of situations where the
                       computed solution can be more accurate than the
                       value of RCOND would suggest.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 308 of file zposvx.f.

308 *
309 * -- LAPACK driver routine (version 3.4.1) --
310 * -- LAPACK is a software package provided by Univ. of Tennessee, --
311 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
312 * April 2012
313 *
314 * .. Scalar Arguments ..
315  CHARACTER equed, fact, uplo
316  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs
317  DOUBLE PRECISION rcond
318 * ..
319 * .. Array Arguments ..
320  DOUBLE PRECISION berr( * ), ferr( * ), rwork( * ), s( * )
321  COMPLEX*16 a( lda, * ), af( ldaf, * ), b( ldb, * ),
322  $ work( * ), x( ldx, * )
323 * ..
324 *
325 * =====================================================================
326 *
327 * .. Parameters ..
328  DOUBLE PRECISION zero, one
329  parameter( zero = 0.0d+0, one = 1.0d+0 )
330 * ..
331 * .. Local Scalars ..
332  LOGICAL equil, nofact, rcequ
333  INTEGER i, infequ, j
334  DOUBLE PRECISION amax, anorm, bignum, scond, smax, smin, smlnum
335 * ..
336 * .. External Functions ..
337  LOGICAL lsame
338  DOUBLE PRECISION dlamch, zlanhe
339  EXTERNAL lsame, dlamch, zlanhe
340 * ..
341 * .. External Subroutines ..
342  EXTERNAL xerbla, zlacpy, zlaqhe, zpocon, zpoequ, zporfs,
343  $ zpotrf, zpotrs
344 * ..
345 * .. Intrinsic Functions ..
346  INTRINSIC max, min
347 * ..
348 * .. Executable Statements ..
349 *
350  info = 0
351  nofact = lsame( fact, 'N' )
352  equil = lsame( fact, 'E' )
353  IF( nofact .OR. equil ) THEN
354  equed = 'N'
355  rcequ = .false.
356  ELSE
357  rcequ = lsame( equed, 'Y' )
358  smlnum = dlamch( 'Safe minimum' )
359  bignum = one / smlnum
360  END IF
361 *
362 * Test the input parameters.
363 *
364  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
365  $ THEN
366  info = -1
367  ELSE IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) )
368  $ THEN
369  info = -2
370  ELSE IF( n.LT.0 ) THEN
371  info = -3
372  ELSE IF( nrhs.LT.0 ) THEN
373  info = -4
374  ELSE IF( lda.LT.max( 1, n ) ) THEN
375  info = -6
376  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
377  info = -8
378  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
379  $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
380  info = -9
381  ELSE
382  IF( rcequ ) THEN
383  smin = bignum
384  smax = zero
385  DO 10 j = 1, n
386  smin = min( smin, s( j ) )
387  smax = max( smax, s( j ) )
388  10 CONTINUE
389  IF( smin.LE.zero ) THEN
390  info = -10
391  ELSE IF( n.GT.0 ) THEN
392  scond = max( smin, smlnum ) / min( smax, bignum )
393  ELSE
394  scond = one
395  END IF
396  END IF
397  IF( info.EQ.0 ) THEN
398  IF( ldb.LT.max( 1, n ) ) THEN
399  info = -12
400  ELSE IF( ldx.LT.max( 1, n ) ) THEN
401  info = -14
402  END IF
403  END IF
404  END IF
405 *
406  IF( info.NE.0 ) THEN
407  CALL xerbla( 'ZPOSVX', -info )
408  RETURN
409  END IF
410 *
411  IF( equil ) THEN
412 *
413 * Compute row and column scalings to equilibrate the matrix A.
414 *
415  CALL zpoequ( n, a, lda, s, scond, amax, infequ )
416  IF( infequ.EQ.0 ) THEN
417 *
418 * Equilibrate the matrix.
419 *
420  CALL zlaqhe( uplo, n, a, lda, s, scond, amax, equed )
421  rcequ = lsame( equed, 'Y' )
422  END IF
423  END IF
424 *
425 * Scale the right hand side.
426 *
427  IF( rcequ ) THEN
428  DO 30 j = 1, nrhs
429  DO 20 i = 1, n
430  b( i, j ) = s( i )*b( i, j )
431  20 CONTINUE
432  30 CONTINUE
433  END IF
434 *
435  IF( nofact .OR. equil ) THEN
436 *
437 * Compute the Cholesky factorization A = U**H *U or A = L*L**H.
438 *
439  CALL zlacpy( uplo, n, n, a, lda, af, ldaf )
440  CALL zpotrf( uplo, n, af, ldaf, info )
441 *
442 * Return if INFO is non-zero.
443 *
444  IF( info.GT.0 )THEN
445  rcond = zero
446  RETURN
447  END IF
448  END IF
449 *
450 * Compute the norm of the matrix A.
451 *
452  anorm = zlanhe( '1', uplo, n, a, lda, rwork )
453 *
454 * Compute the reciprocal of the condition number of A.
455 *
456  CALL zpocon( uplo, n, af, ldaf, anorm, rcond, work, rwork, info )
457 *
458 * Compute the solution matrix X.
459 *
460  CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
461  CALL zpotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
462 *
463 * Use iterative refinement to improve the computed solution and
464 * compute error bounds and backward error estimates for it.
465 *
466  CALL zporfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,
467  $ ferr, berr, work, rwork, info )
468 *
469 * Transform the solution matrix X to a solution of the original
470 * system.
471 *
472  IF( rcequ ) THEN
473  DO 50 j = 1, nrhs
474  DO 40 i = 1, n
475  x( i, j ) = s( i )*x( i, j )
476  40 CONTINUE
477  50 CONTINUE
478  DO 60 j = 1, nrhs
479  ferr( j ) = ferr( j ) / scond
480  60 CONTINUE
481  END IF
482 *
483 * Set INFO = N+1 if the matrix is singular to working precision.
484 *
485  IF( rcond.LT.dlamch( 'Epsilon' ) )
486  $ info = n + 1
487 *
488  RETURN
489 *
490 * End of ZPOSVX
491 *
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
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 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 dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zporfs(UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZPORFS
Definition: zporfs.f:185
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zlaqhe(UPLO, N, A, LDA, S, SCOND, AMAX, EQUED)
ZLAQHE scales a Hermitian matrix.
Definition: zlaqhe.f:136
subroutine zpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
ZPOCON
Definition: zpocon.f:123
subroutine zpoequ(N, A, LDA, S, SCOND, AMAX, INFO)
ZPOEQU
Definition: zpoequ.f:115

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zposvxx ( character  FACT,
character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
character  EQUED,
double precision, dimension( * )  S,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
double precision  RCOND,
double precision  RPVGRW,
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 
)

ZPOSVXX computes the solution to system of linear equations A * X = B for PO matrices

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

Purpose:
    ZPOSVXX uses the Cholesky factorization A = U**T*U or A = L*L**T
    to compute the solution to a complex*16 system of linear equations
    A * X = B, where A is an N-by-N symmetric positive definite matrix
    and X and B are N-by-NRHS matrices.

    If requested, both normwise and maximum componentwise error bounds
    are returned. ZPOSVXX will return a solution with a tiny
    guaranteed error (O(eps) where eps is the working machine
    precision) unless the matrix is very ill-conditioned, in which
    case a warning is returned. Relevant condition numbers also are
    calculated and returned.

    ZPOSVXX accepts user-provided factorizations and equilibration
    factors; see the definitions of the FACT and EQUED options.
    Solving with refinement and using a factorization from a previous
    ZPOSVXX call will also produce a solution with either O(eps)
    errors or warnings, but we cannot make that claim for general
    user-provided factorizations and equilibration factors if they
    differ from what ZPOSVXX would itself produce.
Description:
    The following steps are performed:

    1. If FACT = 'E', double precision scaling factors are computed to equilibrate
    the system:

      diag(S)*A*diag(S)     *inv(diag(S))*X = diag(S)*B

    Whether or not the system will be equilibrated depends on the
    scaling of the matrix A, but if equilibration is used, A is
    overwritten by diag(S)*A*diag(S) and B by diag(S)*B.

    2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
    factor the matrix A (after equilibration if FACT = 'E') as
       A = U**T* U,  if UPLO = 'U', or
       A = L * L**T,  if UPLO = 'L',
    where U is an upper triangular matrix and L is a lower triangular
    matrix.

    3. If the leading i-by-i principal minor is not positive definite,
    then the routine returns with INFO = i. Otherwise, the factored
    form of A is used to estimate the condition number of the matrix
    A (see argument RCOND).  If the reciprocal of the condition number
    is less than machine precision, the routine still goes on to solve
    for X and compute error bounds as described below.

    4. The system of equations is solved for X using the factored form
    of A.

    5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
    the routine will use iterative refinement to try to get a small
    error and error bounds.  Refinement calculates the residual to at
    least twice the working precision.

    6. If equilibration was used, the matrix X is premultiplied by
    diag(S) so that it solves the original system before
    equilibration.
     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]FACT
          FACT is CHARACTER*1
     Specifies whether or not the factored form of the matrix A is
     supplied on entry, and if not, whether the matrix A should be
     equilibrated before it is factored.
       = 'F':  On entry, AF contains the factored form of A.
               If EQUED is not 'N', the matrix A has been
               equilibrated with scaling factors given by S.
               A and AF are not modified.
       = 'N':  The matrix A will be copied to AF and factored.
       = 'E':  The matrix A will be equilibrated if necessary, then
               copied to AF and factored.
[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 matrices B and X.  NRHS >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
     On entry, the symmetric matrix A, except if FACT = 'F' and EQUED =
     'Y', then A must contain the equilibrated matrix
     diag(S)*A*diag(S).  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.  A is
     not modified if FACT = 'F' or 'N', or if FACT = 'E' and EQUED =
     'N' on exit.

     On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
     diag(S)*A*diag(S).
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in,out]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
     If FACT = 'F', then AF is an input argument and on entry
     contains the triangular factor U or L from the Cholesky
     factorization A = U**T*U or A = L*L**T, in the same storage
     format as A.  If EQUED .ne. 'N', then AF is the factored
     form of the equilibrated matrix diag(S)*A*diag(S).

     If FACT = 'N', then AF is an output argument and on exit
     returns the triangular factor U or L from the Cholesky
     factorization A = U**T*U or A = L*L**T of the original
     matrix A.

     If FACT = 'E', then AF is an output argument and on exit
     returns the triangular factor U or L from the Cholesky
     factorization A = U**T*U or A = L*L**T of the equilibrated
     matrix A (see the description of A for the form of the
     equilibrated matrix).
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= max(1,N).
[in,out]EQUED
          EQUED is CHARACTER*1
     Specifies the form of equilibration that was done.
       = 'N':  No equilibration (always true if FACT = 'N').
       = 'Y':  Both row and column equilibration, i.e., A has been
               replaced by diag(S) * A * diag(S).
     EQUED is an input argument if FACT = 'F'; otherwise, it is an
     output argument.
[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,out]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
     On entry, the N-by-NRHS right hand side matrix B.
     On exit,
     if EQUED = 'N', B is not modified;
     if EQUED = 'Y', B is overwritten by diag(S)*B;
[in]LDB
          LDB is INTEGER
     The leading dimension of the array B.  LDB >= max(1,N).
[out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
     If INFO = 0, the N-by-NRHS solution matrix X to the original
     system of equations.  Note that A and B are modified on exit if
     EQUED .ne. 'N', and the solution to the equilibrated system is
     inv(diag(S))*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]RPVGRW
          RPVGRW is DOUBLE PRECISION
     Reciprocal pivot growth.  On exit, this contains the reciprocal
     pivot growth factor norm(A)/norm(U). The "max absolute element"
     norm is used.  If this is much less than 1, then 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. If factorization fails with
     0<INFO<=N, then this contains the reciprocal pivot growth factor
     for the leading INFO columns of A.
[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 extra-precise refinement algorithm.
              (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 495 of file zposvxx.f.

495 *
496 * -- LAPACK driver routine (version 3.4.1) --
497 * -- LAPACK is a software package provided by Univ. of Tennessee, --
498 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
499 * April 2012
500 *
501 * .. Scalar Arguments ..
502  CHARACTER equed, fact, uplo
503  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs, nparams,
504  $ n_err_bnds
505  DOUBLE PRECISION rcond, rpvgrw
506 * ..
507 * .. Array Arguments ..
508  COMPLEX*16 a( lda, * ), af( ldaf, * ), b( ldb, * ),
509  $ work( * ), x( ldx, * )
510  DOUBLE PRECISION s( * ), params( * ), berr( * ), rwork( * ),
511  $ err_bnds_norm( nrhs, * ),
512  $ err_bnds_comp( nrhs, * )
513 * ..
514 *
515 * ==================================================================
516 *
517 * .. Parameters ..
518  DOUBLE PRECISION zero, one
519  parameter( zero = 0.0d+0, one = 1.0d+0 )
520  INTEGER final_nrm_err_i, final_cmp_err_i, berr_i
521  INTEGER rcond_i, nrm_rcond_i, nrm_err_i, cmp_rcond_i
522  INTEGER cmp_err_i, piv_growth_i
523  parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
524  $ berr_i = 3 )
525  parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
526  parameter( cmp_rcond_i = 7, cmp_err_i = 8,
527  $ piv_growth_i = 9 )
528 * ..
529 * .. Local Scalars ..
530  LOGICAL equil, nofact, rcequ
531  INTEGER infequ, j
532  DOUBLE PRECISION amax, bignum, smin, smax, scond, smlnum
533 * ..
534 * .. External Functions ..
535  EXTERNAL lsame, dlamch, zla_porpvgrw
536  LOGICAL lsame
537  DOUBLE PRECISION dlamch, zla_porpvgrw
538 * ..
539 * .. External Subroutines ..
540  EXTERNAL zpocon, zpoequb, zpotrf, zpotrs, zlacpy,
542 * ..
543 * .. Intrinsic Functions ..
544  INTRINSIC max, min
545 * ..
546 * .. Executable Statements ..
547 *
548  info = 0
549  nofact = lsame( fact, 'N' )
550  equil = lsame( fact, 'E' )
551  smlnum = dlamch( 'Safe minimum' )
552  bignum = one / smlnum
553  IF( nofact .OR. equil ) THEN
554  equed = 'N'
555  rcequ = .false.
556  ELSE
557  rcequ = lsame( equed, 'Y' )
558  ENDIF
559 *
560 * Default is failure. If an input parameter is wrong or
561 * factorization fails, make everything look horrible. Only the
562 * pivot growth is set here, the rest is initialized in ZPORFSX.
563 *
564  rpvgrw = zero
565 *
566 * Test the input parameters. PARAMS is not tested until ZPORFSX.
567 *
568  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.
569  $ lsame( fact, 'F' ) ) THEN
570  info = -1
571  ELSE IF( .NOT.lsame( uplo, 'U' ) .AND.
572  $ .NOT.lsame( uplo, 'L' ) ) THEN
573  info = -2
574  ELSE IF( n.LT.0 ) THEN
575  info = -3
576  ELSE IF( nrhs.LT.0 ) THEN
577  info = -4
578  ELSE IF( lda.LT.max( 1, n ) ) THEN
579  info = -6
580  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
581  info = -8
582  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
583  $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
584  info = -9
585  ELSE
586  IF ( rcequ ) THEN
587  smin = bignum
588  smax = zero
589  DO 10 j = 1, n
590  smin = min( smin, s( j ) )
591  smax = max( smax, s( j ) )
592  10 CONTINUE
593  IF( smin.LE.zero ) THEN
594  info = -10
595  ELSE IF( n.GT.0 ) THEN
596  scond = max( smin, smlnum ) / min( smax, bignum )
597  ELSE
598  scond = one
599  END IF
600  END IF
601  IF( info.EQ.0 ) THEN
602  IF( ldb.LT.max( 1, n ) ) THEN
603  info = -12
604  ELSE IF( ldx.LT.max( 1, n ) ) THEN
605  info = -14
606  END IF
607  END IF
608  END IF
609 *
610  IF( info.NE.0 ) THEN
611  CALL xerbla( 'ZPOSVXX', -info )
612  RETURN
613  END IF
614 *
615  IF( equil ) THEN
616 *
617 * Compute row and column scalings to equilibrate the matrix A.
618 *
619  CALL zpoequb( n, a, lda, s, scond, amax, infequ )
620  IF( infequ.EQ.0 ) THEN
621 *
622 * Equilibrate the matrix.
623 *
624  CALL zlaqhe( uplo, n, a, lda, s, scond, amax, equed )
625  rcequ = lsame( equed, 'Y' )
626  END IF
627  END IF
628 *
629 * Scale the right-hand side.
630 *
631  IF( rcequ ) CALL zlascl2( n, nrhs, s, b, ldb )
632 *
633  IF( nofact .OR. equil ) THEN
634 *
635 * Compute the Cholesky factorization of A.
636 *
637  CALL zlacpy( uplo, n, n, a, lda, af, ldaf )
638  CALL zpotrf( uplo, n, af, ldaf, info )
639 *
640 * Return if INFO is non-zero.
641 *
642  IF( info.GT.0 ) THEN
643 *
644 * Pivot in column INFO is exactly 0
645 * Compute the reciprocal pivot growth factor of the
646 * leading rank-deficient INFO columns of A.
647 *
648  rpvgrw = zla_porpvgrw( uplo, n, a, lda, af, ldaf, rwork )
649  RETURN
650  END IF
651  END IF
652 *
653 * Compute the reciprocal pivot growth factor RPVGRW.
654 *
655  rpvgrw = zla_porpvgrw( uplo, n, a, lda, af, ldaf, rwork )
656 *
657 * Compute the solution matrix X.
658 *
659  CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
660  CALL zpotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
661 *
662 * Use iterative refinement to improve the computed solution and
663 * compute error bounds and backward error estimates for it.
664 *
665  CALL zporfsx( uplo, equed, n, nrhs, a, lda, af, ldaf,
666  $ s, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm,
667  $ err_bnds_comp, nparams, params, work, rwork, info )
668 
669 *
670 * Scale solutions.
671 *
672  IF ( rcequ ) THEN
673  CALL zlascl2( n, nrhs, s, x, ldx )
674  END IF
675 *
676  RETURN
677 *
678 * End of ZPOSVXX
679 *
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102
subroutine zlascl2(M, N, D, X, LDX)
ZLASCL2 performs diagonal scaling on a vector.
Definition: zlascl2.f:93
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
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 dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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
subroutine zlaqhe(UPLO, N, A, LDA, S, SCOND, AMAX, EQUED)
ZLAQHE scales a Hermitian matrix.
Definition: zlaqhe.f:136
subroutine zpoequb(N, A, LDA, S, SCOND, AMAX, INFO)
ZPOEQUB
Definition: zpoequb.f:115
subroutine zpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
ZPOCON
Definition: zpocon.f:123
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
Definition: zporfsx.f:395

Here is the call graph for this function:

Here is the caller graph for this function: