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