104 SUBROUTINE dget39( RMAX, LMAX, NINFO, KNT )
112 INTEGER knt, lmax, ninfo
113 DOUBLE PRECISION rmax
120 parameter( ldt = 10, ldt2 = 2*ldt )
121 DOUBLE PRECISION zero, one
122 parameter( zero = 0.0d0, one = 1.0d0 )
125 INTEGER i, info, ivm1, ivm2, ivm3, ivm4, ivm5,
j, k, n,
127 DOUBLE PRECISION bignum, domin, dumm, eps, norm, normtb, resid,
128 $ scale, smlnum, w, xnorm
139 INTRINSIC abs, cos, dble, max, sin, sqrt
142 INTEGER idim( 6 ), ival( 5, 5, 6 )
143 DOUBLE PRECISION 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 /
164 bignum = one / smlnum
165 CALL
dlabad( 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 ) = dble( ival( i,
j, ndim ) )*
219 $ t( i,
j ) = t( i,
j )*vm5( ivm5 )
226 b( i ) = cos( dble( i ) )*vm3( ivm3 )
230 d( i ) = sin( dble( i ) )*vm4( ivm4 )
233 norm =
dlange(
'1', n, n, t, ldt, work )
235 normtb = norm + abs(
b( k ) ) + abs( w )
237 CALL
dcopy( n, d, 1,
x, 1 )
239 CALL
dlaqtr( .false., .true., n, t, ldt, dum,
240 $ dumm, scale,
x, work, info )
247 CALL
dcopy( n, d, 1, y, 1 )
248 CALL
dgemv(
'No transpose', n, n, one, t, ldt,
249 $
x, 1, -scale, y, 1 )
251 resid =
dasum( 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
dcopy( n, d, 1,
x, 1 )
262 CALL
dlaqtr( .true., .true., n, t, ldt, dum,
263 $ dumm, scale,
x, work, info )
270 CALL
dcopy( n, d, 1, y, 1 )
271 CALL
dgemv(
'Transpose', n, n, one, t, ldt,
x,
274 resid =
dasum( 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
dcopy( 2*n, d, 1,
x, 1 )
285 CALL
dlaqtr( .false., .false., n, t, ldt,
b, w,
286 $ scale,
x, work, info )
295 CALL
dcopy( 2*n, d, 1, y, 1 )
296 y( 1 ) =
ddot( n,
b, 1,
x( 1+n ), 1 ) +
299 y( i ) = w*
x( i+n ) + scale*y( i )
301 CALL
dgemv(
'No transpose', n, n, one, t, ldt,
304 y( 1+n ) =
ddot( n,
b, 1,
x, 1 ) -
307 y( i+n ) = w*
x( i ) - scale*y( i+n )
309 CALL
dgemv(
'No transpose', n, n, one, t, ldt,
310 $
x( 1+n ), 1, one, y( 1+n ), 1 )
312 resid =
dasum( 2*n, y, 1 )
313 domin = max( smlnum, ( smlnum / eps )*normtb,
314 $ eps*( normtb*
dasum( 2*n,
x, 1 ) ) )
315 resid = resid / domin
316 IF( resid.GT.rmax )
THEN
321 CALL
dcopy( 2*n, d, 1,
x, 1 )
323 CALL
dlaqtr( .true., .false., n, t, ldt,
b, w,
324 $ scale,
x, work, info )
332 CALL
dcopy( 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
dgemv(
'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
dgemv(
'Transpose', n, n, one, t, ldt,
347 $
x( 1+n ), 1, -one, y( 1+n ), 1 )
349 resid =
dasum( 2*n, y, 1 )
350 domin = max( smlnum, ( smlnum / eps )*normtb,
351 $ eps*( normtb*
dasum( 2*n,
x, 1 ) ) )
352 resid = resid / domin
353 IF( resid.GT.rmax )
THEN
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dget39(RMAX, LMAX, NINFO, KNT)
DGET39
subroutine dlabad(SMALL, LARGE)
DLABAD
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
DOUBLE PRECISION function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dlaqtr(LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
DOUBLE PRECISION function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
DOUBLE PRECISION function dasum(N, DX, INCX)
DASUM
INTEGER function idamax(N, DX, INCX)
IDAMAX