![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DTPT03 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 DTPT03( 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 * DOUBLE PRECISION RESID, SCALE, TSCAL 00018 * .. 00019 * .. Array Arguments .. 00020 * DOUBLE PRECISION AP( * ), B( LDB, * ), CNORM( * ), WORK( * ), 00021 * $ X( LDX, * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> DTPT03 computes the residual for the solution to a scaled triangular 00031 *> system of equations A*x = s*b or A'*x = s*b when the triangular 00032 *> matrix A is stored in packed format. Here A' is the transpose of A, 00033 *> s is a scalar, and x and b are N by NRHS matrices. The test ratio is 00034 *> the 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] AP 00081 *> \verbatim 00082 *> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2) 00083 *> The upper or lower triangular matrix A, packed columnwise in 00084 *> a linear array. The j-th column of A is stored in the array 00085 *> AP as follows: 00086 *> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j; 00087 *> if UPLO = 'L', 00088 *> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] SCALE 00092 *> \verbatim 00093 *> SCALE is DOUBLE PRECISION 00094 *> The scaling factor s used in solving the triangular system. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] CNORM 00098 *> \verbatim 00099 *> CNORM is DOUBLE PRECISION array, dimension (N) 00100 *> The 1-norms of the columns of A, not counting the diagonal. 00101 *> \endverbatim 00102 *> 00103 *> \param[in] TSCAL 00104 *> \verbatim 00105 *> TSCAL is DOUBLE PRECISION 00106 *> The scaling factor used in computing the 1-norms in CNORM. 00107 *> CNORM actually contains the column norms of TSCAL*A. 00108 *> \endverbatim 00109 *> 00110 *> \param[in] X 00111 *> \verbatim 00112 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 00113 *> The computed solution vectors for the system of linear 00114 *> equations. 00115 *> \endverbatim 00116 *> 00117 *> \param[in] LDX 00118 *> \verbatim 00119 *> LDX is INTEGER 00120 *> The leading dimension of the array X. LDX >= max(1,N). 00121 *> \endverbatim 00122 *> 00123 *> \param[in] B 00124 *> \verbatim 00125 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00126 *> The right hand side vectors for the system of linear 00127 *> equations. 00128 *> \endverbatim 00129 *> 00130 *> \param[in] LDB 00131 *> \verbatim 00132 *> LDB is INTEGER 00133 *> The leading dimension of the array B. LDB >= max(1,N). 00134 *> \endverbatim 00135 *> 00136 *> \param[out] WORK 00137 *> \verbatim 00138 *> WORK is DOUBLE PRECISION array, dimension (N) 00139 *> \endverbatim 00140 *> 00141 *> \param[out] RESID 00142 *> \verbatim 00143 *> RESID is DOUBLE PRECISION 00144 *> The maximum over the number of right hand sides of 00145 *> norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ). 00146 *> \endverbatim 00147 * 00148 * Authors: 00149 * ======== 00150 * 00151 *> \author Univ. of Tennessee 00152 *> \author Univ. of California Berkeley 00153 *> \author Univ. of Colorado Denver 00154 *> \author NAG Ltd. 00155 * 00156 *> \date November 2011 00157 * 00158 *> \ingroup double_lin 00159 * 00160 * ===================================================================== 00161 SUBROUTINE DTPT03( UPLO, TRANS, DIAG, N, NRHS, AP, SCALE, CNORM, 00162 $ TSCAL, X, LDX, B, LDB, WORK, RESID ) 00163 * 00164 * -- LAPACK test routine (version 3.4.0) -- 00165 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00167 * November 2011 00168 * 00169 * .. Scalar Arguments .. 00170 CHARACTER DIAG, TRANS, UPLO 00171 INTEGER LDB, LDX, N, NRHS 00172 DOUBLE PRECISION RESID, SCALE, TSCAL 00173 * .. 00174 * .. Array Arguments .. 00175 DOUBLE PRECISION AP( * ), B( LDB, * ), CNORM( * ), WORK( * ), 00176 $ X( LDX, * ) 00177 * .. 00178 * 00179 * ===================================================================== 00180 * 00181 * .. Parameters .. 00182 DOUBLE PRECISION ONE, ZERO 00183 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00184 * .. 00185 * .. Local Scalars .. 00186 INTEGER IX, J, JJ 00187 DOUBLE PRECISION BIGNUM, EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL 00188 * .. 00189 * .. External Functions .. 00190 LOGICAL LSAME 00191 INTEGER IDAMAX 00192 DOUBLE PRECISION DLAMCH 00193 EXTERNAL LSAME, IDAMAX, DLAMCH 00194 * .. 00195 * .. External Subroutines .. 00196 EXTERNAL DAXPY, DCOPY, DLABAD, DSCAL, DTPMV 00197 * .. 00198 * .. Intrinsic Functions .. 00199 INTRINSIC ABS, DBLE, MAX 00200 * .. 00201 * .. Executable Statements .. 00202 * 00203 * Quick exit if N = 0. 00204 * 00205 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 00206 RESID = ZERO 00207 RETURN 00208 END IF 00209 EPS = DLAMCH( 'Epsilon' ) 00210 SMLNUM = DLAMCH( 'Safe minimum' ) 00211 BIGNUM = ONE / SMLNUM 00212 CALL DLABAD( SMLNUM, BIGNUM ) 00213 * 00214 * Compute the norm of the triangular matrix A using the column 00215 * norms already computed by DLATPS. 00216 * 00217 TNORM = ZERO 00218 IF( LSAME( DIAG, 'N' ) ) THEN 00219 IF( LSAME( UPLO, 'U' ) ) THEN 00220 JJ = 1 00221 DO 10 J = 1, N 00222 TNORM = MAX( TNORM, TSCAL*ABS( AP( JJ ) )+CNORM( J ) ) 00223 JJ = JJ + J + 1 00224 10 CONTINUE 00225 ELSE 00226 JJ = 1 00227 DO 20 J = 1, N 00228 TNORM = MAX( TNORM, TSCAL*ABS( AP( JJ ) )+CNORM( J ) ) 00229 JJ = JJ + N - J + 1 00230 20 CONTINUE 00231 END IF 00232 ELSE 00233 DO 30 J = 1, N 00234 TNORM = MAX( TNORM, TSCAL+CNORM( J ) ) 00235 30 CONTINUE 00236 END IF 00237 * 00238 * Compute the maximum over the number of right hand sides of 00239 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ). 00240 * 00241 RESID = ZERO 00242 DO 40 J = 1, NRHS 00243 CALL DCOPY( N, X( 1, J ), 1, WORK, 1 ) 00244 IX = IDAMAX( N, WORK, 1 ) 00245 XNORM = MAX( ONE, ABS( X( IX, J ) ) ) 00246 XSCAL = ( ONE / XNORM ) / DBLE( N ) 00247 CALL DSCAL( N, XSCAL, WORK, 1 ) 00248 CALL DTPMV( UPLO, TRANS, DIAG, N, AP, WORK, 1 ) 00249 CALL DAXPY( N, -SCALE*XSCAL, B( 1, J ), 1, WORK, 1 ) 00250 IX = IDAMAX( N, WORK, 1 ) 00251 ERR = TSCAL*ABS( WORK( IX ) ) 00252 IX = IDAMAX( N, X( 1, J ), 1 ) 00253 XNORM = ABS( X( IX, J ) ) 00254 IF( ERR*SMLNUM.LE.XNORM ) THEN 00255 IF( XNORM.GT.ZERO ) 00256 $ ERR = ERR / XNORM 00257 ELSE 00258 IF( ERR.GT.ZERO ) 00259 $ ERR = ONE / EPS 00260 END IF 00261 IF( ERR*SMLNUM.LE.TNORM ) THEN 00262 IF( TNORM.GT.ZERO ) 00263 $ ERR = ERR / TNORM 00264 ELSE 00265 IF( ERR.GT.ZERO ) 00266 $ ERR = ONE / EPS 00267 END IF 00268 RESID = MAX( RESID, ERR ) 00269 40 CONTINUE 00270 * 00271 RETURN 00272 * 00273 * End of DTPT03 00274 * 00275 END