LAPACK: SRC/dtrrfs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
00002      $                   LDX, FERR, BERR, WORK, IWORK, 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 DLACN2 in place of DLACON, 5 Feb 03, SJH.
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          DIAG, TRANS, UPLO
00013       INTEGER            INFO, LDA, LDB, LDX, N, NRHS
00014 *     ..
00015 *     .. Array Arguments ..
00016       INTEGER            IWORK( * )
00017       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
00018      $                   WORK( * ), X( LDX, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  DTRRFS provides error bounds and backward error estimates for the
00025 *  solution to a system of linear equations with a triangular
00026 *  coefficient matrix.
00027 *
00028 *  The solution matrix X must be computed by DTRTRS or some other
00029 *  means before entering this routine.  DTRRFS does not do iterative
00030 *  refinement because doing so cannot improve the backward error.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  UPLO    (input) CHARACTER*1
00036 *          = 'U':  A is upper triangular;
00037 *          = 'L':  A is lower triangular.
00038 *
00039 *  TRANS   (input) CHARACTER*1
00040 *          Specifies the form of the system of equations:
00041 *          = 'N':  A * X = B  (No transpose)
00042 *          = 'T':  A**T * X = B  (Transpose)
00043 *          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
00044 *
00045 *  DIAG    (input) CHARACTER*1
00046 *          = 'N':  A is non-unit triangular;
00047 *          = 'U':  A is unit triangular.
00048 *
00049 *  N       (input) INTEGER
00050 *          The order of the matrix A.  N >= 0.
00051 *
00052 *  NRHS    (input) INTEGER
00053 *          The number of right hand sides, i.e., the number of columns
00054 *          of the matrices B and X.  NRHS >= 0.
00055 *
00056 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00057 *          The triangular matrix A.  If UPLO = 'U', the leading N-by-N
00058 *          upper triangular part of the array A contains the upper
00059 *          triangular matrix, and the strictly lower triangular part of
00060 *          A is not referenced.  If UPLO = 'L', the leading N-by-N lower
00061 *          triangular part of the array A contains the lower triangular
00062 *          matrix, and the strictly upper triangular part of A is not
00063 *          referenced.  If DIAG = 'U', the diagonal elements of A are
00064 *          also not referenced and are assumed to be 1.
00065 *
00066 *  LDA     (input) INTEGER
00067 *          The leading dimension of the array A.  LDA >= max(1,N).
00068 *
00069 *  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
00070 *          The right hand side matrix B.
00071 *
00072 *  LDB     (input) INTEGER
00073 *          The leading dimension of the array B.  LDB >= max(1,N).
00074 *
00075 *  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
00076 *          The solution matrix X.
00077 *
00078 *  LDX     (input) INTEGER
00079 *          The leading dimension of the array X.  LDX >= max(1,N).
00080 *
00081 *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
00082 *          The estimated forward error bound for each solution vector
00083 *          X(j) (the j-th column of the solution matrix X).
00084 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00085 *          is an estimated upper bound for the magnitude of the largest
00086 *          element in (X(j) - XTRUE) divided by the magnitude of the
00087 *          largest element in X(j).  The estimate is as reliable as
00088 *          the estimate for RCOND, and is almost always a slight
00089 *          overestimate of the true error.
00090 *
00091 *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
00092 *          The componentwise relative backward error of each solution
00093 *          vector X(j) (i.e., the smallest relative change in
00094 *          any element of A or B that makes X(j) an exact solution).
00095 *
00096 *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
00097 *
00098 *  IWORK   (workspace) INTEGER array, dimension (N)
00099 *
00100 *  INFO    (output) INTEGER
00101 *          = 0:  successful exit
00102 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00103 *
00104 *  =====================================================================
00105 *
00106 *     .. Parameters ..
00107       DOUBLE PRECISION   ZERO
00108       PARAMETER          ( ZERO = 0.0D+0 )
00109       DOUBLE PRECISION   ONE
00110       PARAMETER          ( ONE = 1.0D+0 )
00111 *     ..
00112 *     .. Local Scalars ..
00113       LOGICAL            NOTRAN, NOUNIT, UPPER
00114       CHARACTER          TRANST
00115       INTEGER            I, J, K, KASE, NZ
00116       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
00117 *     ..
00118 *     .. Local Arrays ..
00119       INTEGER            ISAVE( 3 )
00120 *     ..
00121 *     .. External Subroutines ..
00122       EXTERNAL           DAXPY, DCOPY, DLACN2, DTRMV, DTRSV, XERBLA
00123 *     ..
00124 *     .. Intrinsic Functions ..
00125       INTRINSIC          ABS, MAX
00126 *     ..
00127 *     .. External Functions ..
00128       LOGICAL            LSAME
00129       DOUBLE PRECISION   DLAMCH
00130       EXTERNAL           LSAME, DLAMCH
00131 *     ..
00132 *     .. Executable Statements ..
00133 *
00134 *     Test the input parameters.
00135 *
00136       INFO = 0
00137       UPPER = LSAME( UPLO, 'U' )
00138       NOTRAN = LSAME( TRANS, 'N' )
00139       NOUNIT = LSAME( DIAG, 'N' )
00140 *
00141       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00142          INFO = -1
00143       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00144      $         LSAME( TRANS, 'C' ) ) THEN
00145          INFO = -2
00146       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00147          INFO = -3
00148       ELSE IF( N.LT.0 ) THEN
00149          INFO = -4
00150       ELSE IF( NRHS.LT.0 ) THEN
00151          INFO = -5
00152       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00153          INFO = -7
00154       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00155          INFO = -9
00156       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00157          INFO = -11
00158       END IF
00159       IF( INFO.NE.0 ) THEN
00160          CALL XERBLA( 'DTRRFS', -INFO )
00161          RETURN
00162       END IF
00163 *
00164 *     Quick return if possible
00165 *
00166       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00167          DO 10 J = 1, NRHS
00168             FERR( J ) = ZERO
00169             BERR( J ) = ZERO
00170    10    CONTINUE
00171          RETURN
00172       END IF
00173 *
00174       IF( NOTRAN ) THEN
00175          TRANST = 'T'
00176       ELSE
00177          TRANST = 'N'
00178       END IF
00179 *
00180 *     NZ = maximum number of nonzero elements in each row of A, plus 1
00181 *
00182       NZ = N + 1
00183       EPS = DLAMCH( 'Epsilon' )
00184       SAFMIN = DLAMCH( 'Safe minimum' )
00185       SAFE1 = NZ*SAFMIN
00186       SAFE2 = SAFE1 / EPS
00187 *
00188 *     Do for each right hand side
00189 *
00190       DO 250 J = 1, NRHS
00191 *
00192 *        Compute residual R = B - op(A) * X,
00193 *        where op(A) = A or A', depending on TRANS.
00194 *
00195          CALL DCOPY( N, X( 1, J ), 1, WORK( N+1 ), 1 )
00196          CALL DTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK( N+1 ), 1 )
00197          CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK( N+1 ), 1 )
00198 *
00199 *        Compute componentwise relative backward error from formula
00200 *
00201 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
00202 *
00203 *        where abs(Z) is the componentwise absolute value of the matrix
00204 *        or vector Z.  If the i-th component of the denominator is less
00205 *        than SAFE2, then SAFE1 is added to the i-th components of the
00206 *        numerator and denominator before dividing.
00207 *
00208          DO 20 I = 1, N
00209             WORK( I ) = ABS( B( I, J ) )
00210    20    CONTINUE
00211 *
00212          IF( NOTRAN ) THEN
00213 *
00214 *           Compute abs(A)*abs(X) + abs(B).
00215 *
00216             IF( UPPER ) THEN
00217                IF( NOUNIT ) THEN
00218                   DO 40 K = 1, N
00219                      XK = ABS( X( K, J ) )
00220                      DO 30 I = 1, K
00221                         WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
00222    30                CONTINUE
00223    40             CONTINUE
00224                ELSE
00225                   DO 60 K = 1, N
00226                      XK = ABS( X( K, J ) )
00227                      DO 50 I = 1, K - 1
00228                         WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
00229    50                CONTINUE
00230                      WORK( K ) = WORK( K ) + XK
00231    60             CONTINUE
00232                END IF
00233             ELSE
00234                IF( NOUNIT ) THEN
00235                   DO 80 K = 1, N
00236                      XK = ABS( X( K, J ) )
00237                      DO 70 I = K, N
00238                         WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
00239    70                CONTINUE
00240    80             CONTINUE
00241                ELSE
00242                   DO 100 K = 1, N
00243                      XK = ABS( X( K, J ) )
00244                      DO 90 I = K + 1, N
00245                         WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
00246    90                CONTINUE
00247                      WORK( K ) = WORK( K ) + XK
00248   100             CONTINUE
00249                END IF
00250             END IF
00251          ELSE
00252 *
00253 *           Compute abs(A')*abs(X) + abs(B).
00254 *
00255             IF( UPPER ) THEN
00256                IF( NOUNIT ) THEN
00257                   DO 120 K = 1, N
00258                      S = ZERO
00259                      DO 110 I = 1, K
00260                         S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
00261   110                CONTINUE
00262                      WORK( K ) = WORK( K ) + S
00263   120             CONTINUE
00264                ELSE
00265                   DO 140 K = 1, N
00266                      S = ABS( X( K, J ) )
00267                      DO 130 I = 1, K - 1
00268                         S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
00269   130                CONTINUE
00270                      WORK( K ) = WORK( K ) + S
00271   140             CONTINUE
00272                END IF
00273             ELSE
00274                IF( NOUNIT ) THEN
00275                   DO 160 K = 1, N
00276                      S = ZERO
00277                      DO 150 I = K, N
00278                         S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
00279   150                CONTINUE
00280                      WORK( K ) = WORK( K ) + S
00281   160             CONTINUE
00282                ELSE
00283                   DO 180 K = 1, N
00284                      S = ABS( X( K, J ) )
00285                      DO 170 I = K + 1, N
00286                         S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
00287   170                CONTINUE
00288                      WORK( K ) = WORK( K ) + S
00289   180             CONTINUE
00290                END IF
00291             END IF
00292          END IF
00293          S = ZERO
00294          DO 190 I = 1, N
00295             IF( WORK( I ).GT.SAFE2 ) THEN
00296                S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
00297             ELSE
00298                S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
00299      $             ( WORK( I )+SAFE1 ) )
00300             END IF
00301   190    CONTINUE
00302          BERR( J ) = S
00303 *
00304 *        Bound error from formula
00305 *
00306 *        norm(X - XTRUE) / norm(X) .le. FERR =
00307 *        norm( abs(inv(op(A)))*
00308 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
00309 *
00310 *        where
00311 *          norm(Z) is the magnitude of the largest component of Z
00312 *          inv(op(A)) is the inverse of op(A)
00313 *          abs(Z) is the componentwise absolute value of the matrix or
00314 *             vector Z
00315 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00316 *          EPS is machine epsilon
00317 *
00318 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
00319 *        is incremented by SAFE1 if the i-th component of
00320 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
00321 *
00322 *        Use DLACN2 to estimate the infinity-norm of the matrix
00323 *           inv(op(A)) * diag(W),
00324 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
00325 *
00326          DO 200 I = 1, N
00327             IF( WORK( I ).GT.SAFE2 ) THEN
00328                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
00329             ELSE
00330                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
00331             END IF
00332   200    CONTINUE
00333 *
00334          KASE = 0
00335   210    CONTINUE
00336          CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
00337      $                KASE, ISAVE )
00338          IF( KASE.NE.0 ) THEN
00339             IF( KASE.EQ.1 ) THEN
00340 *
00341 *              Multiply by diag(W)*inv(op(A)').
00342 *
00343                CALL DTRSV( UPLO, TRANST, DIAG, N, A, LDA, WORK( N+1 ),
00344      $                     1 )
00345                DO 220 I = 1, N
00346                   WORK( N+I ) = WORK( I )*WORK( N+I )
00347   220          CONTINUE
00348             ELSE
00349 *
00350 *              Multiply by inv(op(A))*diag(W).
00351 *
00352                DO 230 I = 1, N
00353                   WORK( N+I ) = WORK( I )*WORK( N+I )
00354   230          CONTINUE
00355                CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, WORK( N+1 ),
00356      $                     1 )
00357             END IF
00358             GO TO 210
00359          END IF
00360 *
00361 *        Normalize error.
00362 *
00363          LSTRES = ZERO
00364          DO 240 I = 1, N
00365             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
00366   240    CONTINUE
00367          IF( LSTRES.NE.ZERO )
00368      $      FERR( J ) = FERR( J ) / LSTRES
00369 *
00370   250 CONTINUE
00371 *
00372       RETURN
00373 *
00374 *     End of DTRRFS
00375 *
00376       END