 LAPACK 3.3.1 Linear Algebra PACKage

# cgesvx.f

Go to the documentation of this file.
```00001       SUBROUTINE CGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
00002      \$                   EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR,
00003      \$                   WORK, RWORK, INFO )
00004 *
00005 *  -- LAPACK driver routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          EQUED, FACT, TRANS
00012       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
00013       REAL               RCOND
00014 *     ..
00015 *     .. Array Arguments ..
00016       INTEGER            IPIV( * )
00017       REAL               BERR( * ), C( * ), FERR( * ), R( * ),
00018      \$                   RWORK( * )
00019       COMPLEX            A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00020      \$                   WORK( * ), X( LDX, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  CGESVX uses the LU factorization to compute the solution to a complex
00027 *  system of linear equations
00028 *     A * X = B,
00029 *  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
00030 *
00031 *  Error bounds on the solution and a condition estimate are also
00032 *  provided.
00033 *
00034 *  Description
00035 *  ===========
00036 *
00037 *  The following steps are performed:
00038 *
00039 *  1. If FACT = 'E', real scaling factors are computed to equilibrate
00040 *     the system:
00041 *        TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
00042 *        TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
00043 *        TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
00044 *     Whether or not the system will be equilibrated depends on the
00045 *     scaling of the matrix A, but if equilibration is used, A is
00046 *     overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
00047 *     or diag(C)*B (if TRANS = 'T' or 'C').
00048 *
00049 *  2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
00050 *     matrix A (after equilibration if FACT = 'E') as
00051 *        A = P * L * U,
00052 *     where P is a permutation matrix, L is a unit lower triangular
00053 *     matrix, and U is upper triangular.
00054 *
00055 *  3. If some U(i,i)=0, so that U is exactly singular, then the routine
00056 *     returns with INFO = i. Otherwise, the factored form of A is used
00057 *     to estimate the condition number of the matrix A.  If the
00058 *     reciprocal of the condition number is less than machine precision,
00059 *     INFO = N+1 is returned as a warning, but the routine still goes on
00060 *     to solve for X and compute error bounds as described below.
00061 *
00062 *  4. The system of equations is solved for X using the factored form
00063 *     of A.
00064 *
00065 *  5. Iterative refinement is applied to improve the computed solution
00066 *     matrix and calculate error bounds and backward error estimates
00067 *     for it.
00068 *
00069 *  6. If equilibration was used, the matrix X is premultiplied by
00070 *     diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
00071 *     that it solves the original system before equilibration.
00072 *
00073 *  Arguments
00074 *  =========
00075 *
00076 *  FACT    (input) CHARACTER*1
00077 *          Specifies whether or not the factored form of the matrix A is
00078 *          supplied on entry, and if not, whether the matrix A should be
00079 *          equilibrated before it is factored.
00080 *          = 'F':  On entry, AF and IPIV contain the factored form of A.
00081 *                  If EQUED is not 'N', the matrix A has been
00082 *                  equilibrated with scaling factors given by R and C.
00083 *                  A, AF, and IPIV are not modified.
00084 *          = 'N':  The matrix A will be copied to AF and factored.
00085 *          = 'E':  The matrix A will be equilibrated if necessary, then
00086 *                  copied to AF and factored.
00087 *
00088 *  TRANS   (input) CHARACTER*1
00089 *          Specifies the form of the system of equations:
00090 *          = 'N':  A * X = B     (No transpose)
00091 *          = 'T':  A**T * X = B  (Transpose)
00092 *          = 'C':  A**H * X = B  (Conjugate transpose)
00093 *
00094 *  N       (input) INTEGER
00095 *          The number of linear equations, i.e., the order of the
00096 *          matrix A.  N >= 0.
00097 *
00098 *  NRHS    (input) INTEGER
00099 *          The number of right hand sides, i.e., the number of columns
00100 *          of the matrices B and X.  NRHS >= 0.
00101 *
00102 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00103 *          On entry, the N-by-N matrix A.  If FACT = 'F' and EQUED is
00104 *          not 'N', then A must have been equilibrated by the scaling
00105 *          factors in R and/or C.  A is not modified if FACT = 'F' or
00106 *          'N', or if FACT = 'E' and EQUED = 'N' on exit.
00107 *
00108 *          On exit, if EQUED .ne. 'N', A is scaled as follows:
00109 *          EQUED = 'R':  A := diag(R) * A
00110 *          EQUED = 'C':  A := A * diag(C)
00111 *          EQUED = 'B':  A := diag(R) * A * diag(C).
00112 *
00113 *  LDA     (input) INTEGER
00114 *          The leading dimension of the array A.  LDA >= max(1,N).
00115 *
00116 *  AF      (input or output) COMPLEX array, dimension (LDAF,N)
00117 *          If FACT = 'F', then AF is an input argument and on entry
00118 *          contains the factors L and U from the factorization
00119 *          A = P*L*U as computed by CGETRF.  If EQUED .ne. 'N', then
00120 *          AF is the factored form of the equilibrated matrix A.
00121 *
00122 *          If FACT = 'N', then AF is an output argument and on exit
00123 *          returns the factors L and U from the factorization A = P*L*U
00124 *          of the original matrix A.
00125 *
00126 *          If FACT = 'E', then AF is an output argument and on exit
00127 *          returns the factors L and U from the factorization A = P*L*U
00128 *          of the equilibrated matrix A (see the description of A for
00129 *          the form of the equilibrated matrix).
00130 *
00131 *  LDAF    (input) INTEGER
00132 *          The leading dimension of the array AF.  LDAF >= max(1,N).
00133 *
00134 *  IPIV    (input or output) INTEGER array, dimension (N)
00135 *          If FACT = 'F', then IPIV is an input argument and on entry
00136 *          contains the pivot indices from the factorization A = P*L*U
00137 *          as computed by CGETRF; row i of the matrix was interchanged
00138 *          with row IPIV(i).
00139 *
00140 *          If FACT = 'N', then IPIV is an output argument and on exit
00141 *          contains the pivot indices from the factorization A = P*L*U
00142 *          of the original matrix A.
00143 *
00144 *          If FACT = 'E', then IPIV is an output argument and on exit
00145 *          contains the pivot indices from the factorization A = P*L*U
00146 *          of the equilibrated matrix A.
00147 *
00148 *  EQUED   (input or output) CHARACTER*1
00149 *          Specifies the form of equilibration that was done.
00150 *          = 'N':  No equilibration (always true if FACT = 'N').
00151 *          = 'R':  Row equilibration, i.e., A has been premultiplied by
00152 *                  diag(R).
00153 *          = 'C':  Column equilibration, i.e., A has been postmultiplied
00154 *                  by diag(C).
00155 *          = 'B':  Both row and column equilibration, i.e., A has been
00156 *                  replaced by diag(R) * A * diag(C).
00157 *          EQUED is an input argument if FACT = 'F'; otherwise, it is an
00158 *          output argument.
00159 *
00160 *  R       (input or output) REAL array, dimension (N)
00161 *          The row scale factors for A.  If EQUED = 'R' or 'B', A is
00162 *          multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
00163 *          is not accessed.  R is an input argument if FACT = 'F';
00164 *          otherwise, R is an output argument.  If FACT = 'F' and
00165 *          EQUED = 'R' or 'B', each element of R must be positive.
00166 *
00167 *  C       (input or output) REAL array, dimension (N)
00168 *          The column scale factors for A.  If EQUED = 'C' or 'B', A is
00169 *          multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
00170 *          is not accessed.  C is an input argument if FACT = 'F';
00171 *          otherwise, C is an output argument.  If FACT = 'F' and
00172 *          EQUED = 'C' or 'B', each element of C must be positive.
00173 *
00174 *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
00175 *          On entry, the N-by-NRHS right hand side matrix B.
00176 *          On exit,
00177 *          if EQUED = 'N', B is not modified;
00178 *          if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
00179 *          diag(R)*B;
00180 *          if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
00181 *          overwritten by diag(C)*B.
00182 *
00183 *  LDB     (input) INTEGER
00184 *          The leading dimension of the array B.  LDB >= max(1,N).
00185 *
00186 *  X       (output) COMPLEX array, dimension (LDX,NRHS)
00187 *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X
00188 *          to the original system of equations.  Note that A and B are
00189 *          modified on exit if EQUED .ne. 'N', and the solution to the
00190 *          equilibrated system is inv(diag(C))*X if TRANS = 'N' and
00191 *          EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
00192 *          and EQUED = 'R' or 'B'.
00193 *
00194 *  LDX     (input) INTEGER
00195 *          The leading dimension of the array X.  LDX >= max(1,N).
00196 *
00197 *  RCOND   (output) REAL
00198 *          The estimate of the reciprocal condition number of the matrix
00199 *          A after equilibration (if done).  If RCOND is less than the
00200 *          machine precision (in particular, if RCOND = 0), the matrix
00201 *          is singular to working precision.  This condition is
00202 *          indicated by a return code of INFO > 0.
00203 *
00204 *  FERR    (output) REAL array, dimension (NRHS)
00205 *          The estimated forward error bound for each solution vector
00206 *          X(j) (the j-th column of the solution matrix X).
00207 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00208 *          is an estimated upper bound for the magnitude of the largest
00209 *          element in (X(j) - XTRUE) divided by the magnitude of the
00210 *          largest element in X(j).  The estimate is as reliable as
00211 *          the estimate for RCOND, and is almost always a slight
00212 *          overestimate of the true error.
00213 *
00214 *  BERR    (output) REAL array, dimension (NRHS)
00215 *          The componentwise relative backward error of each solution
00216 *          vector X(j) (i.e., the smallest relative change in
00217 *          any element of A or B that makes X(j) an exact solution).
00218 *
00219 *  WORK    (workspace) COMPLEX array, dimension (2*N)
00220 *
00221 *  RWORK   (workspace/output) REAL array, dimension (2*N)
00222 *          On exit, RWORK(1) contains the reciprocal pivot growth
00223 *          factor norm(A)/norm(U). The "max absolute element" norm is
00224 *          used. If RWORK(1) is much less than 1, then the stability
00225 *          of the LU factorization of the (equilibrated) matrix A
00226 *          could be poor. This also means that the solution X, condition
00227 *          estimator RCOND, and forward error bound FERR could be
00228 *          unreliable. If factorization fails with 0<INFO<=N, then
00229 *          RWORK(1) contains the reciprocal pivot growth factor for the
00230 *          leading INFO columns of A.
00231 *
00232 *  INFO    (output) INTEGER
00233 *          = 0:  successful exit
00234 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00235 *          > 0:  if INFO = i, and i is
00236 *                <= N:  U(i,i) is exactly zero.  The factorization has
00237 *                       been completed, but the factor U is exactly
00238 *                       singular, so the solution and error bounds
00239 *                       could not be computed. RCOND = 0 is returned.
00240 *                = N+1: U is nonsingular, but RCOND is less than machine
00241 *                       precision, meaning that the matrix is singular
00242 *                       to working precision.  Nevertheless, the
00243 *                       solution and error bounds are computed because
00244 *                       there are a number of situations where the
00245 *                       computed solution can be more accurate than the
00246 *                       value of RCOND would suggest.
00247 *
00248 *  =====================================================================
00249 *
00250 *     .. Parameters ..
00251       REAL               ZERO, ONE
00252       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00253 *     ..
00254 *     .. Local Scalars ..
00255       LOGICAL            COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
00256       CHARACTER          NORM
00257       INTEGER            I, INFEQU, J
00258       REAL               AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
00259      \$                   ROWCND, RPVGRW, SMLNUM
00260 *     ..
00261 *     .. External Functions ..
00262       LOGICAL            LSAME
00263       REAL               CLANGE, CLANTR, SLAMCH
00264       EXTERNAL           LSAME, CLANGE, CLANTR, SLAMCH
00265 *     ..
00266 *     .. External Subroutines ..
00267       EXTERNAL           CGECON, CGEEQU, CGERFS, CGETRF, CGETRS, CLACPY,
00268      \$                   CLAQGE, XERBLA
00269 *     ..
00270 *     .. Intrinsic Functions ..
00271       INTRINSIC          MAX, MIN
00272 *     ..
00273 *     .. Executable Statements ..
00274 *
00275       INFO = 0
00276       NOFACT = LSAME( FACT, 'N' )
00277       EQUIL = LSAME( FACT, 'E' )
00278       NOTRAN = LSAME( TRANS, 'N' )
00279       IF( NOFACT .OR. EQUIL ) THEN
00280          EQUED = 'N'
00281          ROWEQU = .FALSE.
00282          COLEQU = .FALSE.
00283       ELSE
00284          ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
00285          COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
00286          SMLNUM = SLAMCH( 'Safe minimum' )
00287          BIGNUM = ONE / SMLNUM
00288       END IF
00289 *
00290 *     Test the input parameters.
00291 *
00292       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
00293      \$     THEN
00294          INFO = -1
00295       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00296      \$         LSAME( TRANS, 'C' ) ) THEN
00297          INFO = -2
00298       ELSE IF( N.LT.0 ) THEN
00299          INFO = -3
00300       ELSE IF( NRHS.LT.0 ) THEN
00301          INFO = -4
00302       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00303          INFO = -6
00304       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
00305          INFO = -8
00306       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
00307      \$         ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
00308          INFO = -10
00309       ELSE
00310          IF( ROWEQU ) THEN
00311             RCMIN = BIGNUM
00312             RCMAX = ZERO
00313             DO 10 J = 1, N
00314                RCMIN = MIN( RCMIN, R( J ) )
00315                RCMAX = MAX( RCMAX, R( J ) )
00316    10       CONTINUE
00317             IF( RCMIN.LE.ZERO ) THEN
00318                INFO = -11
00319             ELSE IF( N.GT.0 ) THEN
00320                ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
00321             ELSE
00322                ROWCND = ONE
00323             END IF
00324          END IF
00325          IF( COLEQU .AND. INFO.EQ.0 ) THEN
00326             RCMIN = BIGNUM
00327             RCMAX = ZERO
00328             DO 20 J = 1, N
00329                RCMIN = MIN( RCMIN, C( J ) )
00330                RCMAX = MAX( RCMAX, C( J ) )
00331    20       CONTINUE
00332             IF( RCMIN.LE.ZERO ) THEN
00333                INFO = -12
00334             ELSE IF( N.GT.0 ) THEN
00335                COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
00336             ELSE
00337                COLCND = ONE
00338             END IF
00339          END IF
00340          IF( INFO.EQ.0 ) THEN
00341             IF( LDB.LT.MAX( 1, N ) ) THEN
00342                INFO = -14
00343             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00344                INFO = -16
00345             END IF
00346          END IF
00347       END IF
00348 *
00349       IF( INFO.NE.0 ) THEN
00350          CALL XERBLA( 'CGESVX', -INFO )
00351          RETURN
00352       END IF
00353 *
00354       IF( EQUIL ) THEN
00355 *
00356 *        Compute row and column scalings to equilibrate the matrix A.
00357 *
00358          CALL CGEEQU( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFEQU )
00359          IF( INFEQU.EQ.0 ) THEN
00360 *
00361 *           Equilibrate the matrix.
00362 *
00363             CALL CLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
00364      \$                   EQUED )
00365             ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
00366             COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
00367          END IF
00368       END IF
00369 *
00370 *     Scale the right hand side.
00371 *
00372       IF( NOTRAN ) THEN
00373          IF( ROWEQU ) THEN
00374             DO 40 J = 1, NRHS
00375                DO 30 I = 1, N
00376                   B( I, J ) = R( I )*B( I, J )
00377    30          CONTINUE
00378    40       CONTINUE
00379          END IF
00380       ELSE IF( COLEQU ) THEN
00381          DO 60 J = 1, NRHS
00382             DO 50 I = 1, N
00383                B( I, J ) = C( I )*B( I, J )
00384    50       CONTINUE
00385    60    CONTINUE
00386       END IF
00387 *
00388       IF( NOFACT .OR. EQUIL ) THEN
00389 *
00390 *        Compute the LU factorization of A.
00391 *
00392          CALL CLACPY( 'Full', N, N, A, LDA, AF, LDAF )
00393          CALL CGETRF( N, N, AF, LDAF, IPIV, INFO )
00394 *
00395 *        Return if INFO is non-zero.
00396 *
00397          IF( INFO.GT.0 ) THEN
00398 *
00399 *           Compute the reciprocal pivot growth factor of the
00400 *           leading rank-deficient INFO columns of A.
00401 *
00402             RPVGRW = CLANTR( 'M', 'U', 'N', INFO, INFO, AF, LDAF,
00403      \$               RWORK )
00404             IF( RPVGRW.EQ.ZERO ) THEN
00405                RPVGRW = ONE
00406             ELSE
00407                RPVGRW = CLANGE( 'M', N, INFO, A, LDA, RWORK ) /
00408      \$                  RPVGRW
00409             END IF
00410             RWORK( 1 ) = RPVGRW
00411             RCOND = ZERO
00412             RETURN
00413          END IF
00414       END IF
00415 *
00416 *     Compute the norm of the matrix A and the
00417 *     reciprocal pivot growth factor RPVGRW.
00418 *
00419       IF( NOTRAN ) THEN
00420          NORM = '1'
00421       ELSE
00422          NORM = 'I'
00423       END IF
00424       ANORM = CLANGE( NORM, N, N, A, LDA, RWORK )
00425       RPVGRW = CLANTR( 'M', 'U', 'N', N, N, AF, LDAF, RWORK )
00426       IF( RPVGRW.EQ.ZERO ) THEN
00427          RPVGRW = ONE
00428       ELSE
00429          RPVGRW = CLANGE( 'M', N, N, A, LDA, RWORK ) / RPVGRW
00430       END IF
00431 *
00432 *     Compute the reciprocal of the condition number of A.
00433 *
00434       CALL CGECON( NORM, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO )
00435 *
00436 *     Compute the solution matrix X.
00437 *
00438       CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00439       CALL CGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
00440 *
00441 *     Use iterative refinement to improve the computed solution and
00442 *     compute error bounds and backward error estimates for it.
00443 *
00444       CALL CGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
00445      \$             LDX, FERR, BERR, WORK, RWORK, INFO )
00446 *
00447 *     Transform the solution matrix X to a solution of the original
00448 *     system.
00449 *
00450       IF( NOTRAN ) THEN
00451          IF( COLEQU ) THEN
00452             DO 80 J = 1, NRHS
00453                DO 70 I = 1, N
00454                   X( I, J ) = C( I )*X( I, J )
00455    70          CONTINUE
00456    80       CONTINUE
00457             DO 90 J = 1, NRHS
00458                FERR( J ) = FERR( J ) / COLCND
00459    90       CONTINUE
00460          END IF
00461       ELSE IF( ROWEQU ) THEN
00462          DO 110 J = 1, NRHS
00463             DO 100 I = 1, N
00464                X( I, J ) = R( I )*X( I, J )
00465   100       CONTINUE
00466   110    CONTINUE
00467          DO 120 J = 1, NRHS
00468             FERR( J ) = FERR( J ) / ROWCND
00469   120    CONTINUE
00470       END IF
00471 *
00472 *     Set INFO = N+1 if the matrix is singular to working precision.
00473 *
00474       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
00475      \$   INFO = N + 1
00476 *
00477       RWORK( 1 ) = RPVGRW
00478       RETURN
00479 *
00480 *     End of CGESVX
00481 *
00482       END
```