|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT ) 00002 * 00003 * -- LAPACK test routine (version 3.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2006 00006 * 00007 * .. Scalar Arguments .. 00008 INTEGER KNT, LMAX, NINFO 00009 DOUBLE PRECISION RMAX 00010 * .. 00011 * 00012 * Purpose 00013 * ======= 00014 * 00015 * DGET39 tests DLAQTR, a routine for solving the real or 00016 * special complex quasi upper triangular system 00017 * 00018 * op(T)*p = scale*c, 00019 * or 00020 * op(T + iB)*(p+iq) = scale*(c+id), 00021 * 00022 * in real arithmetic. T is upper quasi-triangular. 00023 * If it is complex, then the first diagonal block of T must be 00024 * 1 by 1, B has the special structure 00025 * 00026 * B = [ b(1) b(2) ... b(n) ] 00027 * [ w ] 00028 * [ w ] 00029 * [ . ] 00030 * [ w ] 00031 * 00032 * op(A) = A or A', where A' denotes the conjugate transpose of 00033 * the matrix A. 00034 * 00035 * On input, X = [ c ]. On output, X = [ p ]. 00036 * [ d ] [ q ] 00037 * 00038 * Scale is an output less than or equal to 1, chosen to avoid 00039 * overflow in X. 00040 * This subroutine is specially designed for the condition number 00041 * estimation in the eigenproblem routine DTRSNA. 00042 * 00043 * The test code verifies that the following residual is order 1: 00044 * 00045 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| 00046 * ----------------------------------------- 00047 * max(ulp*(||T||+||B||)*(||x1||+||x2||), 00048 * (||T||+||B||)*smlnum/ulp, 00049 * smlnum) 00050 * 00051 * (The (||T||+||B||)*smlnum/ulp term accounts for possible 00052 * (gradual or nongradual) underflow in x1 and x2.) 00053 * 00054 * Arguments 00055 * ========== 00056 * 00057 * RMAX (output) DOUBLE PRECISION 00058 * Value of the largest test ratio. 00059 * 00060 * LMAX (output) INTEGER 00061 * Example number where largest test ratio achieved. 00062 * 00063 * NINFO (output) INTEGER 00064 * Number of examples where INFO is nonzero. 00065 * 00066 * KNT (output) INTEGER 00067 * Total number of examples tested. 00068 * 00069 * ===================================================================== 00070 * 00071 * .. Parameters .. 00072 INTEGER LDT, LDT2 00073 PARAMETER ( LDT = 10, LDT2 = 2*LDT ) 00074 DOUBLE PRECISION ZERO, ONE 00075 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00076 * .. 00077 * .. Local Scalars .. 00078 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N, 00079 $ NDIM 00080 DOUBLE PRECISION BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID, 00081 $ SCALE, SMLNUM, W, XNORM 00082 * .. 00083 * .. External Functions .. 00084 INTEGER IDAMAX 00085 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE 00086 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE 00087 * .. 00088 * .. External Subroutines .. 00089 EXTERNAL DCOPY, DGEMV, DLABAD, DLAQTR 00090 * .. 00091 * .. Intrinsic Functions .. 00092 INTRINSIC ABS, COS, DBLE, MAX, SIN, SQRT 00093 * .. 00094 * .. Local Arrays .. 00095 INTEGER IDIM( 6 ), IVAL( 5, 5, 6 ) 00096 DOUBLE PRECISION B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ), 00097 $ VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ), 00098 $ VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 ) 00099 * .. 00100 * .. Data statements .. 00101 DATA IDIM / 4, 5*5 / 00102 DATA IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0, 00103 $ 4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4, 00104 $ 0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2, 00105 $ 0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1, 00106 $ 4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2, 00107 $ -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0, 00108 $ 5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1, 00109 $ 4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 / 00110 * .. 00111 * .. Executable Statements .. 00112 * 00113 * Get machine parameters 00114 * 00115 EPS = DLAMCH( 'P' ) 00116 SMLNUM = DLAMCH( 'S' ) 00117 BIGNUM = ONE / SMLNUM 00118 CALL DLABAD( SMLNUM, BIGNUM ) 00119 * 00120 * Set up test case parameters 00121 * 00122 VM1( 1 ) = ONE 00123 VM1( 2 ) = SQRT( SMLNUM ) 00124 VM1( 3 ) = SQRT( VM1( 2 ) ) 00125 VM1( 4 ) = SQRT( BIGNUM ) 00126 VM1( 5 ) = SQRT( VM1( 4 ) ) 00127 * 00128 VM2( 1 ) = ONE 00129 VM2( 2 ) = SQRT( SMLNUM ) 00130 VM2( 3 ) = SQRT( VM2( 2 ) ) 00131 VM2( 4 ) = SQRT( BIGNUM ) 00132 VM2( 5 ) = SQRT( VM2( 4 ) ) 00133 * 00134 VM3( 1 ) = ONE 00135 VM3( 2 ) = SQRT( SMLNUM ) 00136 VM3( 3 ) = SQRT( VM3( 2 ) ) 00137 VM3( 4 ) = SQRT( BIGNUM ) 00138 VM3( 5 ) = SQRT( VM3( 4 ) ) 00139 * 00140 VM4( 1 ) = ONE 00141 VM4( 2 ) = SQRT( SMLNUM ) 00142 VM4( 3 ) = SQRT( VM4( 2 ) ) 00143 VM4( 4 ) = SQRT( BIGNUM ) 00144 VM4( 5 ) = SQRT( VM4( 4 ) ) 00145 * 00146 VM5( 1 ) = ONE 00147 VM5( 2 ) = EPS 00148 VM5( 3 ) = SQRT( SMLNUM ) 00149 * 00150 * Initalization 00151 * 00152 KNT = 0 00153 RMAX = ZERO 00154 NINFO = 0 00155 SMLNUM = SMLNUM / EPS 00156 * 00157 * Begin test loop 00158 * 00159 DO 140 IVM5 = 1, 3 00160 DO 130 IVM4 = 1, 5 00161 DO 120 IVM3 = 1, 5 00162 DO 110 IVM2 = 1, 5 00163 DO 100 IVM1 = 1, 5 00164 DO 90 NDIM = 1, 6 00165 * 00166 N = IDIM( NDIM ) 00167 DO 20 I = 1, N 00168 DO 10 J = 1, N 00169 T( I, J ) = DBLE( IVAL( I, J, NDIM ) )* 00170 $ VM1( IVM1 ) 00171 IF( I.GE.J ) 00172 $ T( I, J ) = T( I, J )*VM5( IVM5 ) 00173 10 CONTINUE 00174 20 CONTINUE 00175 * 00176 W = ONE*VM2( IVM2 ) 00177 * 00178 DO 30 I = 1, N 00179 B( I ) = COS( DBLE( I ) )*VM3( IVM3 ) 00180 30 CONTINUE 00181 * 00182 DO 40 I = 1, 2*N 00183 D( I ) = SIN( DBLE( I ) )*VM4( IVM4 ) 00184 40 CONTINUE 00185 * 00186 NORM = DLANGE( '1', N, N, T, LDT, WORK ) 00187 K = IDAMAX( N, B, 1 ) 00188 NORMTB = NORM + ABS( B( K ) ) + ABS( W ) 00189 * 00190 CALL DCOPY( N, D, 1, X, 1 ) 00191 KNT = KNT + 1 00192 CALL DLAQTR( .FALSE., .TRUE., N, T, LDT, DUM, 00193 $ DUMM, SCALE, X, WORK, INFO ) 00194 IF( INFO.NE.0 ) 00195 $ NINFO = NINFO + 1 00196 * 00197 * || T*x - scale*d || / 00198 * max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum) 00199 * 00200 CALL DCOPY( N, D, 1, Y, 1 ) 00201 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT, 00202 $ X, 1, -SCALE, Y, 1 ) 00203 XNORM = DASUM( N, X, 1 ) 00204 RESID = DASUM( N, Y, 1 ) 00205 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM, 00206 $ ( NORM*EPS )*XNORM ) 00207 RESID = RESID / DOMIN 00208 IF( RESID.GT.RMAX ) THEN 00209 RMAX = RESID 00210 LMAX = KNT 00211 END IF 00212 * 00213 CALL DCOPY( N, D, 1, X, 1 ) 00214 KNT = KNT + 1 00215 CALL DLAQTR( .TRUE., .TRUE., N, T, LDT, DUM, 00216 $ DUMM, SCALE, X, WORK, INFO ) 00217 IF( INFO.NE.0 ) 00218 $ NINFO = NINFO + 1 00219 * 00220 * || T*x - scale*d || / 00221 * max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum) 00222 * 00223 CALL DCOPY( N, D, 1, Y, 1 ) 00224 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X, 00225 $ 1, -SCALE, Y, 1 ) 00226 XNORM = DASUM( N, X, 1 ) 00227 RESID = DASUM( N, Y, 1 ) 00228 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM, 00229 $ ( NORM*EPS )*XNORM ) 00230 RESID = RESID / DOMIN 00231 IF( RESID.GT.RMAX ) THEN 00232 RMAX = RESID 00233 LMAX = KNT 00234 END IF 00235 * 00236 CALL DCOPY( 2*N, D, 1, X, 1 ) 00237 KNT = KNT + 1 00238 CALL DLAQTR( .FALSE., .FALSE., N, T, LDT, B, W, 00239 $ SCALE, X, WORK, INFO ) 00240 IF( INFO.NE.0 ) 00241 $ NINFO = NINFO + 1 00242 * 00243 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| / 00244 * max(ulp*(||T||+||B||)*(||x1||+||x2||), 00245 * smlnum/ulp * (||T||+||B||), smlnum ) 00246 * 00247 * 00248 CALL DCOPY( 2*N, D, 1, Y, 1 ) 00249 Y( 1 ) = DDOT( N, B, 1, X( 1+N ), 1 ) + 00250 $ SCALE*Y( 1 ) 00251 DO 50 I = 2, N 00252 Y( I ) = W*X( I+N ) + SCALE*Y( I ) 00253 50 CONTINUE 00254 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT, 00255 $ X, 1, -ONE, Y, 1 ) 00256 * 00257 Y( 1+N ) = DDOT( N, B, 1, X, 1 ) - 00258 $ SCALE*Y( 1+N ) 00259 DO 60 I = 2, N 00260 Y( I+N ) = W*X( I ) - SCALE*Y( I+N ) 00261 60 CONTINUE 00262 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT, 00263 $ X( 1+N ), 1, ONE, Y( 1+N ), 1 ) 00264 * 00265 RESID = DASUM( 2*N, Y, 1 ) 00266 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB, 00267 $ EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) ) 00268 RESID = RESID / DOMIN 00269 IF( RESID.GT.RMAX ) THEN 00270 RMAX = RESID 00271 LMAX = KNT 00272 END IF 00273 * 00274 CALL DCOPY( 2*N, D, 1, X, 1 ) 00275 KNT = KNT + 1 00276 CALL DLAQTR( .TRUE., .FALSE., N, T, LDT, B, W, 00277 $ SCALE, X, WORK, INFO ) 00278 IF( INFO.NE.0 ) 00279 $ NINFO = NINFO + 1 00280 * 00281 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| / 00282 * max(ulp*(||T||+||B||)*(||x1||+||x2||), 00283 * smlnum/ulp * (||T||+||B||), smlnum ) 00284 * 00285 CALL DCOPY( 2*N, D, 1, Y, 1 ) 00286 Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 ) 00287 DO 70 I = 2, N 00288 Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) - 00289 $ SCALE*Y( I ) 00290 70 CONTINUE 00291 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X, 00292 $ 1, ONE, Y, 1 ) 00293 * 00294 Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N ) 00295 DO 80 I = 2, N 00296 Y( I+N ) = B( I )*X( 1 ) + W*X( I ) + 00297 $ SCALE*Y( I+N ) 00298 80 CONTINUE 00299 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, 00300 $ X( 1+N ), 1, -ONE, Y( 1+N ), 1 ) 00301 * 00302 RESID = DASUM( 2*N, Y, 1 ) 00303 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB, 00304 $ EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) ) 00305 RESID = RESID / DOMIN 00306 IF( RESID.GT.RMAX ) THEN 00307 RMAX = RESID 00308 LMAX = KNT 00309 END IF 00310 * 00311 90 CONTINUE 00312 100 CONTINUE 00313 110 CONTINUE 00314 120 CONTINUE 00315 130 CONTINUE 00316 140 CONTINUE 00317 * 00318 RETURN 00319 * 00320 * End of DGET39 00321 * 00322 END
1.7.3 
