LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
strt03.f
Go to the documentation of this file.
00001 *> \brief \b STRT03
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 STRT03( UPLO, TRANS, DIAG, N, NRHS, A, LDA, SCALE,
00012 *                          CNORM, TSCAL, X, LDX, B, LDB, WORK, RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          DIAG, TRANS, UPLO
00016 *       INTEGER            LDA, LDB, LDX, N, NRHS
00017 *       REAL               RESID, SCALE, TSCAL
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       REAL               A( LDA, * ), B( LDB, * ), CNORM( * ),
00021 *      $                   WORK( * ), X( LDX, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> STRT03 computes the residual for the solution to a scaled triangular
00031 *> system of equations A*x = s*b  or  A'*x = s*b.
00032 *> Here A is a triangular matrix, A' is the transpose of A, s is a
00033 *> scalar, and x and b are N by NRHS matrices.  The test ratio is the
00034 *> maximum over the number of right hand sides of
00035 *>    norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
00036 *> where op(A) denotes A or A' and EPS is the machine epsilon.
00037 *> \endverbatim
00038 *
00039 *  Arguments:
00040 *  ==========
00041 *
00042 *> \param[in] UPLO
00043 *> \verbatim
00044 *>          UPLO is CHARACTER*1
00045 *>          Specifies whether the matrix A is upper or lower triangular.
00046 *>          = 'U':  Upper triangular
00047 *>          = 'L':  Lower triangular
00048 *> \endverbatim
00049 *>
00050 *> \param[in] TRANS
00051 *> \verbatim
00052 *>          TRANS is CHARACTER*1
00053 *>          Specifies the operation applied to A.
00054 *>          = 'N':  A *x = s*b  (No transpose)
00055 *>          = 'T':  A'*x = s*b  (Transpose)
00056 *>          = 'C':  A'*x = s*b  (Conjugate transpose = Transpose)
00057 *> \endverbatim
00058 *>
00059 *> \param[in] DIAG
00060 *> \verbatim
00061 *>          DIAG is CHARACTER*1
00062 *>          Specifies whether or not the matrix A is unit triangular.
00063 *>          = 'N':  Non-unit triangular
00064 *>          = 'U':  Unit triangular
00065 *> \endverbatim
00066 *>
00067 *> \param[in] N
00068 *> \verbatim
00069 *>          N is INTEGER
00070 *>          The order of the matrix A.  N >= 0.
00071 *> \endverbatim
00072 *>
00073 *> \param[in] NRHS
00074 *> \verbatim
00075 *>          NRHS is INTEGER
00076 *>          The number of right hand sides, i.e., the number of columns
00077 *>          of the matrices X and B.  NRHS >= 0.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] A
00081 *> \verbatim
00082 *>          A is REAL array, dimension (LDA,N)
00083 *>          The triangular matrix A.  If UPLO = 'U', the leading n by n
00084 *>          upper triangular part of the array A contains the upper
00085 *>          triangular matrix, and the strictly lower triangular part of
00086 *>          A is not referenced.  If UPLO = 'L', the leading n by n lower
00087 *>          triangular part of the array A contains the lower triangular
00088 *>          matrix, and the strictly upper triangular part of A is not
00089 *>          referenced.  If DIAG = 'U', the diagonal elements of A are
00090 *>          also not referenced and are assumed to be 1.
00091 *> \endverbatim
00092 *>
00093 *> \param[in] LDA
00094 *> \verbatim
00095 *>          LDA is INTEGER
00096 *>          The leading dimension of the array A.  LDA >= max(1,N).
00097 *> \endverbatim
00098 *>
00099 *> \param[in] SCALE
00100 *> \verbatim
00101 *>          SCALE is REAL
00102 *>          The scaling factor s used in solving the triangular system.
00103 *> \endverbatim
00104 *>
00105 *> \param[in] CNORM
00106 *> \verbatim
00107 *>          CNORM is REAL array, dimension (N)
00108 *>          The 1-norms of the columns of A, not counting the diagonal.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] TSCAL
00112 *> \verbatim
00113 *>          TSCAL is REAL
00114 *>          The scaling factor used in computing the 1-norms in CNORM.
00115 *>          CNORM actually contains the column norms of TSCAL*A.
00116 *> \endverbatim
00117 *>
00118 *> \param[in] X
00119 *> \verbatim
00120 *>          X is REAL array, dimension (LDX,NRHS)
00121 *>          The computed solution vectors for the system of linear
00122 *>          equations.
00123 *> \endverbatim
00124 *>
00125 *> \param[in] LDX
00126 *> \verbatim
00127 *>          LDX is INTEGER
00128 *>          The leading dimension of the array X.  LDX >= max(1,N).
00129 *> \endverbatim
00130 *>
00131 *> \param[in] B
00132 *> \verbatim
00133 *>          B is REAL array, dimension (LDB,NRHS)
00134 *>          The right hand side vectors for the system of linear
00135 *>          equations.
00136 *> \endverbatim
00137 *>
00138 *> \param[in] LDB
00139 *> \verbatim
00140 *>          LDB is INTEGER
00141 *>          The leading dimension of the array B.  LDB >= max(1,N).
00142 *> \endverbatim
00143 *>
00144 *> \param[out] WORK
00145 *> \verbatim
00146 *>          WORK is REAL array, dimension (N)
00147 *> \endverbatim
00148 *>
00149 *> \param[out] RESID
00150 *> \verbatim
00151 *>          RESID is REAL
00152 *>          The maximum over the number of right hand sides of
00153 *>          norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
00154 *> \endverbatim
00155 *
00156 *  Authors:
00157 *  ========
00158 *
00159 *> \author Univ. of Tennessee 
00160 *> \author Univ. of California Berkeley 
00161 *> \author Univ. of Colorado Denver 
00162 *> \author NAG Ltd. 
00163 *
00164 *> \date November 2011
00165 *
00166 *> \ingroup single_lin
00167 *
00168 *  =====================================================================
00169       SUBROUTINE STRT03( UPLO, TRANS, DIAG, N, NRHS, A, LDA, SCALE,
00170      $                   CNORM, TSCAL, X, LDX, B, LDB, WORK, RESID )
00171 *
00172 *  -- LAPACK test routine (version 3.4.0) --
00173 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00174 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00175 *     November 2011
00176 *
00177 *     .. Scalar Arguments ..
00178       CHARACTER          DIAG, TRANS, UPLO
00179       INTEGER            LDA, LDB, LDX, N, NRHS
00180       REAL               RESID, SCALE, TSCAL
00181 *     ..
00182 *     .. Array Arguments ..
00183       REAL               A( LDA, * ), B( LDB, * ), CNORM( * ),
00184      $                   WORK( * ), X( LDX, * )
00185 *     ..
00186 *
00187 *  =====================================================================
00188 *
00189 *     .. Parameters ..
00190       REAL               ONE, ZERO
00191       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00192 *     ..
00193 *     .. Local Scalars ..
00194       INTEGER            IX, J
00195       REAL               BIGNUM, EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
00196 *     ..
00197 *     .. External Functions ..
00198       LOGICAL            LSAME
00199       INTEGER            ISAMAX
00200       REAL               SLAMCH
00201       EXTERNAL           LSAME, ISAMAX, SLAMCH
00202 *     ..
00203 *     .. External Subroutines ..
00204       EXTERNAL           SAXPY, SCOPY, SLABAD, SSCAL, STRMV
00205 *     ..
00206 *     .. Intrinsic Functions ..
00207       INTRINSIC          ABS, MAX, REAL
00208 *     ..
00209 *     .. Executable Statements ..
00210 *
00211 *     Quick exit if N = 0
00212 *
00213       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
00214          RESID = ZERO
00215          RETURN
00216       END IF
00217       EPS = SLAMCH( 'Epsilon' )
00218       SMLNUM = SLAMCH( 'Safe minimum' )
00219       BIGNUM = ONE / SMLNUM
00220       CALL SLABAD( SMLNUM, BIGNUM )
00221 *
00222 *     Compute the norm of the triangular matrix A using the column
00223 *     norms already computed by SLATRS.
00224 *
00225       TNORM = ZERO
00226       IF( LSAME( DIAG, 'N' ) ) THEN
00227          DO 10 J = 1, N
00228             TNORM = MAX( TNORM, TSCAL*ABS( A( J, J ) )+CNORM( J ) )
00229    10    CONTINUE
00230       ELSE
00231          DO 20 J = 1, N
00232             TNORM = MAX( TNORM, TSCAL+CNORM( J ) )
00233    20    CONTINUE
00234       END IF
00235 *
00236 *     Compute the maximum over the number of right hand sides of
00237 *        norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
00238 *
00239       RESID = ZERO
00240       DO 30 J = 1, NRHS
00241          CALL SCOPY( N, X( 1, J ), 1, WORK, 1 )
00242          IX = ISAMAX( N, WORK, 1 )
00243          XNORM = MAX( ONE, ABS( X( IX, J ) ) )
00244          XSCAL = ( ONE / XNORM ) / REAL( N )
00245          CALL SSCAL( N, XSCAL, WORK, 1 )
00246          CALL STRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
00247          CALL SAXPY( N, -SCALE*XSCAL, B( 1, J ), 1, WORK, 1 )
00248          IX = ISAMAX( N, WORK, 1 )
00249          ERR = TSCAL*ABS( WORK( IX ) )
00250          IX = ISAMAX( N, X( 1, J ), 1 )
00251          XNORM = ABS( X( IX, J ) )
00252          IF( ERR*SMLNUM.LE.XNORM ) THEN
00253             IF( XNORM.GT.ZERO )
00254      $         ERR = ERR / XNORM
00255          ELSE
00256             IF( ERR.GT.ZERO )
00257      $         ERR = ONE / EPS
00258          END IF
00259          IF( ERR*SMLNUM.LE.TNORM ) THEN
00260             IF( TNORM.GT.ZERO )
00261      $         ERR = ERR / TNORM
00262          ELSE
00263             IF( ERR.GT.ZERO )
00264      $         ERR = ONE / EPS
00265          END IF
00266          RESID = MAX( RESID, ERR )
00267    30 CONTINUE
00268 *
00269       RETURN
00270 *
00271 *     End of STRT03
00272 *
00273       END
 All Files Functions