104 SUBROUTINE sget39( RMAX, LMAX, NINFO, KNT )
112 INTEGER KNT, LMAX, NINFO
120 parameter ( ldt = 10, ldt2 = 2*ldt )
122 parameter ( zero = 0.0, one = 1.0 )
125 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
127 REAL BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
128 $ scale, smlnum, w, xnorm
132 REAL SASUM, SDOT, SLAMCH, SLANGE
133 EXTERNAL isamax, sasum, sdot, slamch, slange
139 INTRINSIC abs, cos, max,
REAL, SIN, SQRT
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 )
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 /
163 smlnum = slamch(
'S' )
164 bignum = one / smlnum
165 CALL slabad( smlnum, bignum )
170 vm1( 2 ) = sqrt( smlnum )
171 vm1( 3 ) = sqrt( vm1( 2 ) )
172 vm1( 4 ) = sqrt( bignum )
173 vm1( 5 ) = sqrt( vm1( 4 ) )
176 vm2( 2 ) = sqrt( smlnum )
177 vm2( 3 ) = sqrt( vm2( 2 ) )
178 vm2( 4 ) = sqrt( bignum )
179 vm2( 5 ) = sqrt( vm2( 4 ) )
182 vm3( 2 ) = sqrt( smlnum )
183 vm3( 3 ) = sqrt( vm3( 2 ) )
184 vm3( 4 ) = sqrt( bignum )
185 vm3( 5 ) = sqrt( vm3( 4 ) )
188 vm4( 2 ) = sqrt( smlnum )
189 vm4( 3 ) = sqrt( vm4( 2 ) )
190 vm4( 4 ) = sqrt( bignum )
191 vm4( 5 ) = sqrt( vm4( 4 ) )
195 vm5( 3 ) = sqrt( smlnum )
202 smlnum = smlnum / eps
216 t( i, j ) =
REAL( IVAL( I, J, NDIM ) )*
219 $ t( i, j ) = t( i, j )*vm5( ivm5 )
226 b( i ) = cos(
REAL( I ) )*vm3( IVM3 )
230 d( i ) = sin(
REAL( I ) )*vm4( IVM4 )
233 norm = slange(
'1', n, n, t, ldt, work )
234 k = isamax( n, b, 1 )
235 normtb = norm + abs( b( k ) ) + abs( w )
237 CALL scopy( n, d, 1, x, 1 )
239 CALL slaqtr( .false., .true., n, t, ldt, dum,
240 $ dumm, scale, x, work, info )
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
260 CALL scopy( n, d, 1, x, 1 )
262 CALL slaqtr( .true., .true., n, t, ldt, dum,
263 $ dumm, scale, x, work, info )
270 CALL scopy( n, d, 1, y, 1 )
271 CALL sgemv(
'Transpose', n, n, one, t, ldt, x,
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
283 CALL scopy( 2*n, d, 1, x, 1 )
285 CALL slaqtr( .false., .false., n, t, ldt, b, w,
286 $ scale, x, work, info )
295 CALL scopy( 2*n, d, 1, y, 1 )
296 y( 1 ) = sdot( n, b, 1, x( 1+n ), 1 ) +
299 y( i ) = w*x( i+n ) + scale*y( i )
301 CALL sgemv(
'No transpose', n, n, one, t, ldt,
304 y( 1+n ) = sdot( n, b, 1, x, 1 ) -
307 y( i+n ) = w*x( i ) - scale*y( i+n )
309 CALL sgemv(
'No transpose', n, n, one, t, ldt,
310 $ x( 1+n ), 1, one, y( 1+n ), 1 )
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
321 CALL scopy( 2*n, d, 1, x, 1 )
323 CALL slaqtr( .true., .false., n, t, ldt, b, w,
324 $ scale, x, work, info )
332 CALL scopy( 2*n, d, 1, y, 1 )
333 y( 1 ) = b( 1 )*x( 1+n ) - scale*y( 1 )
335 y( i ) = b( i )*x( 1+n ) + w*x( i+n ) -
338 CALL sgemv(
'Transpose', n, n, one, t, ldt, x,
341 y( 1+n ) = b( 1 )*x( 1 ) + scale*y( 1+n )
343 y( i+n ) = b( i )*x( 1 ) + w*x( i ) +
346 CALL sgemv(
'Transpose', n, n, one, t, ldt,
347 $ x( 1+n ), 1, -one, y( 1+n ), 1 )
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
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
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 sget39(RMAX, LMAX, NINFO, KNT)
SGET39
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY