1*> \brief \b SGET39 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 SGET39( RMAX, LMAX, NINFO, KNT ) 12* 13* .. Scalar Arguments .. 14* INTEGER KNT, LMAX, NINFO 15* REAL RMAX 16* .. 17* 18* 19*> \par Purpose: 20* ============= 21*> 22*> \verbatim 23*> 24*> SGET39 tests SLAQTR, 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 STRSNA. 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 REAL 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*> \date November 2011 100* 101*> \ingroup single_eig 102* 103* ===================================================================== 104 SUBROUTINE SGET39( RMAX, LMAX, NINFO, KNT ) 105* 106* -- LAPACK test routine (version 3.4.0) -- 107* -- LAPACK is a software package provided by Univ. of Tennessee, -- 108* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 109* November 2011 110* 111* .. Scalar Arguments .. 112 INTEGER KNT, LMAX, NINFO 113 REAL RMAX 114* .. 115* 116* ===================================================================== 117* 118* .. Parameters .. 119 INTEGER LDT, LDT2 120 PARAMETER ( LDT = 10, LDT2 = 2*LDT ) 121 REAL ZERO, ONE 122 PARAMETER ( ZERO = 0.0, ONE = 1.0 ) 123* .. 124* .. Local Scalars .. 125 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N, 126 $ NDIM 127 REAL BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID, 128 $ SCALE, SMLNUM, W, XNORM 129* .. 130* .. External Functions .. 131 INTEGER ISAMAX 132 REAL SASUM, SDOT, SLAMCH, SLANGE 133 EXTERNAL ISAMAX, SASUM, SDOT, SLAMCH, SLANGE 134* .. 135* .. External Subroutines .. 136 EXTERNAL SCOPY, SGEMV, SLABAD, SLAQTR 137* .. 138* .. Intrinsic Functions .. 139 INTRINSIC ABS, COS, MAX, REAL, SIN, SQRT 140* .. 141* .. Local Arrays .. 142 INTEGER IDIM( 6 ), IVAL( 5, 5, 6 ) 143 REAL B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ), 144 $ VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ), 145 $ VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 ) 146* .. 147* .. Data statements .. 148 DATA IDIM / 4, 5*5 / 149 DATA IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0, 150 $ 4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4, 151 $ 0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2, 152 $ 0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1, 153 $ 4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2, 154 $ -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0, 155 $ 5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1, 156 $ 4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 / 157* .. 158* .. Executable Statements .. 159* 160* Get machine parameters 161* 162 EPS = SLAMCH( 'P' ) 163 SMLNUM = SLAMCH( 'S' ) 164 BIGNUM = ONE / SMLNUM 165 CALL SLABAD( SMLNUM, BIGNUM ) 166* 167* Set up test case parameters 168* 169 VM1( 1 ) = ONE 170 VM1( 2 ) = SQRT( SMLNUM ) 171 VM1( 3 ) = SQRT( VM1( 2 ) ) 172 VM1( 4 ) = SQRT( BIGNUM ) 173 VM1( 5 ) = SQRT( VM1( 4 ) ) 174* 175 VM2( 1 ) = ONE 176 VM2( 2 ) = SQRT( SMLNUM ) 177 VM2( 3 ) = SQRT( VM2( 2 ) ) 178 VM2( 4 ) = SQRT( BIGNUM ) 179 VM2( 5 ) = SQRT( VM2( 4 ) ) 180* 181 VM3( 1 ) = ONE 182 VM3( 2 ) = SQRT( SMLNUM ) 183 VM3( 3 ) = SQRT( VM3( 2 ) ) 184 VM3( 4 ) = SQRT( BIGNUM ) 185 VM3( 5 ) = SQRT( VM3( 4 ) ) 186* 187 VM4( 1 ) = ONE 188 VM4( 2 ) = SQRT( SMLNUM ) 189 VM4( 3 ) = SQRT( VM4( 2 ) ) 190 VM4( 4 ) = SQRT( BIGNUM ) 191 VM4( 5 ) = SQRT( VM4( 4 ) ) 192* 193 VM5( 1 ) = ONE 194 VM5( 2 ) = EPS 195 VM5( 3 ) = SQRT( SMLNUM ) 196* 197* Initalization 198* 199 KNT = 0 200 RMAX = ZERO 201 NINFO = 0 202 SMLNUM = SMLNUM / EPS 203* 204* Begin test loop 205* 206 DO 140 IVM5 = 1, 3 207 DO 130 IVM4 = 1, 5 208 DO 120 IVM3 = 1, 5 209 DO 110 IVM2 = 1, 5 210 DO 100 IVM1 = 1, 5 211 DO 90 NDIM = 1, 6 212* 213 N = IDIM( NDIM ) 214 DO 20 I = 1, N 215 DO 10 J = 1, N 216 T( I, J ) = REAL( IVAL( I, J, NDIM ) )* 217 $ VM1( IVM1 ) 218 IF( I.GE.J ) 219 $ T( I, J ) = T( I, J )*VM5( IVM5 ) 220 10 CONTINUE 221 20 CONTINUE 222* 223 W = ONE*VM2( IVM2 ) 224* 225 DO 30 I = 1, N 226 B( I ) = COS( REAL( I ) )*VM3( IVM3 ) 227 30 CONTINUE 228* 229 DO 40 I = 1, 2*N 230 D( I ) = SIN( REAL( I ) )*VM4( IVM4 ) 231 40 CONTINUE 232* 233 NORM = SLANGE( '1', N, N, T, LDT, WORK ) 234 K = ISAMAX( N, B, 1 ) 235 NORMTB = NORM + ABS( B( K ) ) + ABS( W ) 236* 237 CALL SCOPY( N, D, 1, X, 1 ) 238 KNT = KNT + 1 239 CALL SLAQTR( .FALSE., .TRUE., N, T, LDT, DUM, 240 $ DUMM, SCALE, X, WORK, INFO ) 241 IF( INFO.NE.0 ) 242 $ NINFO = NINFO + 1 243* 244* || T*x - scale*d || / 245* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum) 246* 247 CALL SCOPY( N, D, 1, Y, 1 ) 248 CALL SGEMV( 'No transpose', N, N, ONE, T, LDT, 249 $ X, 1, -SCALE, Y, 1 ) 250 XNORM = SASUM( N, X, 1 ) 251 RESID = SASUM( N, Y, 1 ) 252 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM, 253 $ ( NORM*EPS )*XNORM ) 254 RESID = RESID / DOMIN 255 IF( RESID.GT.RMAX ) THEN 256 RMAX = RESID 257 LMAX = KNT 258 END IF 259* 260 CALL SCOPY( N, D, 1, X, 1 ) 261 KNT = KNT + 1 262 CALL SLAQTR( .TRUE., .TRUE., N, T, LDT, DUM, 263 $ DUMM, SCALE, X, WORK, INFO ) 264 IF( INFO.NE.0 ) 265 $ NINFO = NINFO + 1 266* 267* || T*x - scale*d || / 268* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum) 269* 270 CALL SCOPY( N, D, 1, Y, 1 ) 271 CALL SGEMV( 'Transpose', N, N, ONE, T, LDT, X, 272 $ 1, -SCALE, Y, 1 ) 273 XNORM = SASUM( N, X, 1 ) 274 RESID = SASUM( N, Y, 1 ) 275 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM, 276 $ ( NORM*EPS )*XNORM ) 277 RESID = RESID / DOMIN 278 IF( RESID.GT.RMAX ) THEN 279 RMAX = RESID 280 LMAX = KNT 281 END IF 282* 283 CALL SCOPY( 2*N, D, 1, X, 1 ) 284 KNT = KNT + 1 285 CALL SLAQTR( .FALSE., .FALSE., N, T, LDT, B, W, 286 $ SCALE, X, WORK, INFO ) 287 IF( INFO.NE.0 ) 288 $ NINFO = NINFO + 1 289* 290* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| / 291* max(ulp*(||T||+||B||)*(||x1||+||x2||), 292* smlnum/ulp * (||T||+||B||), smlnum ) 293* 294* 295 CALL SCOPY( 2*N, D, 1, Y, 1 ) 296 Y( 1 ) = SDOT( N, B, 1, X( 1+N ), 1 ) + 297 $ SCALE*Y( 1 ) 298 DO 50 I = 2, N 299 Y( I ) = W*X( I+N ) + SCALE*Y( I ) 300 50 CONTINUE 301 CALL SGEMV( 'No transpose', N, N, ONE, T, LDT, 302 $ X, 1, -ONE, Y, 1 ) 303* 304 Y( 1+N ) = SDOT( N, B, 1, X, 1 ) - 305 $ SCALE*Y( 1+N ) 306 DO 60 I = 2, N 307 Y( I+N ) = W*X( I ) - SCALE*Y( I+N ) 308 60 CONTINUE 309 CALL SGEMV( 'No transpose', N, N, ONE, T, LDT, 310 $ X( 1+N ), 1, ONE, Y( 1+N ), 1 ) 311* 312 RESID = SASUM( 2*N, Y, 1 ) 313 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB, 314 $ EPS*( NORMTB*SASUM( 2*N, X, 1 ) ) ) 315 RESID = RESID / DOMIN 316 IF( RESID.GT.RMAX ) THEN 317 RMAX = RESID 318 LMAX = KNT 319 END IF 320* 321 CALL SCOPY( 2*N, D, 1, X, 1 ) 322 KNT = KNT + 1 323 CALL SLAQTR( .TRUE., .FALSE., N, T, LDT, B, W, 324 $ SCALE, X, WORK, INFO ) 325 IF( INFO.NE.0 ) 326 $ NINFO = NINFO + 1 327* 328* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| / 329* max(ulp*(||T||+||B||)*(||x1||+||x2||), 330* smlnum/ulp * (||T||+||B||), smlnum ) 331* 332 CALL SCOPY( 2*N, D, 1, Y, 1 ) 333 Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 ) 334 DO 70 I = 2, N 335 Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) - 336 $ SCALE*Y( I ) 337 70 CONTINUE 338 CALL SGEMV( 'Transpose', N, N, ONE, T, LDT, X, 339 $ 1, ONE, Y, 1 ) 340* 341 Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N ) 342 DO 80 I = 2, N 343 Y( I+N ) = B( I )*X( 1 ) + W*X( I ) + 344 $ SCALE*Y( I+N ) 345 80 CONTINUE 346 CALL SGEMV( 'Transpose', N, N, ONE, T, LDT, 347 $ X( 1+N ), 1, -ONE, Y( 1+N ), 1 ) 348* 349 RESID = SASUM( 2*N, Y, 1 ) 350 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB, 351 $ EPS*( NORMTB*SASUM( 2*N, X, 1 ) ) ) 352 RESID = RESID / DOMIN 353 IF( RESID.GT.RMAX ) THEN 354 RMAX = RESID 355 LMAX = KNT 356 END IF 357* 358 90 CONTINUE 359 100 CONTINUE 360 110 CONTINUE 361 120 CONTINUE 362 130 CONTINUE 363 140 CONTINUE 364* 365 RETURN 366* 367* End of SGET39 368* 369 END 370