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