LAPACK: SRC/ctrrfs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
00002      $                   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          DIAG, TRANS, UPLO
00013       INTEGER            INFO, LDA, LDB, LDX, N, NRHS
00014 *     ..
00015 *     .. Array Arguments ..
00016       REAL               BERR( * ), FERR( * ), RWORK( * )
00017       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * ),
00018      $                   X( LDX, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  CTRRFS 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 CTRTRS or some other
00029 *  means before entering this routine.  CTRRFS 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)
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) COMPLEX 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) COMPLEX 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) COMPLEX 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) REAL 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) REAL 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) COMPLEX array, dimension (2*N)
00097 *
00098 *  RWORK   (workspace) REAL 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       REAL               ZERO
00108       PARAMETER          ( ZERO = 0.0E+0 )
00109       COMPLEX            ONE
00110       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00111 *     ..
00112 *     .. Local Scalars ..
00113       LOGICAL            NOTRAN, NOUNIT, UPPER
00114       CHARACTER          TRANSN, TRANST
00115       INTEGER            I, J, K, KASE, NZ
00116       REAL               EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
00117       COMPLEX            ZDUM
00118 *     ..
00119 *     .. Local Arrays ..
00120       INTEGER            ISAVE( 3 )
00121 *     ..
00122 *     .. External Subroutines ..
00123       EXTERNAL           CAXPY, CCOPY, CLACN2, CTRMV, CTRSV, XERBLA
00124 *     ..
00125 *     .. Intrinsic Functions ..
00126       INTRINSIC          ABS, AIMAG, MAX, REAL
00127 *     ..
00128 *     .. External Functions ..
00129       LOGICAL            LSAME
00130       REAL               SLAMCH
00131       EXTERNAL           LSAME, SLAMCH
00132 *     ..
00133 *     .. Statement Functions ..
00134       REAL               CABS1
00135 *     ..
00136 *     .. Statement Function definitions ..
00137       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00138 *     ..
00139 *     .. Executable Statements ..
00140 *
00141 *     Test the input parameters.
00142 *
00143       INFO = 0
00144       UPPER = LSAME( UPLO, 'U' )
00145       NOTRAN = LSAME( TRANS, 'N' )
00146       NOUNIT = LSAME( DIAG, 'N' )
00147 *
00148       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00149          INFO = -1
00150       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00151      $         LSAME( TRANS, 'C' ) ) THEN
00152          INFO = -2
00153       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00154          INFO = -3
00155       ELSE IF( N.LT.0 ) THEN
00156          INFO = -4
00157       ELSE IF( NRHS.LT.0 ) THEN
00158          INFO = -5
00159       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00160          INFO = -7
00161       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00162          INFO = -9
00163       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00164          INFO = -11
00165       END IF
00166       IF( INFO.NE.0 ) THEN
00167          CALL XERBLA( 'CTRRFS', -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 250 J = 1, NRHS
00200 *
00201 *        Compute residual R = B - op(A) * X,
00202 *        where op(A) = A, A**T, or A**H, depending on TRANS.
00203 *
00204          CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
00205          CALL CTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
00206          CALL CAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
00207 *
00208 *        Compute componentwise relative backward error from formula
00209 *
00210 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
00211 *
00212 *        where abs(Z) is the componentwise absolute value of the matrix
00213 *        or vector Z.  If the i-th component of the denominator is less
00214 *        than SAFE2, then SAFE1 is added to the i-th components of the
00215 *        numerator and denominator before dividing.
00216 *
00217          DO 20 I = 1, N
00218             RWORK( I ) = CABS1( B( I, J ) )
00219    20    CONTINUE
00220 *
00221          IF( NOTRAN ) THEN
00222 *
00223 *           Compute abs(A)*abs(X) + abs(B).
00224 *
00225             IF( UPPER ) THEN
00226                IF( NOUNIT ) THEN
00227                   DO 40 K = 1, N
00228                      XK = CABS1( X( K, J ) )
00229                      DO 30 I = 1, K
00230                         RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00231    30                CONTINUE
00232    40             CONTINUE
00233                ELSE
00234                   DO 60 K = 1, N
00235                      XK = CABS1( X( K, J ) )
00236                      DO 50 I = 1, K - 1
00237                         RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00238    50                CONTINUE
00239                      RWORK( K ) = RWORK( K ) + XK
00240    60             CONTINUE
00241                END IF
00242             ELSE
00243                IF( NOUNIT ) THEN
00244                   DO 80 K = 1, N
00245                      XK = CABS1( X( K, J ) )
00246                      DO 70 I = K, N
00247                         RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00248    70                CONTINUE
00249    80             CONTINUE
00250                ELSE
00251                   DO 100 K = 1, N
00252                      XK = CABS1( X( K, J ) )
00253                      DO 90 I = K + 1, N
00254                         RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00255    90                CONTINUE
00256                      RWORK( K ) = RWORK( K ) + XK
00257   100             CONTINUE
00258                END IF
00259             END IF
00260          ELSE
00261 *
00262 *           Compute abs(A**H)*abs(X) + abs(B).
00263 *
00264             IF( UPPER ) THEN
00265                IF( NOUNIT ) THEN
00266                   DO 120 K = 1, N
00267                      S = ZERO
00268                      DO 110 I = 1, K
00269                         S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00270   110                CONTINUE
00271                      RWORK( K ) = RWORK( K ) + S
00272   120             CONTINUE
00273                ELSE
00274                   DO 140 K = 1, N
00275                      S = CABS1( X( K, J ) )
00276                      DO 130 I = 1, K - 1
00277                         S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00278   130                CONTINUE
00279                      RWORK( K ) = RWORK( K ) + S
00280   140             CONTINUE
00281                END IF
00282             ELSE
00283                IF( NOUNIT ) THEN
00284                   DO 160 K = 1, N
00285                      S = ZERO
00286                      DO 150 I = K, N
00287                         S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00288   150                CONTINUE
00289                      RWORK( K ) = RWORK( K ) + S
00290   160             CONTINUE
00291                ELSE
00292                   DO 180 K = 1, N
00293                      S = CABS1( X( K, J ) )
00294                      DO 170 I = K + 1, N
00295                         S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00296   170                CONTINUE
00297                      RWORK( K ) = RWORK( K ) + S
00298   180             CONTINUE
00299                END IF
00300             END IF
00301          END IF
00302          S = ZERO
00303          DO 190 I = 1, N
00304             IF( RWORK( I ).GT.SAFE2 ) THEN
00305                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
00306             ELSE
00307                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
00308      $             ( RWORK( I )+SAFE1 ) )
00309             END IF
00310   190    CONTINUE
00311          BERR( J ) = S
00312 *
00313 *        Bound error from formula
00314 *
00315 *        norm(X - XTRUE) / norm(X) .le. FERR =
00316 *        norm( abs(inv(op(A)))*
00317 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
00318 *
00319 *        where
00320 *          norm(Z) is the magnitude of the largest component of Z
00321 *          inv(op(A)) is the inverse of op(A)
00322 *          abs(Z) is the componentwise absolute value of the matrix or
00323 *             vector Z
00324 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00325 *          EPS is machine epsilon
00326 *
00327 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
00328 *        is incremented by SAFE1 if the i-th component of
00329 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
00330 *
00331 *        Use CLACN2 to estimate the infinity-norm of the matrix
00332 *           inv(op(A)) * diag(W),
00333 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
00334 *
00335          DO 200 I = 1, N
00336             IF( RWORK( I ).GT.SAFE2 ) THEN
00337                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
00338             ELSE
00339                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
00340      $                      SAFE1
00341             END IF
00342   200    CONTINUE
00343 *
00344          KASE = 0
00345   210    CONTINUE
00346          CALL CLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
00347          IF( KASE.NE.0 ) THEN
00348             IF( KASE.EQ.1 ) THEN
00349 *
00350 *              Multiply by diag(W)*inv(op(A)**H).
00351 *
00352                CALL CTRSV( UPLO, TRANST, DIAG, N, A, LDA, WORK, 1 )
00353                DO 220 I = 1, N
00354                   WORK( I ) = RWORK( I )*WORK( I )
00355   220          CONTINUE
00356             ELSE
00357 *
00358 *              Multiply by inv(op(A))*diag(W).
00359 *
00360                DO 230 I = 1, N
00361                   WORK( I ) = RWORK( I )*WORK( I )
00362   230          CONTINUE
00363                CALL CTRSV( UPLO, TRANSN, DIAG, N, A, LDA, WORK, 1 )
00364             END IF
00365             GO TO 210
00366          END IF
00367 *
00368 *        Normalize error.
00369 *
00370          LSTRES = ZERO
00371          DO 240 I = 1, N
00372             LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
00373   240    CONTINUE
00374          IF( LSTRES.NE.ZERO )
00375      $      FERR( J ) = FERR( J ) / LSTRES
00376 *
00377   250 CONTINUE
00378 *
00379       RETURN
00380 *
00381 *     End of CTRRFS
00382 *
00383       END