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