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