LAPACK 3.3.1
Linear Algebra PACKage

dsgesv.f

Go to the documentation of this file.
00001       SUBROUTINE DSGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
00002      $                   SWORK, ITER, INFO )
00003 *
00004 *  -- LAPACK PROTOTYPE driver routine (version 3.3.1) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *  -- April 2011                                                      --
00008 *
00009 *     ..
00010 *     .. Scalar Arguments ..
00011       INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            IPIV( * )
00015       REAL               SWORK( * )
00016       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( N, * ),
00017      $                   X( LDX, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  DSGESV computes the solution to a real system of linear equations
00024 *     A * X = B,
00025 *  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
00026 *
00027 *  DSGESV first attempts to factorize the matrix in SINGLE PRECISION
00028 *  and use this factorization within an iterative refinement procedure
00029 *  to produce a solution with DOUBLE PRECISION normwise backward error
00030 *  quality (see below). If the approach fails the method switches to a
00031 *  DOUBLE PRECISION factorization and solve.
00032 *
00033 *  The iterative refinement is not going to be a winning strategy if
00034 *  the ratio SINGLE PRECISION performance over DOUBLE PRECISION
00035 *  performance is too small. A reasonable strategy should take the
00036 *  number of right-hand sides and the size of the matrix into account.
00037 *  This might be done with a call to ILAENV in the future. Up to now, we
00038 *  always try iterative refinement.
00039 *
00040 *  The iterative refinement process is stopped if
00041 *      ITER > ITERMAX
00042 *  or for all the RHS we have:
00043 *      RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
00044 *  where
00045 *      o ITER is the number of the current iteration in the iterative
00046 *        refinement process
00047 *      o RNRM is the infinity-norm of the residual
00048 *      o XNRM is the infinity-norm of the solution
00049 *      o ANRM is the infinity-operator-norm of the matrix A
00050 *      o EPS is the machine epsilon returned by DLAMCH('Epsilon')
00051 *  The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
00052 *  respectively.
00053 *
00054 *  Arguments
00055 *  =========
00056 *
00057 *  N       (input) INTEGER
00058 *          The number of linear equations, i.e., the order of the
00059 *          matrix A.  N >= 0.
00060 *
00061 *  NRHS    (input) INTEGER
00062 *          The number of right hand sides, i.e., the number of columns
00063 *          of the matrix B.  NRHS >= 0.
00064 *
00065 *  A       (input/output) DOUBLE PRECISION array,
00066 *          dimension (LDA,N)
00067 *          On entry, the N-by-N coefficient matrix A.
00068 *          On exit, if iterative refinement has been successfully used
00069 *          (INFO.EQ.0 and ITER.GE.0, see description below), then A is
00070 *          unchanged, if double precision factorization has been used
00071 *          (INFO.EQ.0 and ITER.LT.0, see description below), then the
00072 *          array A contains the factors L and U from the factorization
00073 *          A = P*L*U; the unit diagonal elements of L are not stored.
00074 *
00075 *  LDA     (input) INTEGER
00076 *          The leading dimension of the array A.  LDA >= max(1,N).
00077 *
00078 *  IPIV    (output) INTEGER array, dimension (N)
00079 *          The pivot indices that define the permutation matrix P;
00080 *          row i of the matrix was interchanged with row IPIV(i).
00081 *          Corresponds either to the single precision factorization
00082 *          (if INFO.EQ.0 and ITER.GE.0) or the double precision
00083 *          factorization (if INFO.EQ.0 and ITER.LT.0).
00084 *
00085 *  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
00086 *          The N-by-NRHS right hand side matrix B.
00087 *
00088 *  LDB     (input) INTEGER
00089 *          The leading dimension of the array B.  LDB >= max(1,N).
00090 *
00091 *  X       (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
00092 *          If INFO = 0, the N-by-NRHS solution matrix X.
00093 *
00094 *  LDX     (input) INTEGER
00095 *          The leading dimension of the array X.  LDX >= max(1,N).
00096 *
00097 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N,NRHS)
00098 *          This array is used to hold the residual vectors.
00099 *
00100 *  SWORK   (workspace) REAL array, dimension (N*(N+NRHS))
00101 *          This array is used to use the single precision matrix and the
00102 *          right-hand sides or solutions in single precision.
00103 *
00104 *  ITER    (output) INTEGER
00105 *          < 0: iterative refinement has failed, double precision
00106 *               factorization has been performed
00107 *               -1 : the routine fell back to full precision for
00108 *                    implementation- or machine-specific reasons
00109 *               -2 : narrowing the precision induced an overflow,
00110 *                    the routine fell back to full precision
00111 *               -3 : failure of SGETRF
00112 *               -31: stop the iterative refinement after the 30th
00113 *                    iterations
00114 *          > 0: iterative refinement has been sucessfully used.
00115 *               Returns the number of iterations
00116 *
00117 *  INFO    (output) INTEGER
00118 *          = 0:  successful exit
00119 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00120 *          > 0:  if INFO = i, U(i,i) computed in DOUBLE PRECISION is
00121 *                exactly zero.  The factorization has been completed,
00122 *                but the factor U is exactly singular, so the solution
00123 *                could not be computed.
00124 *
00125 *  =====================================================================
00126 *
00127 *     .. Parameters ..
00128       LOGICAL            DOITREF
00129       PARAMETER          ( DOITREF = .TRUE. )
00130 *
00131       INTEGER            ITERMAX
00132       PARAMETER          ( ITERMAX = 30 )
00133 *
00134       DOUBLE PRECISION   BWDMAX
00135       PARAMETER          ( BWDMAX = 1.0E+00 )
00136 *
00137       DOUBLE PRECISION   NEGONE, ONE
00138       PARAMETER          ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
00139 *
00140 *     .. Local Scalars ..
00141       INTEGER            I, IITER, PTSA, PTSX
00142       DOUBLE PRECISION   ANRM, CTE, EPS, RNRM, XNRM
00143 *
00144 *     .. External Subroutines ..
00145       EXTERNAL           DAXPY, DGEMM, DLACPY, DLAG2S, SLAG2D, SGETRF,
00146      $                   SGETRS, XERBLA
00147 *     ..
00148 *     .. External Functions ..
00149       INTEGER            IDAMAX
00150       DOUBLE PRECISION   DLAMCH, DLANGE
00151       EXTERNAL           IDAMAX, DLAMCH, DLANGE
00152 *     ..
00153 *     .. Intrinsic Functions ..
00154       INTRINSIC          ABS, DBLE, MAX, SQRT
00155 *     ..
00156 *     .. Executable Statements ..
00157 *
00158       INFO = 0
00159       ITER = 0
00160 *
00161 *     Test the input parameters.
00162 *
00163       IF( N.LT.0 ) THEN
00164          INFO = -1
00165       ELSE IF( NRHS.LT.0 ) THEN
00166          INFO = -2
00167       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00168          INFO = -4
00169       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00170          INFO = -7
00171       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00172          INFO = -9
00173       END IF
00174       IF( INFO.NE.0 ) THEN
00175          CALL XERBLA( 'DSGESV', -INFO )
00176          RETURN
00177       END IF
00178 *
00179 *     Quick return if (N.EQ.0).
00180 *
00181       IF( N.EQ.0 )
00182      $   RETURN
00183 *
00184 *     Skip single precision iterative refinement if a priori slower
00185 *     than double precision factorization.
00186 *
00187       IF( .NOT.DOITREF ) THEN
00188          ITER = -1
00189          GO TO 40
00190       END IF
00191 *
00192 *     Compute some constants.
00193 *
00194       ANRM = DLANGE( 'I', N, N, A, LDA, WORK )
00195       EPS = DLAMCH( 'Epsilon' )
00196       CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
00197 *
00198 *     Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
00199 *
00200       PTSA = 1
00201       PTSX = PTSA + N*N
00202 *
00203 *     Convert B from double precision to single precision and store the
00204 *     result in SX.
00205 *
00206       CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
00207 *
00208       IF( INFO.NE.0 ) THEN
00209          ITER = -2
00210          GO TO 40
00211       END IF
00212 *
00213 *     Convert A from double precision to single precision and store the
00214 *     result in SA.
00215 *
00216       CALL DLAG2S( N, N, A, LDA, SWORK( PTSA ), N, INFO )
00217 *
00218       IF( INFO.NE.0 ) THEN
00219          ITER = -2
00220          GO TO 40
00221       END IF
00222 *
00223 *     Compute the LU factorization of SA.
00224 *
00225       CALL SGETRF( N, N, SWORK( PTSA ), N, IPIV, INFO )
00226 *
00227       IF( INFO.NE.0 ) THEN
00228          ITER = -3
00229          GO TO 40
00230       END IF
00231 *
00232 *     Solve the system SA*SX = SB.
00233 *
00234       CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
00235      $             SWORK( PTSX ), N, INFO )
00236 *
00237 *     Convert SX back to double precision
00238 *
00239       CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
00240 *
00241 *     Compute R = B - AX (R is WORK).
00242 *
00243       CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
00244 *
00245       CALL DGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE, A,
00246      $            LDA, X, LDX, ONE, WORK, N )
00247 *
00248 *     Check whether the NRHS normwise backward errors satisfy the
00249 *     stopping criterion. If yes, set ITER=0 and return.
00250 *
00251       DO I = 1, NRHS
00252          XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
00253          RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
00254          IF( RNRM.GT.XNRM*CTE )
00255      $      GO TO 10
00256       END DO
00257 *
00258 *     If we are here, the NRHS normwise backward errors satisfy the
00259 *     stopping criterion. We are good to exit.
00260 *
00261       ITER = 0
00262       RETURN
00263 *
00264    10 CONTINUE
00265 *
00266       DO 30 IITER = 1, ITERMAX
00267 *
00268 *        Convert R (in WORK) from double precision to single precision
00269 *        and store the result in SX.
00270 *
00271          CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
00272 *
00273          IF( INFO.NE.0 ) THEN
00274             ITER = -2
00275             GO TO 40
00276          END IF
00277 *
00278 *        Solve the system SA*SX = SR.
00279 *
00280          CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
00281      $                SWORK( PTSX ), N, INFO )
00282 *
00283 *        Convert SX back to double precision and update the current
00284 *        iterate.
00285 *
00286          CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
00287 *
00288          DO I = 1, NRHS
00289             CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
00290          END DO
00291 *
00292 *        Compute R = B - AX (R is WORK).
00293 *
00294          CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
00295 *
00296          CALL DGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE,
00297      $               A, LDA, X, LDX, ONE, WORK, N )
00298 *
00299 *        Check whether the NRHS normwise backward errors satisfy the
00300 *        stopping criterion. If yes, set ITER=IITER>0 and return.
00301 *
00302          DO I = 1, NRHS
00303             XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
00304             RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
00305             IF( RNRM.GT.XNRM*CTE )
00306      $         GO TO 20
00307          END DO
00308 *
00309 *        If we are here, the NRHS normwise backward errors satisfy the
00310 *        stopping criterion, we are good to exit.
00311 *
00312          ITER = IITER
00313 *
00314          RETURN
00315 *
00316    20    CONTINUE
00317 *
00318    30 CONTINUE
00319 *
00320 *     If we are at this place of the code, this is because we have
00321 *     performed ITER=ITERMAX iterations and never satisified the
00322 *     stopping criterion, set up the ITER flag accordingly and follow up
00323 *     on double precision routine.
00324 *
00325       ITER = -ITERMAX - 1
00326 *
00327    40 CONTINUE
00328 *
00329 *     Single-precision iterative refinement failed to converge to a
00330 *     satisfactory solution, so we resort to double precision.
00331 *
00332       CALL DGETRF( N, N, A, LDA, IPIV, INFO )
00333 *
00334       IF( INFO.NE.0 )
00335      $   RETURN
00336 *
00337       CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
00338       CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, X, LDX,
00339      $             INFO )
00340 *
00341       RETURN
00342 *
00343 *     End of DSGESV.
00344 *
00345       END
 All Files Functions