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