 LAPACK 3.3.1 Linear Algebra PACKage

# dgesvx.f

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