82
83
84
85
86
87
88 INTEGER KNT, LMAX, NINFO
89 DOUBLE PRECISION RMAX
90
91
92
93
94
95 DOUBLE PRECISION ZERO, ONE
96 parameter( zero = 0.0d0, one = 1.0d0 )
97 DOUBLE PRECISION TWO, FOUR, EIGHT
98 parameter( two = 2.0d0, four = 4.0d0, eight = 8.0d0 )
99
100
101 LOGICAL LTRANL, LTRANR
102 INTEGER IB, IB1, IB2, IB3, INFO, ISGN, ITL, ITLSCL,
103 $ ITR, ITRANL, ITRANR, ITRSCL, N1, N2
104 DOUBLE PRECISION BIGNUM, DEN, EPS, RES, SCALE, SGN, SMLNUM, TMP,
105 $ TNRM, XNORM, XNRM
106
107
108 INTEGER ITVAL( 2, 2, 8 )
109 DOUBLE PRECISION B( 2, 2 ), TL( 2, 2 ), TR( 2, 2 ), VAL( 3 ),
110 $ X( 2, 2 )
111
112
113 DOUBLE PRECISION DLAMCH
115
116
118
119
120 INTRINSIC abs, max, min, sqrt
121
122
123 DATA itval / 8, 4, 2, 1, 4, 8, 1, 2, 2, 1, 8, 4, 1,
124 $ 2, 4, 8, 9, 4, 2, 1, 4, 9, 1, 2, 2, 1, 9, 4, 1,
125 $ 2, 4, 9 /
126
127
128
129
130
132 smlnum =
dlamch(
'S' ) / eps
133 bignum = one / smlnum
134
135
136
137 val( 1 ) = sqrt( smlnum )
138 val( 2 ) = one
139 val( 3 ) = sqrt( bignum )
140
141 knt = 0
142 ninfo = 0
143 lmax = 0
144 rmax = zero
145
146
147
148 DO 230 itranl = 0, 1
149 DO 220 itranr = 0, 1
150 DO 210 isgn = -1, 1, 2
151 sgn = isgn
152 ltranl = itranl.EQ.1
153 ltranr = itranr.EQ.1
154
155 n1 = 1
156 n2 = 1
157 DO 30 itl = 1, 3
158 DO 20 itr = 1, 3
159 DO 10 ib = 1, 3
160 tl( 1, 1 ) = val( itl )
161 tr( 1, 1 ) = val( itr )
162 b( 1, 1 ) = val( ib )
163 knt = knt + 1
164 CALL dlasy2( ltranl, ltranr, isgn, n1, n2, tl,
165 $ 2, tr, 2, b, 2, scale, x, 2, xnorm,
166 $ info )
167 IF( info.NE.0 )
168 $ ninfo = ninfo + 1
169 res = abs( ( tl( 1, 1 )+sgn*tr( 1, 1 ) )*
170 $ x( 1, 1 )-scale*b( 1, 1 ) )
171 IF( info.EQ.0 ) THEN
172 den = max( eps*( ( abs( tr( 1,
173 $ 1 ) )+abs( tl( 1, 1 ) ) )*abs( x( 1,
174 $ 1 ) ) ), smlnum )
175 ELSE
176 den = smlnum*max( abs( x( 1, 1 ) ), one )
177 END IF
178 res = res / den
179 IF( scale.GT.one )
180 $ res = res + one / eps
181 res = res + abs( xnorm-abs( x( 1, 1 ) ) ) /
182 $ max( smlnum, xnorm ) / eps
183 IF( info.NE.0 .AND. info.NE.1 )
184 $ res = res + one / eps
185 IF( res.GT.rmax ) THEN
186 lmax = knt
187 rmax = res
188 END IF
189 10 CONTINUE
190 20 CONTINUE
191 30 CONTINUE
192
193 n1 = 2
194 n2 = 1
195 DO 80 itl = 1, 8
196 DO 70 itlscl = 1, 3
197 DO 60 itr = 1, 3
198 DO 50 ib1 = 1, 3
199 DO 40 ib2 = 1, 3
200 b( 1, 1 ) = val( ib1 )
201 b( 2, 1 ) = -four*val( ib2 )
202 tl( 1, 1 ) = itval( 1, 1, itl )*
203 $ val( itlscl )
204 tl( 2, 1 ) = itval( 2, 1, itl )*
205 $ val( itlscl )
206 tl( 1, 2 ) = itval( 1, 2, itl )*
207 $ val( itlscl )
208 tl( 2, 2 ) = itval( 2, 2, itl )*
209 $ val( itlscl )
210 tr( 1, 1 ) = val( itr )
211 knt = knt + 1
212 CALL dlasy2( ltranl, ltranr, isgn, n1, n2,
213 $ tl, 2, tr, 2, b, 2, scale, x,
214 $ 2, xnorm, info )
215 IF( info.NE.0 )
216 $ ninfo = ninfo + 1
217 IF( ltranl ) THEN
218 tmp = tl( 1, 2 )
219 tl( 1, 2 ) = tl( 2, 1 )
220 tl( 2, 1 ) = tmp
221 END IF
222 res = abs( ( tl( 1, 1 )+sgn*tr( 1, 1 ) )*
223 $ x( 1, 1 )+tl( 1, 2 )*x( 2, 1 )-
224 $ scale*b( 1, 1 ) )
225 res = res + abs( ( tl( 2, 2 )+sgn*tr( 1,
226 $ 1 ) )*x( 2, 1 )+tl( 2, 1 )*
227 $ x( 1, 1 )-scale*b( 2, 1 ) )
228 tnrm = abs( tr( 1, 1 ) ) +
229 $ abs( tl( 1, 1 ) ) +
230 $ abs( tl( 1, 2 ) ) +
231 $ abs( tl( 2, 1 ) ) +
232 $ abs( tl( 2, 2 ) )
233 xnrm = max( abs( x( 1, 1 ) ),
234 $ abs( x( 2, 1 ) ) )
235 den = max( smlnum, smlnum*xnrm,
236 $ ( tnrm*eps )*xnrm )
237 res = res / den
238 IF( scale.GT.one )
239 $ res = res + one / eps
240 res = res + abs( xnorm-xnrm ) /
241 $ max( smlnum, xnorm ) / eps
242 IF( res.GT.rmax ) THEN
243 lmax = knt
244 rmax = res
245 END IF
246 40 CONTINUE
247 50 CONTINUE
248 60 CONTINUE
249 70 CONTINUE
250 80 CONTINUE
251
252 n1 = 1
253 n2 = 2
254 DO 130 itr = 1, 8
255 DO 120 itrscl = 1, 3
256 DO 110 itl = 1, 3
257 DO 100 ib1 = 1, 3
258 DO 90 ib2 = 1, 3
259 b( 1, 1 ) = val( ib1 )
260 b( 1, 2 ) = -two*val( ib2 )
261 tr( 1, 1 ) = itval( 1, 1, itr )*
262 $ val( itrscl )
263 tr( 2, 1 ) = itval( 2, 1, itr )*
264 $ val( itrscl )
265 tr( 1, 2 ) = itval( 1, 2, itr )*
266 $ val( itrscl )
267 tr( 2, 2 ) = itval( 2, 2, itr )*
268 $ val( itrscl )
269 tl( 1, 1 ) = val( itl )
270 knt = knt + 1
271 CALL dlasy2( ltranl, ltranr, isgn, n1, n2,
272 $ tl, 2, tr, 2, b, 2, scale, x,
273 $ 2, xnorm, info )
274 IF( info.NE.0 )
275 $ ninfo = ninfo + 1
276 IF( ltranr ) THEN
277 tmp = tr( 1, 2 )
278 tr( 1, 2 ) = tr( 2, 1 )
279 tr( 2, 1 ) = tmp
280 END IF
281 tnrm = abs( tl( 1, 1 ) ) +
282 $ abs( tr( 1, 1 ) ) +
283 $ abs( tr( 1, 2 ) ) +
284 $ abs( tr( 2, 2 ) ) +
285 $ abs( tr( 2, 1 ) )
286 xnrm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
287 res = abs( ( ( tl( 1, 1 )+sgn*tr( 1,
288 $ 1 ) ) )*( x( 1, 1 ) )+
289 $ ( sgn*tr( 2, 1 ) )*( x( 1, 2 ) )-
290 $ ( scale*b( 1, 1 ) ) )
291 res = res + abs( ( ( tl( 1, 1 )+sgn*tr( 2,
292 $ 2 ) ) )*( x( 1, 2 ) )+
293 $ ( sgn*tr( 1, 2 ) )*( x( 1, 1 ) )-
294 $ ( scale*b( 1, 2 ) ) )
295 den = max( smlnum, smlnum*xnrm,
296 $ ( tnrm*eps )*xnrm )
297 res = res / den
298 IF( scale.GT.one )
299 $ res = res + one / eps
300 res = res + abs( xnorm-xnrm ) /
301 $ max( smlnum, xnorm ) / eps
302 IF( res.GT.rmax ) THEN
303 lmax = knt
304 rmax = res
305 END IF
306 90 CONTINUE
307 100 CONTINUE
308 110 CONTINUE
309 120 CONTINUE
310 130 CONTINUE
311
312 n1 = 2
313 n2 = 2
314 DO 200 itr = 1, 8
315 DO 190 itrscl = 1, 3
316 DO 180 itl = 1, 8
317 DO 170 itlscl = 1, 3
318 DO 160 ib1 = 1, 3
319 DO 150 ib2 = 1, 3
320 DO 140 ib3 = 1, 3
321 b( 1, 1 ) = val( ib1 )
322 b( 2, 1 ) = -four*val( ib2 )
323 b( 1, 2 ) = -two*val( ib3 )
324 b( 2, 2 ) = eight*
325 $ min( val( ib1 ), val
326 $ ( ib2 ), val( ib3 ) )
327 tr( 1, 1 ) = itval( 1, 1, itr )*
328 $ val( itrscl )
329 tr( 2, 1 ) = itval( 2, 1, itr )*
330 $ val( itrscl )
331 tr( 1, 2 ) = itval( 1, 2, itr )*
332 $ val( itrscl )
333 tr( 2, 2 ) = itval( 2, 2, itr )*
334 $ val( itrscl )
335 tl( 1, 1 ) = itval( 1, 1, itl )*
336 $ val( itlscl )
337 tl( 2, 1 ) = itval( 2, 1, itl )*
338 $ val( itlscl )
339 tl( 1, 2 ) = itval( 1, 2, itl )*
340 $ val( itlscl )
341 tl( 2, 2 ) = itval( 2, 2, itl )*
342 $ val( itlscl )
343 knt = knt + 1
344 CALL dlasy2( ltranl, ltranr, isgn,
345 $ n1, n2, tl, 2, tr, 2,
346 $ b, 2, scale, x, 2,
347 $ xnorm, info )
348 IF( info.NE.0 )
349 $ ninfo = ninfo + 1
350 IF( ltranr ) THEN
351 tmp = tr( 1, 2 )
352 tr( 1, 2 ) = tr( 2, 1 )
353 tr( 2, 1 ) = tmp
354 END IF
355 IF( ltranl ) THEN
356 tmp = tl( 1, 2 )
357 tl( 1, 2 ) = tl( 2, 1 )
358 tl( 2, 1 ) = tmp
359 END IF
360 tnrm = abs( tr( 1, 1 ) ) +
361 $ abs( tr( 2, 1 ) ) +
362 $ abs( tr( 1, 2 ) ) +
363 $ abs( tr( 2, 2 ) ) +
364 $ abs( tl( 1, 1 ) ) +
365 $ abs( tl( 2, 1 ) ) +
366 $ abs( tl( 1, 2 ) ) +
367 $ abs( tl( 2, 2 ) )
368 xnrm = max( abs( x( 1, 1 ) )+
369 $ abs( x( 1, 2 ) ),
370 $ abs( x( 2, 1 ) )+
371 $ abs( x( 2, 2 ) ) )
372 res = abs( ( ( tl( 1, 1 )+sgn*tr( 1,
373 $ 1 ) ) )*( x( 1, 1 ) )+
374 $ ( sgn*tr( 2, 1 ) )*
375 $ ( x( 1, 2 ) )+( tl( 1, 2 ) )*
376 $ ( x( 2, 1 ) )-
377 $ ( scale*b( 1, 1 ) ) )
378 res = res + abs( ( tl( 1, 1 ) )*
379 $ ( x( 1, 2 ) )+
380 $ ( sgn*tr( 1, 2 ) )*
381 $ ( x( 1, 1 ) )+
382 $ ( sgn*tr( 2, 2 ) )*
383 $ ( x( 1, 2 ) )+( tl( 1, 2 ) )*
384 $ ( x( 2, 2 ) )-
385 $ ( scale*b( 1, 2 ) ) )
386 res = res + abs( ( tl( 2, 1 ) )*
387 $ ( x( 1, 1 ) )+
388 $ ( sgn*tr( 1, 1 ) )*
389 $ ( x( 2, 1 ) )+
390 $ ( sgn*tr( 2, 1 ) )*
391 $ ( x( 2, 2 ) )+( tl( 2, 2 ) )*
392 $ ( x( 2, 1 ) )-
393 $ ( scale*b( 2, 1 ) ) )
394 res = res + abs( ( ( tl( 2,
395 $ 2 )+sgn*tr( 2, 2 ) ) )*
396 $ ( x( 2, 2 ) )+
397 $ ( sgn*tr( 1, 2 ) )*
398 $ ( x( 2, 1 ) )+( tl( 2, 1 ) )*
399 $ ( x( 1, 2 ) )-
400 $ ( scale*b( 2, 2 ) ) )
401 den = max( smlnum, smlnum*xnrm,
402 $ ( tnrm*eps )*xnrm )
403 res = res / den
404 IF( scale.GT.one )
405 $ res = res + one / eps
406 res = res + abs( xnorm-xnrm ) /
407 $ max( smlnum, xnorm ) / eps
408 IF( res.GT.rmax ) THEN
409 lmax = knt
410 rmax = res
411 END IF
412 140 CONTINUE
413 150 CONTINUE
414 160 CONTINUE
415 170 CONTINUE
416 180 CONTINUE
417 190 CONTINUE
418 200 CONTINUE
419 210 CONTINUE
420 220 CONTINUE
421 230 CONTINUE
422
423 RETURN
424
425
426
double precision function dlamch(cmach)
DLAMCH
subroutine dlasy2(ltranl, ltranr, isgn, n1, n2, tl, ldtl, tr, ldtr, b, ldb, scale, x, ldx, xnorm, info)
DLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.