LAPACK: SRC/stprfs.f Source File

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