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