1*> \brief \b DGET39 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT ) 12* 13* .. Scalar Arguments .. 14* INTEGER KNT, LMAX, NINFO 15* DOUBLE PRECISION RMAX 16* .. 17* 18* 19*> \par Purpose: 20* ============= 21*> 22*> \verbatim 23*> 24*> DGET39 tests DLAQTR, a routine for solving the real or 25*> special complex quasi upper triangular system 26*> 27*> op(T)*p = scale*c, 28*> or 29*> op(T + iB)*(p+iq) = scale*(c+id), 30*> 31*> in real arithmetic. T is upper quasi-triangular. 32*> If it is complex, then the first diagonal block of T must be 33*> 1 by 1, B has the special structure 34*> 35*> B = [ b(1) b(2) ... b(n) ] 36*> [ w ] 37*> [ w ] 38*> [ . ] 39*> [ w ] 40*> 41*> op(A) = A or A', where A' denotes the conjugate transpose of 42*> the matrix A. 43*> 44*> On input, X = [ c ]. On output, X = [ p ]. 45*> [ d ] [ q ] 46*> 47*> Scale is an output less than or equal to 1, chosen to avoid 48*> overflow in X. 49*> This subroutine is specially designed for the condition number 50*> estimation in the eigenproblem routine DTRSNA. 51*> 52*> The test code verifies that the following residual is order 1: 53*> 54*> ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| 55*> ----------------------------------------- 56*> max(ulp*(||T||+||B||)*(||x1||+||x2||), 57*> (||T||+||B||)*smlnum/ulp, 58*> smlnum) 59*> 60*> (The (||T||+||B||)*smlnum/ulp term accounts for possible 61*> (gradual or nongradual) underflow in x1 and x2.) 62*> \endverbatim 63* 64* Arguments: 65* ========== 66* 67*> \param[out] RMAX 68*> \verbatim 69*> RMAX is DOUBLE PRECISION 70*> Value of the largest test ratio. 71*> \endverbatim 72*> 73*> \param[out] LMAX 74*> \verbatim 75*> LMAX is INTEGER 76*> Example number where largest test ratio achieved. 77*> \endverbatim 78*> 79*> \param[out] NINFO 80*> \verbatim 81*> NINFO is INTEGER 82*> Number of examples where INFO is nonzero. 83*> \endverbatim 84*> 85*> \param[out] KNT 86*> \verbatim 87*> KNT is INTEGER 88*> Total number of examples tested. 89*> \endverbatim 90* 91* Authors: 92* ======== 93* 94*> \author Univ. of Tennessee 95*> \author Univ. of California Berkeley 96*> \author Univ. of Colorado Denver 97*> \author NAG Ltd. 98* 99*> \ingroup double_eig 100* 101* ===================================================================== 102 SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT ) 103* 104* -- LAPACK test routine -- 105* -- LAPACK is a software package provided by Univ. of Tennessee, -- 106* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 107* 108* .. Scalar Arguments .. 109 INTEGER KNT, LMAX, NINFO 110 DOUBLE PRECISION RMAX 111* .. 112* 113* ===================================================================== 114* 115* .. Parameters .. 116 INTEGER LDT, LDT2 117 PARAMETER ( LDT = 10, LDT2 = 2*LDT ) 118 DOUBLE PRECISION ZERO, ONE 119 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 120* .. 121* .. Local Scalars .. 122 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N, 123 $ NDIM 124 DOUBLE PRECISION BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID, 125 $ SCALE, SMLNUM, W, XNORM 126* .. 127* .. External Functions .. 128 INTEGER IDAMAX 129 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE 130 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE 131* .. 132* .. External Subroutines .. 133 EXTERNAL DCOPY, DGEMV, DLABAD, DLAQTR 134* .. 135* .. Intrinsic Functions .. 136 INTRINSIC ABS, COS, DBLE, MAX, SIN, SQRT 137* .. 138* .. Local Arrays .. 139 INTEGER IDIM( 6 ), IVAL( 5, 5, 6 ) 140 DOUBLE PRECISION B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ), 141 $ VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ), 142 $ VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 ) 143* .. 144* .. Data statements .. 145 DATA IDIM / 4, 5*5 / 146 DATA IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0, 147 $ 4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4, 148 $ 0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2, 149 $ 0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1, 150 $ 4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2, 151 $ -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0, 152 $ 5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1, 153 $ 4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 / 154* .. 155* .. Executable Statements .. 156* 157* Get machine parameters 158* 159 EPS = DLAMCH( 'P' ) 160 SMLNUM = DLAMCH( 'S' ) 161 BIGNUM = ONE / SMLNUM 162 CALL DLABAD( SMLNUM, BIGNUM ) 163* 164* Set up test case parameters 165* 166 VM1( 1 ) = ONE 167 VM1( 2 ) = SQRT( SMLNUM ) 168 VM1( 3 ) = SQRT( VM1( 2 ) ) 169 VM1( 4 ) = SQRT( BIGNUM ) 170 VM1( 5 ) = SQRT( VM1( 4 ) ) 171* 172 VM2( 1 ) = ONE 173 VM2( 2 ) = SQRT( SMLNUM ) 174 VM2( 3 ) = SQRT( VM2( 2 ) ) 175 VM2( 4 ) = SQRT( BIGNUM ) 176 VM2( 5 ) = SQRT( VM2( 4 ) ) 177* 178 VM3( 1 ) = ONE 179 VM3( 2 ) = SQRT( SMLNUM ) 180 VM3( 3 ) = SQRT( VM3( 2 ) ) 181 VM3( 4 ) = SQRT( BIGNUM ) 182 VM3( 5 ) = SQRT( VM3( 4 ) ) 183* 184 VM4( 1 ) = ONE 185 VM4( 2 ) = SQRT( SMLNUM ) 186 VM4( 3 ) = SQRT( VM4( 2 ) ) 187 VM4( 4 ) = SQRT( BIGNUM ) 188 VM4( 5 ) = SQRT( VM4( 4 ) ) 189* 190 VM5( 1 ) = ONE 191 VM5( 2 ) = EPS 192 VM5( 3 ) = SQRT( SMLNUM ) 193* 194* Initialization 195* 196 KNT = 0 197 RMAX = ZERO 198 NINFO = 0 199 SMLNUM = SMLNUM / EPS 200* 201* Begin test loop 202* 203 DO 140 IVM5 = 1, 3 204 DO 130 IVM4 = 1, 5 205 DO 120 IVM3 = 1, 5 206 DO 110 IVM2 = 1, 5 207 DO 100 IVM1 = 1, 5 208 DO 90 NDIM = 1, 6 209* 210 N = IDIM( NDIM ) 211 DO 20 I = 1, N 212 DO 10 J = 1, N 213 T( I, J ) = DBLE( IVAL( I, J, NDIM ) )* 214 $ VM1( IVM1 ) 215 IF( I.GE.J ) 216 $ T( I, J ) = T( I, J )*VM5( IVM5 ) 217 10 CONTINUE 218 20 CONTINUE 219* 220 W = ONE*VM2( IVM2 ) 221* 222 DO 30 I = 1, N 223 B( I ) = COS( DBLE( I ) )*VM3( IVM3 ) 224 30 CONTINUE 225* 226 DO 40 I = 1, 2*N 227 D( I ) = SIN( DBLE( I ) )*VM4( IVM4 ) 228 40 CONTINUE 229* 230 NORM = DLANGE( '1', N, N, T, LDT, WORK ) 231 K = IDAMAX( N, B, 1 ) 232 NORMTB = NORM + ABS( B( K ) ) + ABS( W ) 233* 234 CALL DCOPY( N, D, 1, X, 1 ) 235 KNT = KNT + 1 236 CALL DLAQTR( .FALSE., .TRUE., N, T, LDT, DUM, 237 $ DUMM, SCALE, X, WORK, INFO ) 238 IF( INFO.NE.0 ) 239 $ NINFO = NINFO + 1 240* 241* || T*x - scale*d || / 242* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum) 243* 244 CALL DCOPY( N, D, 1, Y, 1 ) 245 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT, 246 $ X, 1, -SCALE, Y, 1 ) 247 XNORM = DASUM( N, X, 1 ) 248 RESID = DASUM( N, Y, 1 ) 249 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM, 250 $ ( NORM*EPS )*XNORM ) 251 RESID = RESID / DOMIN 252 IF( RESID.GT.RMAX ) THEN 253 RMAX = RESID 254 LMAX = KNT 255 END IF 256* 257 CALL DCOPY( N, D, 1, X, 1 ) 258 KNT = KNT + 1 259 CALL DLAQTR( .TRUE., .TRUE., N, T, LDT, DUM, 260 $ DUMM, SCALE, X, WORK, INFO ) 261 IF( INFO.NE.0 ) 262 $ NINFO = NINFO + 1 263* 264* || T*x - scale*d || / 265* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum) 266* 267 CALL DCOPY( N, D, 1, Y, 1 ) 268 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X, 269 $ 1, -SCALE, Y, 1 ) 270 XNORM = DASUM( N, X, 1 ) 271 RESID = DASUM( N, Y, 1 ) 272 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM, 273 $ ( NORM*EPS )*XNORM ) 274 RESID = RESID / DOMIN 275 IF( RESID.GT.RMAX ) THEN 276 RMAX = RESID 277 LMAX = KNT 278 END IF 279* 280 CALL DCOPY( 2*N, D, 1, X, 1 ) 281 KNT = KNT + 1 282 CALL DLAQTR( .FALSE., .FALSE., N, T, LDT, B, W, 283 $ SCALE, X, WORK, INFO ) 284 IF( INFO.NE.0 ) 285 $ NINFO = NINFO + 1 286* 287* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| / 288* max(ulp*(||T||+||B||)*(||x1||+||x2||), 289* smlnum/ulp * (||T||+||B||), smlnum ) 290* 291* 292 CALL DCOPY( 2*N, D, 1, Y, 1 ) 293 Y( 1 ) = DDOT( N, B, 1, X( 1+N ), 1 ) + 294 $ SCALE*Y( 1 ) 295 DO 50 I = 2, N 296 Y( I ) = W*X( I+N ) + SCALE*Y( I ) 297 50 CONTINUE 298 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT, 299 $ X, 1, -ONE, Y, 1 ) 300* 301 Y( 1+N ) = DDOT( N, B, 1, X, 1 ) - 302 $ SCALE*Y( 1+N ) 303 DO 60 I = 2, N 304 Y( I+N ) = W*X( I ) - SCALE*Y( I+N ) 305 60 CONTINUE 306 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT, 307 $ X( 1+N ), 1, ONE, Y( 1+N ), 1 ) 308* 309 RESID = DASUM( 2*N, Y, 1 ) 310 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB, 311 $ EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) ) 312 RESID = RESID / DOMIN 313 IF( RESID.GT.RMAX ) THEN 314 RMAX = RESID 315 LMAX = KNT 316 END IF 317* 318 CALL DCOPY( 2*N, D, 1, X, 1 ) 319 KNT = KNT + 1 320 CALL DLAQTR( .TRUE., .FALSE., N, T, LDT, B, W, 321 $ SCALE, X, WORK, INFO ) 322 IF( INFO.NE.0 ) 323 $ NINFO = NINFO + 1 324* 325* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| / 326* max(ulp*(||T||+||B||)*(||x1||+||x2||), 327* smlnum/ulp * (||T||+||B||), smlnum ) 328* 329 CALL DCOPY( 2*N, D, 1, Y, 1 ) 330 Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 ) 331 DO 70 I = 2, N 332 Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) - 333 $ SCALE*Y( I ) 334 70 CONTINUE 335 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X, 336 $ 1, ONE, Y, 1 ) 337* 338 Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N ) 339 DO 80 I = 2, N 340 Y( I+N ) = B( I )*X( 1 ) + W*X( I ) + 341 $ SCALE*Y( I+N ) 342 80 CONTINUE 343 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, 344 $ X( 1+N ), 1, -ONE, Y( 1+N ), 1 ) 345* 346 RESID = DASUM( 2*N, Y, 1 ) 347 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB, 348 $ EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) ) 349 RESID = RESID / DOMIN 350 IF( RESID.GT.RMAX ) THEN 351 RMAX = RESID 352 LMAX = KNT 353 END IF 354* 355 90 CONTINUE 356 100 CONTINUE 357 110 CONTINUE 358 120 CONTINUE 359 130 CONTINUE 360 140 CONTINUE 361* 362 RETURN 363* 364* End of DGET39 365* 366 END 367