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