LAPACK: SRC/zgbrfs.f Source File

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