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