![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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