Return to dsgesv.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE DSGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
    2:      +                   SWORK, ITER, INFO )
    3: *
    4: *  -- LAPACK PROTOTYPE driver routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     February 2007
    8: *
    9: *     ..
   10: *     .. Scalar Arguments ..
   11:       INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
   12: *     ..
   13: *     .. Array Arguments ..
   14:       INTEGER            IPIV( * )
   15:       REAL               SWORK( * )
   16:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( N, * ),
   17:      +                   X( LDX, * )
   18: *     ..
   19: *
   20: *  Purpose
   21: *  =======
   22: *
   23: *  DSGESV computes the solution to a real system of linear equations
   24: *     A * X = B,
   25: *  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
   26: *
   27: *  DSGESV first attempts to factorize the matrix in SINGLE PRECISION
   28: *  and use this factorization within an iterative refinement procedure
   29: *  to produce a solution with DOUBLE PRECISION normwise backward error
   30: *  quality (see below). If the approach fails the method switches to a
   31: *  DOUBLE PRECISION factorization and solve.
   32: *
   33: *  The iterative refinement is not going to be a winning strategy if
   34: *  the ratio SINGLE PRECISION performance over DOUBLE PRECISION
   35: *  performance is too small. A reasonable strategy should take the
   36: *  number of right-hand sides and the size of the matrix into account.
   37: *  This might be done with a call to ILAENV in the future. Up to now, we
   38: *  always try iterative refinement.
   39: *
   40: *  The iterative refinement process is stopped if
   41: *      ITER > ITERMAX
   42: *  or for all the RHS we have:
   43: *      RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
   44: *  where
   45: *      o ITER is the number of the current iteration in the iterative
   46: *        refinement process
   47: *      o RNRM is the infinity-norm of the residual
   48: *      o XNRM is the infinity-norm of the solution
   49: *      o ANRM is the infinity-operator-norm of the matrix A
   50: *      o EPS is the machine epsilon returned by DLAMCH('Epsilon')
   51: *  The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
   52: *  respectively.
   53: *
   54: *  Arguments
   55: *  =========
   56: *
   57: *  N       (input) INTEGER
   58: *          The number of linear equations, i.e., the order of the
   59: *          matrix A.  N >= 0.
   60: *
   61: *  NRHS    (input) INTEGER
   62: *          The number of right hand sides, i.e., the number of columns
   63: *          of the matrix B.  NRHS >= 0.
   64: *
   65: *  A       (input or input/ouptut) DOUBLE PRECISION array,
   66: *          dimension (LDA,N)
   67: *          On entry, the N-by-N coefficient matrix A.
   68: *          On exit, if iterative refinement has been successfully used
   69: *          (INFO.EQ.0 and ITER.GE.0, see description below), then A is
   70: *          unchanged, if double precision factorization has been used
   71: *          (INFO.EQ.0 and ITER.LT.0, see description below), then the
   72: *          array A contains the factors L and U from the factorization
   73: *          A = P*L*U; the unit diagonal elements of L are not stored.
   74: *
   75: *  LDA     (input) INTEGER
   76: *          The leading dimension of the array A.  LDA >= max(1,N).
   77: *
   78: *  IPIV    (output) INTEGER array, dimension (N)
   79: *          The pivot indices that define the permutation matrix P;
   80: *          row i of the matrix was interchanged with row IPIV(i).
   81: *          Corresponds either to the single precision factorization
   82: *          (if INFO.EQ.0 and ITER.GE.0) or the double precision
   83: *          factorization (if INFO.EQ.0 and ITER.LT.0).
   84: *
   85: *  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
   86: *          The N-by-NRHS right hand side matrix B.
   87: *
   88: *  LDB     (input) INTEGER
   89: *          The leading dimension of the array B.  LDB >= max(1,N).
   90: *
   91: *  X       (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
   92: *          If INFO = 0, the N-by-NRHS solution matrix X.
   93: *
   94: *  LDX     (input) INTEGER
   95: *          The leading dimension of the array X.  LDX >= max(1,N).
   96: *
   97: *  WORK    (workspace) DOUBLE PRECISION array, dimension (N*NRHS)
   98: *          This array is used to hold the residual vectors.
   99: *
  100: *  SWORK   (workspace) REAL array, dimension (N*(N+NRHS))
  101: *          This array is used to use the single precision matrix and the
  102: *          right-hand sides or solutions in single precision.
  103: *
  104: *  ITER    (output) INTEGER
  105: *          < 0: iterative refinement has failed, double precision
  106: *               factorization has been performed
  107: *               -1 : the routine fell back to full precision for
  108: *                    implementation- or machine-specific reasons
  109: *               -2 : narrowing the precision induced an overflow,
  110: *                    the routine fell back to full precision
  111: *               -3 : failure of SGETRF
  112: *               -31: stop the iterative refinement after the 30th
  113: *                    iterations
  114: *          > 0: iterative refinement has been sucessfully used.
  115: *               Returns the number of iterations
  116: *
  117: *  INFO    (output) INTEGER
  118: *          = 0:  successful exit
  119: *          < 0:  if INFO = -i, the i-th argument had an illegal value
  120: *          > 0:  if INFO = i, U(i,i) computed in DOUBLE PRECISION is
  121: *                exactly zero.  The factorization has been completed,
  122: *                but the factor U is exactly singular, so the solution
  123: *                could not be computed.
  124: *
  125: *  =========
  126: *
  127: *     .. Parameters ..
  128:       LOGICAL            DOITREF
  129:       PARAMETER          ( DOITREF = .TRUE. )
  130: *
  131:       INTEGER            ITERMAX
  132:       PARAMETER          ( ITERMAX = 30 )
  133: *
  134:       DOUBLE PRECISION   BWDMAX
  135:       PARAMETER          ( BWDMAX = 1.0E+00 )
  136: *
  137:       DOUBLE PRECISION   NEGONE, ONE
  138:       PARAMETER          ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
  139: *
  140: *     .. Local Scalars ..
  141:       INTEGER            I, IITER, PTSA, PTSX
  142:       DOUBLE PRECISION   ANRM, CTE, EPS, RNRM, XNRM
  143: *
  144: *     .. External Subroutines ..
  145:       EXTERNAL           DAXPY, DGEMM, DLACPY, DLAG2S, SLAG2D, SGETRF,
  146:      +                   SGETRS, XERBLA
  147: *     ..
  148: *     .. External Functions ..
  149:       INTEGER            IDAMAX
  150:       DOUBLE PRECISION   DLAMCH, DLANGE
  151:       EXTERNAL           IDAMAX, DLAMCH, DLANGE
  152: *     ..
  153: *     .. Intrinsic Functions ..
  154:       INTRINSIC          ABS, DBLE, MAX, SQRT
  155: *     ..
  156: *     .. Executable Statements ..
  157: *
  158:       INFO = 0
  159:       ITER = 0
  160: *
  161: *     Test the input parameters.
  162: *
  163:       IF( N.LT.0 ) THEN
  164:          INFO = -1
  165:       ELSE IF( NRHS.LT.0 ) THEN
  166:          INFO = -2
  167:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  168:          INFO = -4
  169:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  170:          INFO = -7
  171:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
  172:          INFO = -9
  173:       END IF
  174:       IF( INFO.NE.0 ) THEN
  175:          CALL XERBLA( 'DSGESV', -INFO )
  176:          RETURN
  177:       END IF
  178: *
  179: *     Quick return if (N.EQ.0).
  180: *
  181:       IF( N.EQ.0 )
  182:      +   RETURN
  183: *
  184: *     Skip single precision iterative refinement if a priori slower
  185: *     than double precision factorization.
  186: *
  187:       IF( .NOT.DOITREF ) THEN
  188:          ITER = -1
  189:          GO TO 40
  190:       END IF
  191: *
  192: *     Compute some constants.
  193: *
  194:       ANRM = DLANGE( 'I', N, N, A, LDA, WORK )
  195:       EPS = DLAMCH( 'Epsilon' )
  196:       CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
  197: *
  198: *     Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
  199: *
  200:       PTSA = 1
  201:       PTSX = PTSA + N*N
  202: *
  203: *     Convert B from double precision to single precision and store the
  204: *     result in SX.
  205: *
  206:       CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
  207: *
  208:       IF( INFO.NE.0 ) THEN
  209:          ITER = -2
  210:          GO TO 40
  211:       END IF
  212: *
  213: *     Convert A from double precision to single precision and store the
  214: *     result in SA.
  215: *
  216:       CALL DLAG2S( N, N, A, LDA, SWORK( PTSA ), N, INFO )
  217: *
  218:       IF( INFO.NE.0 ) THEN
  219:          ITER = -2
  220:          GO TO 40
  221:       END IF
  222: *
  223: *     Compute the LU factorization of SA.
  224: *
  225:       CALL SGETRF( N, N, SWORK( PTSA ), N, IPIV, INFO )
  226: *
  227:       IF( INFO.NE.0 ) THEN
  228:          ITER = -3
  229:          GO TO 40
  230:       END IF
  231: *
  232: *     Solve the system SA*SX = SB.
  233: *
  234:       CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
  235:      +             SWORK( PTSX ), N, INFO )
  236: *
  237: *     Convert SX back to double precision
  238: *
  239:       CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
  240: *
  241: *     Compute R = B - AX (R is WORK).
  242: *
  243:       CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
  244: *
  245:       CALL DGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE, A,
  246:      +            LDA, X, LDX, ONE, WORK, N )
  247: *
  248: *     Check whether the NRHS normwise backward errors satisfy the
  249: *     stopping criterion. If yes, set ITER=0 and return.
  250: *
  251:       DO I = 1, NRHS
  252:          XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
  253:          RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
  254:          IF( RNRM.GT.XNRM*CTE )
  255:      +      GO TO 10
  256:       END DO
  257: *
  258: *     If we are here, the NRHS normwise backward errors satisfy the
  259: *     stopping criterion. We are good to exit.
  260: *
  261:       ITER = 0
  262:       RETURN
  263: *
  264:    10 CONTINUE
  265: *
  266:       DO 30 IITER = 1, ITERMAX
  267: *
  268: *        Convert R (in WORK) from double precision to single precision
  269: *        and store the result in SX.
  270: *
  271:          CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
  272: *
  273:          IF( INFO.NE.0 ) THEN
  274:             ITER = -2
  275:             GO TO 40
  276:          END IF
  277: *
  278: *        Solve the system SA*SX = SR.
  279: *
  280:          CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
  281:      +                SWORK( PTSX ), N, INFO )
  282: *
  283: *        Convert SX back to double precision and update the current
  284: *        iterate.
  285: *
  286:          CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
  287: *
  288:          DO I = 1, NRHS
  289:             CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
  290:          END DO
  291: *
  292: *        Compute R = B - AX (R is WORK).
  293: *
  294:          CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
  295: *
  296:          CALL DGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE,
  297:      +               A, LDA, X, LDX, ONE, WORK, N )
  298: *
  299: *        Check whether the NRHS normwise backward errors satisfy the
  300: *        stopping criterion. If yes, set ITER=IITER>0 and return.
  301: *
  302:          DO I = 1, NRHS
  303:             XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
  304:             RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
  305:             IF( RNRM.GT.XNRM*CTE )
  306:      +         GO TO 20
  307:          END DO
  308: *
  309: *        If we are here, the NRHS normwise backward errors satisfy the
  310: *        stopping criterion, we are good to exit.
  311: *
  312:          ITER = IITER
  313: *
  314:          RETURN
  315: *
  316:    20    CONTINUE
  317: *
  318:    30 CONTINUE
  319: *
  320: *     If we are at this place of the code, this is because we have
  321: *     performed ITER=ITERMAX iterations and never satisified the
  322: *     stopping criterion, set up the ITER flag accordingly and follow up
  323: *     on double precision routine.
  324: *
  325:       ITER = -ITERMAX - 1
  326: *
  327:    40 CONTINUE
  328: *
  329: *     Single-precision iterative refinement failed to converge to a
  330: *     satisfactory solution, so we resort to double precision.
  331: *
  332:       CALL DGETRF( N, N, A, LDA, IPIV, INFO )
  333: *
  334:       IF( INFO.NE.0 )
  335:      +   RETURN
  336: *
  337:       CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
  338:       CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, X, LDX,
  339:      +             INFO )
  340: *
  341:       RETURN
  342: *
  343: *     End of DSGESV.
  344: *
  345:       END