102 SUBROUTINE sget39( RMAX, LMAX, NINFO, KNT )
109 INTEGER KNT, LMAX, NINFO
117 parameter( ldt = 10, ldt2 = 2*ldt )
119 parameter( zero = 0.0, one = 1.0 )
122 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
124 REAL BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
125 $ SCALE, SMLNUM, W, XNORM
129 REAL SASUM, SDOT, SLAMCH, SLANGE
130 EXTERNAL isamax, sasum, sdot, slamch, slange
136 INTRINSIC abs, cos, max, real, sin, sqrt
139 INTEGER IDIM( 6 ), IVAL( 5, 5, 6 )
140 REAL 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 )
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 /
160 smlnum = slamch(
'S' )
161 bignum = one / smlnum
162 CALL slabad( smlnum, bignum )
167 vm1( 2 ) = sqrt( smlnum )
168 vm1( 3 ) = sqrt( vm1( 2 ) )
169 vm1( 4 ) = sqrt( bignum )
170 vm1( 5 ) = sqrt( vm1( 4 ) )
173 vm2( 2 ) = sqrt( smlnum )
174 vm2( 3 ) = sqrt( vm2( 2 ) )
175 vm2( 4 ) = sqrt( bignum )
176 vm2( 5 ) = sqrt( vm2( 4 ) )
179 vm3( 2 ) = sqrt( smlnum )
180 vm3( 3 ) = sqrt( vm3( 2 ) )
181 vm3( 4 ) = sqrt( bignum )
182 vm3( 5 ) = sqrt( vm3( 4 ) )
185 vm4( 2 ) = sqrt( smlnum )
186 vm4( 3 ) = sqrt( vm4( 2 ) )
187 vm4( 4 ) = sqrt( bignum )
188 vm4( 5 ) = sqrt( vm4( 4 ) )
192 vm5( 3 ) = sqrt( smlnum )
199 smlnum = smlnum / eps
213 t( i, j ) = real( ival( i, j, ndim ) )*
216 $ t( i, j ) = t( i, j )*vm5( ivm5 )
223 b( i ) = cos( real( i ) )*vm3( ivm3 )
227 d( i ) = sin( real( i ) )*vm4( ivm4 )
230 norm = slange(
'1', n, n, t, ldt, work )
231 k = isamax( n, b, 1 )
232 normtb = norm + abs( b( k ) ) + abs( w )
234 CALL scopy( n, d, 1, x, 1 )
236 CALL slaqtr( .false., .true., n, t, ldt, dum,
237 $ dumm, scale, x, work, info )
244 CALL scopy( n, d, 1, y, 1 )
245 CALL sgemv(
'No transpose', n, n, one, t, ldt,
246 $ x, 1, -scale, y, 1 )
247 xnorm = sasum( n, x, 1 )
248 resid = sasum( n, y, 1 )
249 domin = max( smlnum, ( smlnum / eps )*norm,
250 $ ( norm*eps )*xnorm )
251 resid = resid / domin
252 IF( resid.GT.rmax )
THEN
257 CALL scopy( n, d, 1, x, 1 )
259 CALL slaqtr( .true., .true., n, t, ldt, dum,
260 $ dumm, scale, x, work, info )
267 CALL scopy( n, d, 1, y, 1 )
268 CALL sgemv(
'Transpose', n, n, one, t, ldt, x,
270 xnorm = sasum( n, x, 1 )
271 resid = sasum( n, y, 1 )
272 domin = max( smlnum, ( smlnum / eps )*norm,
273 $ ( norm*eps )*xnorm )
274 resid = resid / domin
275 IF( resid.GT.rmax )
THEN
280 CALL scopy( 2*n, d, 1, x, 1 )
282 CALL slaqtr( .false., .false., n, t, ldt, b, w,
283 $ scale, x, work, info )
292 CALL scopy( 2*n, d, 1, y, 1 )
293 y( 1 ) = sdot( n, b, 1, x( 1+n ), 1 ) +
296 y( i ) = w*x( i+n ) + scale*y( i )
298 CALL sgemv(
'No transpose', n, n, one, t, ldt,
301 y( 1+n ) = sdot( n, b, 1, x, 1 ) -
304 y( i+n ) = w*x( i ) - scale*y( i+n )
306 CALL sgemv(
'No transpose', n, n, one, t, ldt,
307 $ x( 1+n ), 1, one, y( 1+n ), 1 )
309 resid = sasum( 2*n, y, 1 )
310 domin = max( smlnum, ( smlnum / eps )*normtb,
311 $ eps*( normtb*sasum( 2*n, x, 1 ) ) )
312 resid = resid / domin
313 IF( resid.GT.rmax )
THEN
318 CALL scopy( 2*n, d, 1, x, 1 )
320 CALL slaqtr( .true., .false., n, t, ldt, b, w,
321 $ scale, x, work, info )
329 CALL scopy( 2*n, d, 1, y, 1 )
330 y( 1 ) = b( 1 )*x( 1+n ) - scale*y( 1 )
332 y( i ) = b( i )*x( 1+n ) + w*x( i+n ) -
335 CALL sgemv(
'Transpose', n, n, one, t, ldt, x,
338 y( 1+n ) = b( 1 )*x( 1 ) + scale*y( 1+n )
340 y( i+n ) = b( i )*x( 1 ) + w*x( i ) +
343 CALL sgemv(
'Transpose', n, n, one, t, ldt,
344 $ x( 1+n ), 1, -one, y( 1+n ), 1 )
346 resid = sasum( 2*n, y, 1 )
347 domin = max( smlnum, ( smlnum / eps )*normtb,
348 $ eps*( normtb*sasum( 2*n, x, 1 ) ) )
349 resid = resid / domin
350 IF( resid.GT.rmax )
THEN
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slaqtr(LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
SLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine sget39(RMAX, LMAX, NINFO, KNT)
SGET39