LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cgerfs.f
Go to the documentation of this file.
00001 *> \brief \b CGERFS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CGERFS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgerfs.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgerfs.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgerfs.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
00022 *                          X, LDX, FERR, BERR, WORK, RWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          TRANS
00026 *       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IPIV( * )
00030 *       REAL               BERR( * ), FERR( * ), RWORK( * )
00031 *       COMPLEX            A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00032 *      $                   WORK( * ), X( LDX, * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> CGERFS improves the computed solution to a system of linear
00042 *> equations and provides error bounds and backward error estimates for
00043 *> the solution.
00044 *> \endverbatim
00045 *
00046 *  Arguments:
00047 *  ==========
00048 *
00049 *> \param[in] TRANS
00050 *> \verbatim
00051 *>          TRANS is CHARACTER*1
00052 *>          Specifies the form of the system of equations:
00053 *>          = 'N':  A * X = B     (No transpose)
00054 *>          = 'T':  A**T * X = B  (Transpose)
00055 *>          = 'C':  A**H * X = B  (Conjugate transpose)
00056 *> \endverbatim
00057 *>
00058 *> \param[in] N
00059 *> \verbatim
00060 *>          N is INTEGER
00061 *>          The order of the matrix A.  N >= 0.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] NRHS
00065 *> \verbatim
00066 *>          NRHS is INTEGER
00067 *>          The number of right hand sides, i.e., the number of columns
00068 *>          of the matrices B and X.  NRHS >= 0.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] A
00072 *> \verbatim
00073 *>          A is COMPLEX array, dimension (LDA,N)
00074 *>          The original N-by-N matrix A.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] LDA
00078 *> \verbatim
00079 *>          LDA is INTEGER
00080 *>          The leading dimension of the array A.  LDA >= max(1,N).
00081 *> \endverbatim
00082 *>
00083 *> \param[in] AF
00084 *> \verbatim
00085 *>          AF is COMPLEX array, dimension (LDAF,N)
00086 *>          The factors L and U from the factorization A = P*L*U
00087 *>          as computed by CGETRF.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] LDAF
00091 *> \verbatim
00092 *>          LDAF is INTEGER
00093 *>          The leading dimension of the array AF.  LDAF >= max(1,N).
00094 *> \endverbatim
00095 *>
00096 *> \param[in] IPIV
00097 *> \verbatim
00098 *>          IPIV is INTEGER array, dimension (N)
00099 *>          The pivot indices from CGETRF; for 1<=i<=N, row i of the
00100 *>          matrix was interchanged with row IPIV(i).
00101 *> \endverbatim
00102 *>
00103 *> \param[in] B
00104 *> \verbatim
00105 *>          B is COMPLEX array, dimension (LDB,NRHS)
00106 *>          The right hand side matrix B.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] LDB
00110 *> \verbatim
00111 *>          LDB is INTEGER
00112 *>          The leading dimension of the array B.  LDB >= max(1,N).
00113 *> \endverbatim
00114 *>
00115 *> \param[in,out] X
00116 *> \verbatim
00117 *>          X is COMPLEX array, dimension (LDX,NRHS)
00118 *>          On entry, the solution matrix X, as computed by CGETRS.
00119 *>          On exit, the improved solution matrix X.
00120 *> \endverbatim
00121 *>
00122 *> \param[in] LDX
00123 *> \verbatim
00124 *>          LDX is INTEGER
00125 *>          The leading dimension of the array X.  LDX >= max(1,N).
00126 *> \endverbatim
00127 *>
00128 *> \param[out] FERR
00129 *> \verbatim
00130 *>          FERR is REAL array, dimension (NRHS)
00131 *>          The estimated forward error bound for each solution vector
00132 *>          X(j) (the j-th column of the solution matrix X).
00133 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
00134 *>          is an estimated upper bound for the magnitude of the largest
00135 *>          element in (X(j) - XTRUE) divided by the magnitude of the
00136 *>          largest element in X(j).  The estimate is as reliable as
00137 *>          the estimate for RCOND, and is almost always a slight
00138 *>          overestimate of the true error.
00139 *> \endverbatim
00140 *>
00141 *> \param[out] BERR
00142 *> \verbatim
00143 *>          BERR is REAL array, dimension (NRHS)
00144 *>          The componentwise relative backward error of each solution
00145 *>          vector X(j) (i.e., the smallest relative change in
00146 *>          any element of A or B that makes X(j) an exact solution).
00147 *> \endverbatim
00148 *>
00149 *> \param[out] WORK
00150 *> \verbatim
00151 *>          WORK is COMPLEX array, dimension (2*N)
00152 *> \endverbatim
00153 *>
00154 *> \param[out] RWORK
00155 *> \verbatim
00156 *>          RWORK is REAL array, dimension (N)
00157 *> \endverbatim
00158 *>
00159 *> \param[out] INFO
00160 *> \verbatim
00161 *>          INFO is INTEGER
00162 *>          = 0:  successful exit
00163 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00164 *> \endverbatim
00165 *
00166 *> \par Internal Parameters:
00167 *  =========================
00168 *>
00169 *> \verbatim
00170 *>  ITMAX is the maximum number of steps of iterative refinement.
00171 *> \endverbatim
00172 *
00173 *  Authors:
00174 *  ========
00175 *
00176 *> \author Univ. of Tennessee 
00177 *> \author Univ. of California Berkeley 
00178 *> \author Univ. of Colorado Denver 
00179 *> \author NAG Ltd. 
00180 *
00181 *> \date November 2011
00182 *
00183 *> \ingroup complexGEcomputational
00184 *
00185 *  =====================================================================
00186       SUBROUTINE CGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
00187      $                   X, LDX, FERR, BERR, WORK, RWORK, INFO )
00188 *
00189 *  -- LAPACK computational routine (version 3.4.0) --
00190 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00191 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00192 *     November 2011
00193 *
00194 *     .. Scalar Arguments ..
00195       CHARACTER          TRANS
00196       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
00197 *     ..
00198 *     .. Array Arguments ..
00199       INTEGER            IPIV( * )
00200       REAL               BERR( * ), FERR( * ), RWORK( * )
00201       COMPLEX            A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00202      $                   WORK( * ), X( LDX, * )
00203 *     ..
00204 *
00205 *  =====================================================================
00206 *
00207 *     .. Parameters ..
00208       INTEGER            ITMAX
00209       PARAMETER          ( ITMAX = 5 )
00210       REAL               ZERO
00211       PARAMETER          ( ZERO = 0.0E+0 )
00212       COMPLEX            ONE
00213       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00214       REAL               TWO
00215       PARAMETER          ( TWO = 2.0E+0 )
00216       REAL               THREE
00217       PARAMETER          ( THREE = 3.0E+0 )
00218 *     ..
00219 *     .. Local Scalars ..
00220       LOGICAL            NOTRAN
00221       CHARACTER          TRANSN, TRANST
00222       INTEGER            COUNT, I, J, K, KASE, NZ
00223       REAL               EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
00224       COMPLEX            ZDUM
00225 *     ..
00226 *     .. Local Arrays ..
00227       INTEGER            ISAVE( 3 )
00228 *     ..
00229 *     .. External Functions ..
00230       LOGICAL            LSAME
00231       REAL               SLAMCH
00232       EXTERNAL           LSAME, SLAMCH
00233 *     ..
00234 *     .. External Subroutines ..
00235       EXTERNAL           CAXPY, CCOPY, CGEMV, CGETRS, CLACN2, XERBLA
00236 *     ..
00237 *     .. Intrinsic Functions ..
00238       INTRINSIC          ABS, AIMAG, MAX, REAL
00239 *     ..
00240 *     .. Statement Functions ..
00241       REAL               CABS1
00242 *     ..
00243 *     .. Statement Function definitions ..
00244       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00245 *     ..
00246 *     .. Executable Statements ..
00247 *
00248 *     Test the input parameters.
00249 *
00250       INFO = 0
00251       NOTRAN = LSAME( TRANS, 'N' )
00252       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00253      $    LSAME( TRANS, 'C' ) ) THEN
00254          INFO = -1
00255       ELSE IF( N.LT.0 ) THEN
00256          INFO = -2
00257       ELSE IF( NRHS.LT.0 ) THEN
00258          INFO = -3
00259       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00260          INFO = -5
00261       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
00262          INFO = -7
00263       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00264          INFO = -10
00265       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00266          INFO = -12
00267       END IF
00268       IF( INFO.NE.0 ) THEN
00269          CALL XERBLA( 'CGERFS', -INFO )
00270          RETURN
00271       END IF
00272 *
00273 *     Quick return if possible
00274 *
00275       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00276          DO 10 J = 1, NRHS
00277             FERR( J ) = ZERO
00278             BERR( J ) = ZERO
00279    10    CONTINUE
00280          RETURN
00281       END IF
00282 *
00283       IF( NOTRAN ) THEN
00284          TRANSN = 'N'
00285          TRANST = 'C'
00286       ELSE
00287          TRANSN = 'C'
00288          TRANST = 'N'
00289       END IF
00290 *
00291 *     NZ = maximum number of nonzero elements in each row of A, plus 1
00292 *
00293       NZ = N + 1
00294       EPS = SLAMCH( 'Epsilon' )
00295       SAFMIN = SLAMCH( 'Safe minimum' )
00296       SAFE1 = NZ*SAFMIN
00297       SAFE2 = SAFE1 / EPS
00298 *
00299 *     Do for each right hand side
00300 *
00301       DO 140 J = 1, NRHS
00302 *
00303          COUNT = 1
00304          LSTRES = THREE
00305    20    CONTINUE
00306 *
00307 *        Loop until stopping criterion is satisfied.
00308 *
00309 *        Compute residual R = B - op(A) * X,
00310 *        where op(A) = A, A**T, or A**H, depending on TRANS.
00311 *
00312          CALL CCOPY( N, B( 1, J ), 1, WORK, 1 )
00313          CALL CGEMV( TRANS, N, N, -ONE, A, LDA, X( 1, J ), 1, ONE, WORK,
00314      $               1 )
00315 *
00316 *        Compute componentwise relative backward error from formula
00317 *
00318 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
00319 *
00320 *        where abs(Z) is the componentwise absolute value of the matrix
00321 *        or vector Z.  If the i-th component of the denominator is less
00322 *        than SAFE2, then SAFE1 is added to the i-th components of the
00323 *        numerator and denominator before dividing.
00324 *
00325          DO 30 I = 1, N
00326             RWORK( I ) = CABS1( B( I, J ) )
00327    30    CONTINUE
00328 *
00329 *        Compute abs(op(A))*abs(X) + abs(B).
00330 *
00331          IF( NOTRAN ) THEN
00332             DO 50 K = 1, N
00333                XK = CABS1( X( K, J ) )
00334                DO 40 I = 1, N
00335                   RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00336    40          CONTINUE
00337    50       CONTINUE
00338          ELSE
00339             DO 70 K = 1, N
00340                S = ZERO
00341                DO 60 I = 1, N
00342                   S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00343    60          CONTINUE
00344                RWORK( K ) = RWORK( K ) + S
00345    70       CONTINUE
00346          END IF
00347          S = ZERO
00348          DO 80 I = 1, N
00349             IF( RWORK( I ).GT.SAFE2 ) THEN
00350                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
00351             ELSE
00352                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
00353      $             ( RWORK( I )+SAFE1 ) )
00354             END IF
00355    80    CONTINUE
00356          BERR( J ) = S
00357 *
00358 *        Test stopping criterion. Continue iterating if
00359 *           1) The residual BERR(J) is larger than machine epsilon, and
00360 *           2) BERR(J) decreased by at least a factor of 2 during the
00361 *              last iteration, and
00362 *           3) At most ITMAX iterations tried.
00363 *
00364          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
00365      $       COUNT.LE.ITMAX ) THEN
00366 *
00367 *           Update solution and try again.
00368 *
00369             CALL CGETRS( TRANS, N, 1, AF, LDAF, IPIV, WORK, N, INFO )
00370             CALL CAXPY( N, ONE, WORK, 1, X( 1, J ), 1 )
00371             LSTRES = BERR( J )
00372             COUNT = COUNT + 1
00373             GO TO 20
00374          END IF
00375 *
00376 *        Bound error from formula
00377 *
00378 *        norm(X - XTRUE) / norm(X) .le. FERR =
00379 *        norm( abs(inv(op(A)))*
00380 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
00381 *
00382 *        where
00383 *          norm(Z) is the magnitude of the largest component of Z
00384 *          inv(op(A)) is the inverse of op(A)
00385 *          abs(Z) is the componentwise absolute value of the matrix or
00386 *             vector Z
00387 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00388 *          EPS is machine epsilon
00389 *
00390 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
00391 *        is incremented by SAFE1 if the i-th component of
00392 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
00393 *
00394 *        Use CLACN2 to estimate the infinity-norm of the matrix
00395 *           inv(op(A)) * diag(W),
00396 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
00397 *
00398          DO 90 I = 1, N
00399             IF( RWORK( I ).GT.SAFE2 ) THEN
00400                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
00401             ELSE
00402                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
00403      $                      SAFE1
00404             END IF
00405    90    CONTINUE
00406 *
00407          KASE = 0
00408   100    CONTINUE
00409          CALL CLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
00410          IF( KASE.NE.0 ) THEN
00411             IF( KASE.EQ.1 ) THEN
00412 *
00413 *              Multiply by diag(W)*inv(op(A)**H).
00414 *
00415                CALL CGETRS( TRANST, N, 1, AF, LDAF, IPIV, WORK, N,
00416      $                      INFO )
00417                DO 110 I = 1, N
00418                   WORK( I ) = RWORK( I )*WORK( I )
00419   110          CONTINUE
00420             ELSE
00421 *
00422 *              Multiply by inv(op(A))*diag(W).
00423 *
00424                DO 120 I = 1, N
00425                   WORK( I ) = RWORK( I )*WORK( I )
00426   120          CONTINUE
00427                CALL CGETRS( TRANSN, N, 1, AF, LDAF, IPIV, WORK, N,
00428      $                      INFO )
00429             END IF
00430             GO TO 100
00431          END IF
00432 *
00433 *        Normalize error.
00434 *
00435          LSTRES = ZERO
00436          DO 130 I = 1, N
00437             LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
00438   130    CONTINUE
00439          IF( LSTRES.NE.ZERO )
00440      $      FERR( J ) = FERR( J ) / LSTRES
00441 *
00442   140 CONTINUE
00443 *
00444       RETURN
00445 *
00446 *     End of CGERFS
00447 *
00448       END
 All Files Functions