LAPACK: SRC/ctprfs.f Source File

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