LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zposvx()

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.7.0) --
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 zpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
ZPOTRS
Definition: zpotrs.f:112
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zporfs(UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZPORFS
Definition: zporfs.f:185
subroutine zpoequ(N, A, LDA, S, SCOND, AMAX, INFO)
ZPOEQU
Definition: zpoequ.f:115
subroutine zpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
ZPOCON
Definition: zpocon.f:123
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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 zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102
subroutine zlaqhe(UPLO, N, A, LDA, S, SCOND, AMAX, EQUED)
ZLAQHE scales a Hermitian matrix.
Definition: zlaqhe.f:136
Here is the call graph for this function:
Here is the caller graph for this function: