LAPACK 3.3.1
Linear Algebra PACKage

zcgesv.f

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