LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zgesvx()

subroutine zgesvx ( character fact,
character trans,
integer n,
integer nrhs,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldaf, * ) af,
integer ldaf,
integer, dimension( * ) ipiv,
character equed,
double precision, dimension( * ) r,
double precision, dimension( * ) c,
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 )

ZGESVX computes the solution to system of linear equations A * X = B for GE matrices

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

Purpose:
!>
!> ZGESVX uses the LU factorization to compute the solution to a complex
!> system of linear equations
!>    A * X = B,
!> where A is an N-by-N 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:
!>       TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
!>       TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
!>       TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*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(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
!>    or diag(C)*B (if TRANS = 'T' or 'C').
!>
!> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
!>    matrix A (after equilibration if FACT = 'E') as
!>       A = P * L * U,
!>    where P is a permutation matrix, L is a unit lower triangular
!>    matrix, and U is upper triangular.
!>
!> 3. If some U(i,i)=0, so that U is exactly singular, 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(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') 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 and IPIV contain the factored form of A.
!>                  If EQUED is not 'N', the matrix A has been
!>                  equilibrated with scaling factors given by R and C.
!>                  A, AF, and IPIV 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]TRANS
!>          TRANS is CHARACTER*1
!>          Specifies the form of the system of equations:
!>          = 'N':  A * X = B     (No transpose)
!>          = 'T':  A**T * X = B  (Transpose)
!>          = 'C':  A**H * X = B  (Conjugate transpose)
!>
[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 N-by-N matrix A.  If FACT = 'F' and EQUED is
!>          not 'N', then A must have been equilibrated by the scaling
!>          factors in R and/or C.  A is not modified if FACT = 'F' or
!>          'N', or if FACT = 'E' and EQUED = 'N' on exit.
!>
!>          On exit, if EQUED .ne. 'N', A is scaled as follows:
!>          EQUED = 'R':  A := diag(R) * A
!>          EQUED = 'C':  A := A * diag(C)
!>          EQUED = 'B':  A := diag(R) * A * diag(C).
!>
[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 factors L and U from the factorization
!>          A = P*L*U as computed by ZGETRF.  If EQUED .ne. 'N', then
!>          AF is the factored form of the equilibrated matrix A.
!>
!>          If FACT = 'N', then AF is an output argument and on exit
!>          returns the factors L and U from the factorization A = P*L*U
!>          of the original matrix A.
!>
!>          If FACT = 'E', then AF is an output argument and on exit
!>          returns the factors L and U from the factorization A = P*L*U
!>          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]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          If FACT = 'F', then IPIV is an input argument and on entry
!>          contains the pivot indices from the factorization A = P*L*U
!>          as computed by ZGETRF; row i of the matrix was interchanged
!>          with row IPIV(i).
!>
!>          If FACT = 'N', then IPIV is an output argument and on exit
!>          contains the pivot indices from the factorization A = P*L*U
!>          of the original matrix A.
!>
!>          If FACT = 'E', then IPIV is an output argument and on exit
!>          contains the pivot indices from the factorization A = P*L*U
!>          of the equilibrated matrix A.
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>          Specifies the form of equilibration that was done.
!>          = 'N':  No equilibration (always true if FACT = 'N').
!>          = 'R':  Row equilibration, i.e., A has been premultiplied by
!>                  diag(R).
!>          = 'C':  Column equilibration, i.e., A has been postmultiplied
!>                  by diag(C).
!>          = 'B':  Both row and column equilibration, i.e., A has been
!>                  replaced by diag(R) * A * diag(C).
!>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>          output argument.
!>
[in,out]R
!>          R is DOUBLE PRECISION array, dimension (N)
!>          The row scale factors for A.  If EQUED = 'R' or 'B', A is
!>          multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
!>          is not accessed.  R is an input argument if FACT = 'F';
!>          otherwise, R is an output argument.  If FACT = 'F' and
!>          EQUED = 'R' or 'B', each element of R must be positive.
!>
[in,out]C
!>          C is DOUBLE PRECISION array, dimension (N)
!>          The column scale factors for A.  If EQUED = 'C' or 'B', A is
!>          multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
!>          is not accessed.  C is an input argument if FACT = 'F';
!>          otherwise, C is an output argument.  If FACT = 'F' and
!>          EQUED = 'C' or 'B', each element of C must be positive.
!>
[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 TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
!>          diag(R)*B;
!>          if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
!>          overwritten by diag(C)*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 A and B are
!>          modified on exit if EQUED .ne. 'N', and the solution to the
!>          equilibrated system is inv(diag(C))*X if TRANS = 'N' and
!>          EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
!>          and EQUED = 'R' or 'B'.
!>
[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 (MAX(1,2*N))
!>          On exit, RWORK(1) contains the reciprocal pivot growth
!>          factor norm(A)/norm(U). The  norm is
!>          used. If RWORK(1) 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, condition
!>          estimator RCOND, and forward error bound FERR could be
!>          unreliable. If factorization fails with 0<INFO<=N, then
!>          RWORK(1) contains the reciprocal pivot growth factor for the
!>          leading INFO columns of A.
!>
[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:  U(i,i) 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+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.

Definition at line 345 of file zgesvx.f.

349*
350* -- LAPACK driver routine --
351* -- LAPACK is a software package provided by Univ. of Tennessee, --
352* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
353*
354* .. Scalar Arguments ..
355 CHARACTER EQUED, FACT, TRANS
356 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
357 DOUBLE PRECISION RCOND
358* ..
359* .. Array Arguments ..
360 INTEGER IPIV( * )
361 DOUBLE PRECISION BERR( * ), C( * ), FERR( * ), R( * ),
362 $ RWORK( * )
363 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
364 $ WORK( * ), X( LDX, * )
365* ..
366*
367* =====================================================================
368*
369* .. Parameters ..
370 DOUBLE PRECISION ZERO, ONE
371 parameter( zero = 0.0d+0, one = 1.0d+0 )
372* ..
373* .. Local Scalars ..
374 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
375 CHARACTER NORM
376 INTEGER I, INFEQU, J
377 DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
378 $ ROWCND, RPVGRW, SMLNUM
379* ..
380* .. External Functions ..
381 LOGICAL LSAME
382 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANTR
383 EXTERNAL lsame, dlamch, zlange, zlantr
384* ..
385* .. External Subroutines ..
386 EXTERNAL xerbla, zgecon, zgeequ, zgerfs, zgetrf,
387 $ zgetrs,
388 $ zlacpy, zlaqge
389* ..
390* .. Intrinsic Functions ..
391 INTRINSIC max, min
392* ..
393* .. Executable Statements ..
394*
395 info = 0
396 nofact = lsame( fact, 'N' )
397 equil = lsame( fact, 'E' )
398 notran = lsame( trans, 'N' )
399 IF( nofact .OR. equil ) THEN
400 equed = 'N'
401 rowequ = .false.
402 colequ = .false.
403 ELSE
404 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
405 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
406 smlnum = dlamch( 'Safe minimum' )
407 bignum = one / smlnum
408 END IF
409*
410* Test the input parameters.
411*
412 IF( .NOT.nofact .AND.
413 $ .NOT.equil .AND.
414 $ .NOT.lsame( fact, 'F' ) )
415 $ THEN
416 info = -1
417 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
418 $ lsame( trans, 'C' ) ) THEN
419 info = -2
420 ELSE IF( n.LT.0 ) THEN
421 info = -3
422 ELSE IF( nrhs.LT.0 ) THEN
423 info = -4
424 ELSE IF( lda.LT.max( 1, n ) ) THEN
425 info = -6
426 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
427 info = -8
428 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
429 $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
430 info = -10
431 ELSE
432 IF( rowequ ) THEN
433 rcmin = bignum
434 rcmax = zero
435 DO 10 j = 1, n
436 rcmin = min( rcmin, r( j ) )
437 rcmax = max( rcmax, r( j ) )
438 10 CONTINUE
439 IF( rcmin.LE.zero ) THEN
440 info = -11
441 ELSE IF( n.GT.0 ) THEN
442 rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
443 ELSE
444 rowcnd = one
445 END IF
446 END IF
447 IF( colequ .AND. info.EQ.0 ) THEN
448 rcmin = bignum
449 rcmax = zero
450 DO 20 j = 1, n
451 rcmin = min( rcmin, c( j ) )
452 rcmax = max( rcmax, c( j ) )
453 20 CONTINUE
454 IF( rcmin.LE.zero ) THEN
455 info = -12
456 ELSE IF( n.GT.0 ) THEN
457 colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
458 ELSE
459 colcnd = one
460 END IF
461 END IF
462 IF( info.EQ.0 ) THEN
463 IF( ldb.LT.max( 1, n ) ) THEN
464 info = -14
465 ELSE IF( ldx.LT.max( 1, n ) ) THEN
466 info = -16
467 END IF
468 END IF
469 END IF
470*
471 IF( info.NE.0 ) THEN
472 CALL xerbla( 'ZGESVX', -info )
473 RETURN
474 END IF
475*
476 IF( equil ) THEN
477*
478* Compute row and column scalings to equilibrate the matrix A.
479*
480 CALL zgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax,
481 $ infequ )
482 IF( infequ.EQ.0 ) THEN
483*
484* Equilibrate the matrix.
485*
486 CALL zlaqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,
487 $ equed )
488 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
489 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
490 END IF
491 END IF
492*
493* Scale the right hand side.
494*
495 IF( notran ) THEN
496 IF( rowequ ) THEN
497 DO 40 j = 1, nrhs
498 DO 30 i = 1, n
499 b( i, j ) = r( i )*b( i, j )
500 30 CONTINUE
501 40 CONTINUE
502 END IF
503 ELSE IF( colequ ) THEN
504 DO 60 j = 1, nrhs
505 DO 50 i = 1, n
506 b( i, j ) = c( i )*b( i, j )
507 50 CONTINUE
508 60 CONTINUE
509 END IF
510*
511 IF( nofact .OR. equil ) THEN
512*
513* Compute the LU factorization of A.
514*
515 CALL zlacpy( 'Full', n, n, a, lda, af, ldaf )
516 CALL zgetrf( n, n, af, ldaf, ipiv, info )
517*
518* Return if INFO is non-zero.
519*
520 IF( info.GT.0 ) THEN
521*
522* Compute the reciprocal pivot growth factor of the
523* leading rank-deficient INFO columns of A.
524*
525 rpvgrw = zlantr( 'M', 'U', 'N', info, info, af, ldaf,
526 $ rwork )
527 IF( rpvgrw.EQ.zero ) THEN
528 rpvgrw = one
529 ELSE
530 rpvgrw = zlange( 'M', n, info, a, lda, rwork ) /
531 $ rpvgrw
532 END IF
533 rwork( 1 ) = rpvgrw
534 rcond = zero
535 RETURN
536 END IF
537 END IF
538*
539* Compute the norm of the matrix A and the
540* reciprocal pivot growth factor RPVGRW.
541*
542 IF( notran ) THEN
543 norm = '1'
544 ELSE
545 norm = 'I'
546 END IF
547 anorm = zlange( norm, n, n, a, lda, rwork )
548 rpvgrw = zlantr( 'M', 'U', 'N', n, n, af, ldaf, rwork )
549 IF( rpvgrw.EQ.zero ) THEN
550 rpvgrw = one
551 ELSE
552 rpvgrw = zlange( 'M', n, n, a, lda, rwork ) / rpvgrw
553 END IF
554*
555* Compute the reciprocal of the condition number of A.
556*
557 CALL zgecon( norm, n, af, ldaf, anorm, rcond, work, rwork,
558 $ info )
559*
560* Compute the solution matrix X.
561*
562 CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
563 CALL zgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
564*
565* Use iterative refinement to improve the computed solution and
566* compute error bounds and backward error estimates for it.
567*
568 CALL zgerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,
569 $ ldx, ferr, berr, work, rwork, info )
570*
571* Transform the solution matrix X to a solution of the original
572* system.
573*
574 IF( notran ) THEN
575 IF( colequ ) THEN
576 DO 80 j = 1, nrhs
577 DO 70 i = 1, n
578 x( i, j ) = c( i )*x( i, j )
579 70 CONTINUE
580 80 CONTINUE
581 DO 90 j = 1, nrhs
582 ferr( j ) = ferr( j ) / colcnd
583 90 CONTINUE
584 END IF
585 ELSE IF( rowequ ) THEN
586 DO 110 j = 1, nrhs
587 DO 100 i = 1, n
588 x( i, j ) = r( i )*x( i, j )
589 100 CONTINUE
590 110 CONTINUE
591 DO 120 j = 1, nrhs
592 ferr( j ) = ferr( j ) / rowcnd
593 120 CONTINUE
594 END IF
595*
596* Set INFO = N+1 if the matrix is singular to working precision.
597*
598 IF( rcond.LT.dlamch( 'Epsilon' ) )
599 $ info = n + 1
600*
601 rwork( 1 ) = rpvgrw
602 RETURN
603*
604* End of ZGESVX
605*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgecon
subroutine zgecon(norm, n, a, lda, anorm, rcond, work, rwork, info)
ZGECON
Definition zgecon.f:130
zgeequ
subroutine zgeequ(m, n, a, lda, r, c, rowcnd, colcnd, amax, info)
ZGEEQU
Definition zgeequ.f:138
zgerfs
subroutine zgerfs(trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZGERFS
Definition zgerfs.f:185
zgetrf
subroutine zgetrf(m, n, a, lda, ipiv, info)
ZGETRF
Definition zgetrf.f:106
zgetrs
subroutine zgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
ZGETRS
Definition zgetrs.f:119
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlange
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:113
zlantr
double precision function zlantr(norm, uplo, diag, m, n, a, lda, work)
ZLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlantr.f:141
zlaqge
subroutine zlaqge(m, n, a, lda, r, c, rowcnd, colcnd, amax, equed)
ZLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ.
Definition zlaqge.f:141
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: