LAPACK: SRC/cgerfs.f Source File

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