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