LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sppt05.f
Go to the documentation of this file.
00001 *> \brief \b SPPT05
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE SPPT05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT,
00012 *                          LDXACT, FERR, BERR, RESLTS )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            LDB, LDX, LDXACT, N, NRHS
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       REAL               AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
00020 *      $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> SPPT05 tests the error bounds from iterative refinement for the
00030 *> computed solution to a system of equations A*X = B, where A is a
00031 *> symmetric matrix in packed storage format.
00032 *>
00033 *> RESLTS(1) = test of the error bound
00034 *>           = norm(X - XACT) / ( norm(X) * FERR )
00035 *>
00036 *> A large value is returned if this ratio is not less than one.
00037 *>
00038 *> RESLTS(2) = residual from the iterative refinement routine
00039 *>           = the maximum of BERR / ( (n+1)*EPS + (*) ), where
00040 *>             (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
00041 *> \endverbatim
00042 *
00043 *  Arguments:
00044 *  ==========
00045 *
00046 *> \param[in] UPLO
00047 *> \verbatim
00048 *>          UPLO is CHARACTER*1
00049 *>          Specifies whether the upper or lower triangular part of the
00050 *>          symmetric matrix A is stored.
00051 *>          = 'U':  Upper triangular
00052 *>          = 'L':  Lower triangular
00053 *> \endverbatim
00054 *>
00055 *> \param[in] N
00056 *> \verbatim
00057 *>          N is INTEGER
00058 *>          The number of rows of the matrices X, B, and XACT, and the
00059 *>          order of the matrix A.  N >= 0.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] NRHS
00063 *> \verbatim
00064 *>          NRHS is INTEGER
00065 *>          The number of columns of the matrices X, B, and XACT.
00066 *>          NRHS >= 0.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] AP
00070 *> \verbatim
00071 *>          AP is REAL array, dimension (N*(N+1)/2)
00072 *>          The upper or lower triangle of the symmetric matrix A, packed
00073 *>          columnwise in a linear array.  The j-th column of A is stored
00074 *>          in the array AP as follows:
00075 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00076 *>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] B
00080 *> \verbatim
00081 *>          B is REAL array, dimension (LDB,NRHS)
00082 *>          The right hand side vectors for the system of linear
00083 *>          equations.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] LDB
00087 *> \verbatim
00088 *>          LDB is INTEGER
00089 *>          The leading dimension of the array B.  LDB >= max(1,N).
00090 *> \endverbatim
00091 *>
00092 *> \param[in] X
00093 *> \verbatim
00094 *>          X is REAL array, dimension (LDX,NRHS)
00095 *>          The computed solution vectors.  Each vector is stored as a
00096 *>          column of the matrix X.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] LDX
00100 *> \verbatim
00101 *>          LDX is INTEGER
00102 *>          The leading dimension of the array X.  LDX >= max(1,N).
00103 *> \endverbatim
00104 *>
00105 *> \param[in] XACT
00106 *> \verbatim
00107 *>          XACT is REAL array, dimension (LDX,NRHS)
00108 *>          The exact solution vectors.  Each vector is stored as a
00109 *>          column of the matrix XACT.
00110 *> \endverbatim
00111 *>
00112 *> \param[in] LDXACT
00113 *> \verbatim
00114 *>          LDXACT is INTEGER
00115 *>          The leading dimension of the array XACT.  LDXACT >= max(1,N).
00116 *> \endverbatim
00117 *>
00118 *> \param[in] FERR
00119 *> \verbatim
00120 *>          FERR is REAL array, dimension (NRHS)
00121 *>          The estimated forward error bounds for each solution vector
00122 *>          X.  If XTRUE is the true solution, FERR bounds the magnitude
00123 *>          of the largest entry in (X - XTRUE) divided by the magnitude
00124 *>          of the largest entry in X.
00125 *> \endverbatim
00126 *>
00127 *> \param[in] BERR
00128 *> \verbatim
00129 *>          BERR is REAL array, dimension (NRHS)
00130 *>          The componentwise relative backward error of each solution
00131 *>          vector (i.e., the smallest relative change in any entry of A
00132 *>          or B that makes X an exact solution).
00133 *> \endverbatim
00134 *>
00135 *> \param[out] RESLTS
00136 *> \verbatim
00137 *>          RESLTS is REAL array, dimension (2)
00138 *>          The maximum over the NRHS solution vectors of the ratios:
00139 *>          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
00140 *>          RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
00141 *> \endverbatim
00142 *
00143 *  Authors:
00144 *  ========
00145 *
00146 *> \author Univ. of Tennessee 
00147 *> \author Univ. of California Berkeley 
00148 *> \author Univ. of Colorado Denver 
00149 *> \author NAG Ltd. 
00150 *
00151 *> \date November 2011
00152 *
00153 *> \ingroup single_lin
00154 *
00155 *  =====================================================================
00156       SUBROUTINE SPPT05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT,
00157      $                   LDXACT, FERR, BERR, RESLTS )
00158 *
00159 *  -- LAPACK test routine (version 3.4.0) --
00160 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00161 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00162 *     November 2011
00163 *
00164 *     .. Scalar Arguments ..
00165       CHARACTER          UPLO
00166       INTEGER            LDB, LDX, LDXACT, N, NRHS
00167 *     ..
00168 *     .. Array Arguments ..
00169       REAL               AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
00170      $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
00171 *     ..
00172 *
00173 *  =====================================================================
00174 *
00175 *     .. Parameters ..
00176       REAL               ZERO, ONE
00177       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00178 *     ..
00179 *     .. Local Scalars ..
00180       LOGICAL            UPPER
00181       INTEGER            I, IMAX, J, JC, K
00182       REAL               AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
00183 *     ..
00184 *     .. External Functions ..
00185       LOGICAL            LSAME
00186       INTEGER            ISAMAX
00187       REAL               SLAMCH
00188       EXTERNAL           LSAME, ISAMAX, SLAMCH
00189 *     ..
00190 *     .. Intrinsic Functions ..
00191       INTRINSIC          ABS, MAX, MIN
00192 *     ..
00193 *     .. Executable Statements ..
00194 *
00195 *     Quick exit if N = 0 or NRHS = 0.
00196 *
00197       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
00198          RESLTS( 1 ) = ZERO
00199          RESLTS( 2 ) = ZERO
00200          RETURN
00201       END IF
00202 *
00203       EPS = SLAMCH( 'Epsilon' )
00204       UNFL = SLAMCH( 'Safe minimum' )
00205       OVFL = ONE / UNFL
00206       UPPER = LSAME( UPLO, 'U' )
00207 *
00208 *     Test 1:  Compute the maximum of
00209 *        norm(X - XACT) / ( norm(X) * FERR )
00210 *     over all the vectors X and XACT using the infinity-norm.
00211 *
00212       ERRBND = ZERO
00213       DO 30 J = 1, NRHS
00214          IMAX = ISAMAX( N, X( 1, J ), 1 )
00215          XNORM = MAX( ABS( X( IMAX, J ) ), UNFL )
00216          DIFF = ZERO
00217          DO 10 I = 1, N
00218             DIFF = MAX( DIFF, ABS( X( I, J )-XACT( I, J ) ) )
00219    10    CONTINUE
00220 *
00221          IF( XNORM.GT.ONE ) THEN
00222             GO TO 20
00223          ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
00224             GO TO 20
00225          ELSE
00226             ERRBND = ONE / EPS
00227             GO TO 30
00228          END IF
00229 *
00230    20    CONTINUE
00231          IF( DIFF / XNORM.LE.FERR( J ) ) THEN
00232             ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
00233          ELSE
00234             ERRBND = ONE / EPS
00235          END IF
00236    30 CONTINUE
00237       RESLTS( 1 ) = ERRBND
00238 *
00239 *     Test 2:  Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
00240 *     (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
00241 *
00242       DO 90 K = 1, NRHS
00243          DO 80 I = 1, N
00244             TMP = ABS( B( I, K ) )
00245             IF( UPPER ) THEN
00246                JC = ( ( I-1 )*I ) / 2
00247                DO 40 J = 1, I
00248                   TMP = TMP + ABS( AP( JC+J ) )*ABS( X( J, K ) )
00249    40          CONTINUE
00250                JC = JC + I
00251                DO 50 J = I + 1, N
00252                   TMP = TMP + ABS( AP( JC ) )*ABS( X( J, K ) )
00253                   JC = JC + J
00254    50          CONTINUE
00255             ELSE
00256                JC = I
00257                DO 60 J = 1, I - 1
00258                   TMP = TMP + ABS( AP( JC ) )*ABS( X( J, K ) )
00259                   JC = JC + N - J
00260    60          CONTINUE
00261                DO 70 J = I, N
00262                   TMP = TMP + ABS( AP( JC+J-I ) )*ABS( X( J, K ) )
00263    70          CONTINUE
00264             END IF
00265             IF( I.EQ.1 ) THEN
00266                AXBI = TMP
00267             ELSE
00268                AXBI = MIN( AXBI, TMP )
00269             END IF
00270    80    CONTINUE
00271          TMP = BERR( K ) / ( ( N+1 )*EPS+( N+1 )*UNFL /
00272      $         MAX( AXBI, ( N+1 )*UNFL ) )
00273          IF( K.EQ.1 ) THEN
00274             RESLTS( 2 ) = TMP
00275          ELSE
00276             RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
00277          END IF
00278    90 CONTINUE
00279 *
00280       RETURN
00281 *
00282 *     End of SPPT05
00283 *
00284       END
 All Files Functions