Contents
zgesvx - use the LU factorization to compute the solution to
a complex system of linear equations A * X = B,
SUBROUTINE ZGESVX(FACT, TRANSA, N, NRHS, A, LDA, AF, LDAF, IPIVOT,
EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
WORK2, INFO)
CHARACTER * 1 FACT, TRANSA, EQUED
DOUBLE COMPLEX A(LDA,*), AF(LDAF,*), B(LDB,*), X(LDX,*),
WORK(*)
INTEGER N, NRHS, LDA, LDAF, LDB, LDX, INFO
INTEGER IPIVOT(*)
DOUBLE PRECISION RCOND
DOUBLE PRECISION R(*), C(*), FERR(*), BERR(*), WORK2(*)
SUBROUTINE ZGESVX_64(FACT, TRANSA, N, NRHS, A, LDA, AF, LDAF, IPIVOT,
EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
WORK2, INFO)
CHARACTER * 1 FACT, TRANSA, EQUED
DOUBLE COMPLEX A(LDA,*), AF(LDAF,*), B(LDB,*), X(LDX,*),
WORK(*)
INTEGER*8 N, NRHS, LDA, LDAF, LDB, LDX, INFO
INTEGER*8 IPIVOT(*)
DOUBLE PRECISION RCOND
DOUBLE PRECISION R(*), C(*), FERR(*), BERR(*), WORK2(*)
F95 INTERFACE
SUBROUTINE GESVX(FACT, [TRANSA], [N], [NRHS], A, [LDA], AF, [LDAF],
IPIVOT, EQUED, R, C, B, [LDB], X, [LDX], RCOND, FERR,
BERR, [WORK], [WORK2], [INFO])
CHARACTER(LEN=1) :: FACT, TRANSA, EQUED
COMPLEX(8), DIMENSION(:) :: WORK
COMPLEX(8), DIMENSION(:,:) :: A, AF, B, X
INTEGER :: N, NRHS, LDA, LDAF, LDB, LDX, INFO
INTEGER, DIMENSION(:) :: IPIVOT
REAL(8) :: RCOND
REAL(8), DIMENSION(:) :: R, C, FERR, BERR, WORK2
SUBROUTINE GESVX_64(FACT, [TRANSA], [N], [NRHS], A, [LDA], AF, [LDAF],
IPIVOT, EQUED, R, C, B, [LDB], X, [LDX], RCOND, FERR,
BERR, [WORK], [WORK2], [INFO])
CHARACTER(LEN=1) :: FACT, TRANSA, EQUED
COMPLEX(8), DIMENSION(:) :: WORK
COMPLEX(8), DIMENSION(:,:) :: A, AF, B, X
INTEGER(8) :: N, NRHS, LDA, LDAF, LDB, LDX, INFO
INTEGER(8), DIMENSION(:) :: IPIVOT
REAL(8) :: RCOND
REAL(8), DIMENSION(:) :: R, C, FERR, BERR, WORK2
C INTERFACE
#include <sunperf.h>
void zgesvx(char fact, char transa, int n, int nrhs, doub-
lecomplex *a, int lda, doublecomplex *af, int
ldaf, int *ipivot, char equed, double *r, double
*c, doublecomplex *b, int ldb, doublecomplex *x,
int ldx, double *rcond, double *ferr, double
*berr, int *info);
void zgesvx_64(char fact, char transa, long n, long nrhs,
doublecomplex *a, long lda, doublecomplex *af,
long ldaf, long *ipivot, char equed, double *r,
double *c, doublecomplex *b, long ldb, doublecom-
plex *x, long ldx, double *rcond, double *ferr,
double *berr, long *info);
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.
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 tri-
angular
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 fac-
tored 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.
FACT (input)
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 IPIVOT contain
the factored form of A. If EQUED is not 'N', the
matrix A has been equilibrated with scaling fac-
tors given by R and C. A, AF, and IPIVOT 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.
TRANSA (input)
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)
N (input) The number of linear equations, i.e., the order of
the matrix A. N >= 0.
NRHS (input)
The number of right hand sides, i.e., the number
of columns of the matrices B and X. NRHS >= 0.
A (input/output)
On entry, the N-by-N matrix A. If FACT = 'F' and
EQUED is not 'N', then A must have been equili-
brated 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 fol-
lows: EQUED = 'R': A := diag(R) * A
EQUED = 'C': A := A * diag(C)
EQUED = 'B': A := diag(R) * A * diag(C).
LDA (input)
The leading dimension of the array A. LDA >=
max(1,N).
AF (input or output)
If FACT = 'F', then AF is an input argument and on
entry contains the factors L and U from the fac-
torization 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 fac-
torization 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 fac-
torization A = P*L*U of the equilibrated matrix A
(see the description of A for the form of the
equilibrated matrix).
LDAF (input)
The leading dimension of the array AF. LDAF >=
max(1,N).
IPIVOT (input or output)
If FACT = 'F', then IPIVOT 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
IPIVOT(i).
If FACT = 'N', then IPIVOT 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 IPIVOT is an output argument
and on exit contains the pivot indices from the
factorization A = P*L*U of the equilibrated matrix
A.
EQUED (input or output)
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 equili-
bration, i.e., A has been postmultiplied by
diag(C). = 'B': Both row and column equilibra-
tion, 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.
R (input or output)
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.
C (input or output)
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.
B (input/output)
On entry, the N-by-NRHS right hand side matrix B.
On exit, if EQUED = 'N', B is not modified; if
TRANSA = 'N' and EQUED = 'R' or 'B', B is
overwritten by diag(R)*B; if TRANSA = 'T' or 'C'
and EQUED = 'C' or 'B', B is overwritten by
diag(C)*B.
LDB (input)
The leading dimension of the array B. LDB >=
max(1,N).
X (output)
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 TRANSA = 'N' and EQUED
= 'C' or 'B', or inv(diag(R))*X if TRANSA = 'T' or
'C' and EQUED = 'R' or 'B'.
LDX (input)
The leading dimension of the array X. LDX >=
max(1,N).
RCOND (output)
The estimate of the reciprocal condition number of
the matrix A after equilibration (if done). If
RCOND is less than the machine precision (in par-
ticular, if RCOND = 0), the matrix is singular to
working precision. This condition is indicated by
a return code of INFO > 0.
FERR (output)
The estimated forward error bound for each solu-
tion 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 ele-
ment 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.
BERR (output)
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).
WORK (workspace)
dimension(2*N)
WORK2 (workspace)
dimension(2*N) On exit, WORK2(1) contains the
reciprocal pivot growth factor norm(A)/norm(U).
The "max absolute element" norm is used. If
WORK2(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 factori-
zation fails with 0<INFO<=N, then WORK2(1) con-
tains the reciprocal pivot growth factor for the
leading INFO columns of A.
INFO (output)
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an ille-
gal 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.