LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ztrrfs.f
Go to the documentation of this file.
00001 *> \brief \b ZTRRFS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZTRRFS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrrfs.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrrfs.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrrfs.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
00022 *                          LDX, FERR, BERR, WORK, RWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          DIAG, TRANS, UPLO
00026 *       INTEGER            INFO, LDA, LDB, LDX, N, NRHS
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
00030 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * ),
00031 *      $                   X( LDX, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> ZTRRFS 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 ZTRTRS or some other
00045 *> means before entering this routine.  ZTRRFS 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)
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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 COMPLEX*16 array, dimension (2*N)
00155 *> \endverbatim
00156 *>
00157 *> \param[out] RWORK
00158 *> \verbatim
00159 *>          RWORK is DOUBLE PRECISION 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 complex16OTHERcomputational
00180 *
00181 *  =====================================================================
00182       SUBROUTINE ZTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
00183      $                   LDX, FERR, BERR, WORK, RWORK, 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       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
00196       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * ),
00197      $                   X( LDX, * )
00198 *     ..
00199 *
00200 *  =====================================================================
00201 *
00202 *     .. Parameters ..
00203       DOUBLE PRECISION   ZERO
00204       PARAMETER          ( ZERO = 0.0D+0 )
00205       COMPLEX*16         ONE
00206       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
00207 *     ..
00208 *     .. Local Scalars ..
00209       LOGICAL            NOTRAN, NOUNIT, UPPER
00210       CHARACTER          TRANSN, TRANST
00211       INTEGER            I, J, K, KASE, NZ
00212       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
00213       COMPLEX*16         ZDUM
00214 *     ..
00215 *     .. Local Arrays ..
00216       INTEGER            ISAVE( 3 )
00217 *     ..
00218 *     .. External Subroutines ..
00219       EXTERNAL           XERBLA, ZAXPY, ZCOPY, ZLACN2, ZTRMV, ZTRSV
00220 *     ..
00221 *     .. Intrinsic Functions ..
00222       INTRINSIC          ABS, DBLE, DIMAG, MAX
00223 *     ..
00224 *     .. External Functions ..
00225       LOGICAL            LSAME
00226       DOUBLE PRECISION   DLAMCH
00227       EXTERNAL           LSAME, DLAMCH
00228 *     ..
00229 *     .. Statement Functions ..
00230       DOUBLE PRECISION   CABS1
00231 *     ..
00232 *     .. Statement Function definitions ..
00233       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
00234 *     ..
00235 *     .. Executable Statements ..
00236 *
00237 *     Test the input parameters.
00238 *
00239       INFO = 0
00240       UPPER = LSAME( UPLO, 'U' )
00241       NOTRAN = LSAME( TRANS, 'N' )
00242       NOUNIT = LSAME( DIAG, 'N' )
00243 *
00244       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00245          INFO = -1
00246       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00247      $         LSAME( TRANS, 'C' ) ) THEN
00248          INFO = -2
00249       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00250          INFO = -3
00251       ELSE IF( N.LT.0 ) THEN
00252          INFO = -4
00253       ELSE IF( NRHS.LT.0 ) THEN
00254          INFO = -5
00255       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00256          INFO = -7
00257       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00258          INFO = -9
00259       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00260          INFO = -11
00261       END IF
00262       IF( INFO.NE.0 ) THEN
00263          CALL XERBLA( 'ZTRRFS', -INFO )
00264          RETURN
00265       END IF
00266 *
00267 *     Quick return if possible
00268 *
00269       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00270          DO 10 J = 1, NRHS
00271             FERR( J ) = ZERO
00272             BERR( J ) = ZERO
00273    10    CONTINUE
00274          RETURN
00275       END IF
00276 *
00277       IF( NOTRAN ) THEN
00278          TRANSN = 'N'
00279          TRANST = 'C'
00280       ELSE
00281          TRANSN = 'C'
00282          TRANST = 'N'
00283       END IF
00284 *
00285 *     NZ = maximum number of nonzero elements in each row of A, plus 1
00286 *
00287       NZ = N + 1
00288       EPS = DLAMCH( 'Epsilon' )
00289       SAFMIN = DLAMCH( 'Safe minimum' )
00290       SAFE1 = NZ*SAFMIN
00291       SAFE2 = SAFE1 / EPS
00292 *
00293 *     Do for each right hand side
00294 *
00295       DO 250 J = 1, NRHS
00296 *
00297 *        Compute residual R = B - op(A) * X,
00298 *        where op(A) = A, A**T, or A**H, depending on TRANS.
00299 *
00300          CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 )
00301          CALL ZTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
00302          CALL ZAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
00303 *
00304 *        Compute componentwise relative backward error from formula
00305 *
00306 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
00307 *
00308 *        where abs(Z) is the componentwise absolute value of the matrix
00309 *        or vector Z.  If the i-th component of the denominator is less
00310 *        than SAFE2, then SAFE1 is added to the i-th components of the
00311 *        numerator and denominator before dividing.
00312 *
00313          DO 20 I = 1, N
00314             RWORK( I ) = CABS1( B( I, J ) )
00315    20    CONTINUE
00316 *
00317          IF( NOTRAN ) THEN
00318 *
00319 *           Compute abs(A)*abs(X) + abs(B).
00320 *
00321             IF( UPPER ) THEN
00322                IF( NOUNIT ) THEN
00323                   DO 40 K = 1, N
00324                      XK = CABS1( X( K, J ) )
00325                      DO 30 I = 1, K
00326                         RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00327    30                CONTINUE
00328    40             CONTINUE
00329                ELSE
00330                   DO 60 K = 1, N
00331                      XK = CABS1( X( K, J ) )
00332                      DO 50 I = 1, K - 1
00333                         RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00334    50                CONTINUE
00335                      RWORK( K ) = RWORK( K ) + XK
00336    60             CONTINUE
00337                END IF
00338             ELSE
00339                IF( NOUNIT ) THEN
00340                   DO 80 K = 1, N
00341                      XK = CABS1( X( K, J ) )
00342                      DO 70 I = K, N
00343                         RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00344    70                CONTINUE
00345    80             CONTINUE
00346                ELSE
00347                   DO 100 K = 1, N
00348                      XK = CABS1( X( K, J ) )
00349                      DO 90 I = K + 1, N
00350                         RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
00351    90                CONTINUE
00352                      RWORK( K ) = RWORK( K ) + XK
00353   100             CONTINUE
00354                END IF
00355             END IF
00356          ELSE
00357 *
00358 *           Compute abs(A**H)*abs(X) + abs(B).
00359 *
00360             IF( UPPER ) THEN
00361                IF( NOUNIT ) THEN
00362                   DO 120 K = 1, N
00363                      S = ZERO
00364                      DO 110 I = 1, K
00365                         S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00366   110                CONTINUE
00367                      RWORK( K ) = RWORK( K ) + S
00368   120             CONTINUE
00369                ELSE
00370                   DO 140 K = 1, N
00371                      S = CABS1( X( K, J ) )
00372                      DO 130 I = 1, K - 1
00373                         S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00374   130                CONTINUE
00375                      RWORK( K ) = RWORK( K ) + S
00376   140             CONTINUE
00377                END IF
00378             ELSE
00379                IF( NOUNIT ) THEN
00380                   DO 160 K = 1, N
00381                      S = ZERO
00382                      DO 150 I = K, N
00383                         S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00384   150                CONTINUE
00385                      RWORK( K ) = RWORK( K ) + S
00386   160             CONTINUE
00387                ELSE
00388                   DO 180 K = 1, N
00389                      S = CABS1( X( K, J ) )
00390                      DO 170 I = K + 1, N
00391                         S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
00392   170                CONTINUE
00393                      RWORK( K ) = RWORK( K ) + S
00394   180             CONTINUE
00395                END IF
00396             END IF
00397          END IF
00398          S = ZERO
00399          DO 190 I = 1, N
00400             IF( RWORK( I ).GT.SAFE2 ) THEN
00401                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
00402             ELSE
00403                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
00404      $             ( RWORK( I )+SAFE1 ) )
00405             END IF
00406   190    CONTINUE
00407          BERR( J ) = S
00408 *
00409 *        Bound error from formula
00410 *
00411 *        norm(X - XTRUE) / norm(X) .le. FERR =
00412 *        norm( abs(inv(op(A)))*
00413 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
00414 *
00415 *        where
00416 *          norm(Z) is the magnitude of the largest component of Z
00417 *          inv(op(A)) is the inverse of op(A)
00418 *          abs(Z) is the componentwise absolute value of the matrix or
00419 *             vector Z
00420 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00421 *          EPS is machine epsilon
00422 *
00423 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
00424 *        is incremented by SAFE1 if the i-th component of
00425 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
00426 *
00427 *        Use ZLACN2 to estimate the infinity-norm of the matrix
00428 *           inv(op(A)) * diag(W),
00429 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
00430 *
00431          DO 200 I = 1, N
00432             IF( RWORK( I ).GT.SAFE2 ) THEN
00433                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
00434             ELSE
00435                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
00436      $                      SAFE1
00437             END IF
00438   200    CONTINUE
00439 *
00440          KASE = 0
00441   210    CONTINUE
00442          CALL ZLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
00443          IF( KASE.NE.0 ) THEN
00444             IF( KASE.EQ.1 ) THEN
00445 *
00446 *              Multiply by diag(W)*inv(op(A)**H).
00447 *
00448                CALL ZTRSV( UPLO, TRANST, DIAG, N, A, LDA, WORK, 1 )
00449                DO 220 I = 1, N
00450                   WORK( I ) = RWORK( I )*WORK( I )
00451   220          CONTINUE
00452             ELSE
00453 *
00454 *              Multiply by inv(op(A))*diag(W).
00455 *
00456                DO 230 I = 1, N
00457                   WORK( I ) = RWORK( I )*WORK( I )
00458   230          CONTINUE
00459                CALL ZTRSV( UPLO, TRANSN, DIAG, N, A, LDA, WORK, 1 )
00460             END IF
00461             GO TO 210
00462          END IF
00463 *
00464 *        Normalize error.
00465 *
00466          LSTRES = ZERO
00467          DO 240 I = 1, N
00468             LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
00469   240    CONTINUE
00470          IF( LSTRES.NE.ZERO )
00471      $      FERR( J ) = FERR( J ) / LSTRES
00472 *
00473   250 CONTINUE
00474 *
00475       RETURN
00476 *
00477 *     End of ZTRRFS
00478 *
00479       END
 All Files Functions