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