91
92
93
94
95
96
97 INTEGER KNT, LMAX
98 REAL RMAX
99
100
101 INTEGER NINFO( 2 )
102
103
104
105
106
107 REAL ZERO, HALF, ONE
108 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0 )
109 REAL TWO, THREE, FOUR
110 parameter( two = 2.0e0, three = 3.0e0, four = 4.0e0 )
111 REAL SEVEN, TEN
112 parameter( seven = 7.0e0, ten = 10.0e0 )
113 REAL TWNONE
114 parameter( twnone = 21.0e0 )
115
116
117 INTEGER IA, IB, ICA, ID1, ID2, INFO, ISMIN, ITRANS,
118 $ IWI, IWR, NA, NW
119 REAL BIGNUM, CA, D1, D2, DEN, EPS, RES, SCALE, SMIN,
120 $ SMLNUM, TMP, UNFL, WI, WR, XNORM
121
122
123 LOGICAL LTRANS( 0: 1 )
124 REAL A( 2, 2 ), B( 2, 2 ), VAB( 3 ), VCA( 5 ),
125 $ VDD( 4 ), VSMIN( 4 ), VWI( 4 ), VWR( 4 ),
126 $ X( 2, 2 )
127
128
129 REAL SLAMCH
131
132
134
135
136 INTRINSIC abs, max, sqrt
137
138
139 DATA ltrans / .false., .true. /
140
141
142
143
144
147 smlnum =
slamch(
'S' ) / eps
148 bignum = one / smlnum
149
150
151
152 vsmin( 1 ) = smlnum
153 vsmin( 2 ) = eps
154 vsmin( 3 ) = one / ( ten*ten )
155 vsmin( 4 ) = one / eps
156 vab( 1 ) = sqrt( smlnum )
157 vab( 2 ) = one
158 vab( 3 ) = sqrt( bignum )
159 vwr( 1 ) = zero
160 vwr( 2 ) = half
161 vwr( 3 ) = two
162 vwr( 4 ) = one
163 vwi( 1 ) = smlnum
164 vwi( 2 ) = eps
165 vwi( 3 ) = one
166 vwi( 4 ) = two
167 vdd( 1 ) = sqrt( smlnum )
168 vdd( 2 ) = one
169 vdd( 3 ) = two
170 vdd( 4 ) = sqrt( bignum )
171 vca( 1 ) = zero
172 vca( 2 ) = sqrt( smlnum )
173 vca( 3 ) = eps
174 vca( 4 ) = half
175 vca( 5 ) = one
176
177 knt = 0
178 ninfo( 1 ) = 0
179 ninfo( 2 ) = 0
180 lmax = 0
181 rmax = zero
182
183
184
185 DO 190 id1 = 1, 4
186 d1 = vdd( id1 )
187 DO 180 id2 = 1, 4
188 d2 = vdd( id2 )
189 DO 170 ica = 1, 5
190 ca = vca( ica )
191 DO 160 itrans = 0, 1
192 DO 150 ismin = 1, 4
193 smin = vsmin( ismin )
194
195 na = 1
196 nw = 1
197 DO 30 ia = 1, 3
198 a( 1, 1 ) = vab( ia )
199 DO 20 ib = 1, 3
200 b( 1, 1 ) = vab( ib )
201 DO 10 iwr = 1, 4
202 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
203 $ one ) THEN
204 wr = vwr( iwr )*a( 1, 1 )
205 ELSE
206 wr = vwr( iwr )
207 END IF
208 wi = zero
209 CALL slaln2( ltrans( itrans ), na, nw,
210 $ smin, ca, a, 2, d1, d2, b, 2,
211 $ wr, wi, x, 2, scale, xnorm,
212 $ info )
213 IF( info.LT.0 )
214 $ ninfo( 1 ) = ninfo( 1 ) + 1
215 IF( info.GT.0 )
216 $ ninfo( 2 ) = ninfo( 2 ) + 1
217 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
218 $ x( 1, 1 )-scale*b( 1, 1 ) )
219 IF( info.EQ.0 ) THEN
220 den = max( eps*( abs( ( ca*a( 1,
221 $ 1 )-wr*d1 )*x( 1, 1 ) ) ),
222 $ smlnum )
223 ELSE
224 den = max( smin*abs( x( 1, 1 ) ),
225 $ smlnum )
226 END IF
227 res = res / den
228 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
229 $ abs( b( 1, 1 ) ).LE.smlnum*
230 $ abs( ca*a( 1, 1 )-wr*d1 ) )res = zero
231 IF( scale.GT.one )
232 $ res = res + one / eps
233 res = res + abs( xnorm-abs( x( 1, 1 ) ) )
234 $ / max( smlnum, xnorm ) / eps
235 IF( info.NE.0 .AND. info.NE.1 )
236 $ res = res + one / eps
237 knt = knt + 1
238 IF( res.GT.rmax ) THEN
239 lmax = knt
240 rmax = res
241 END IF
242 10 CONTINUE
243 20 CONTINUE
244 30 CONTINUE
245
246 na = 1
247 nw = 2
248 DO 70 ia = 1, 3
249 a( 1, 1 ) = vab( ia )
250 DO 60 ib = 1, 3
251 b( 1, 1 ) = vab( ib )
252 b( 1, 2 ) = -half*vab( ib )
253 DO 50 iwr = 1, 4
254 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
255 $ one ) THEN
256 wr = vwr( iwr )*a( 1, 1 )
257 ELSE
258 wr = vwr( iwr )
259 END IF
260 DO 40 iwi = 1, 4
261 IF( d1.EQ.one .AND. d2.EQ.one .AND.
262 $ ca.EQ.one ) THEN
263 wi = vwi( iwi )*a( 1, 1 )
264 ELSE
265 wi = vwi( iwi )
266 END IF
267 CALL slaln2( ltrans( itrans ), na, nw,
268 $ smin, ca, a, 2, d1, d2, b,
269 $ 2, wr, wi, x, 2, scale,
270 $ xnorm, info )
271 IF( info.LT.0 )
272 $ ninfo( 1 ) = ninfo( 1 ) + 1
273 IF( info.GT.0 )
274 $ ninfo( 2 ) = ninfo( 2 ) + 1
275 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
276 $ x( 1, 1 )+( wi*d1 )*x( 1, 2 )-
277 $ scale*b( 1, 1 ) )
278 res = res + abs( ( -wi*d1 )*x( 1, 1 )+
279 $ ( ca*a( 1, 1 )-wr*d1 )*x( 1, 2 )-
280 $ scale*b( 1, 2 ) )
281 IF( info.EQ.0 ) THEN
282 den = max( eps*( max( abs( ca*a( 1,
283 $ 1 )-wr*d1 ), abs( d1*wi ) )*
284 $ ( abs( x( 1, 1 ) )+abs( x( 1,
285 $ 2 ) ) ) ), smlnum )
286 ELSE
287 den = max( smin*( abs( x( 1,
288 $ 1 ) )+abs( x( 1, 2 ) ) ),
289 $ smlnum )
290 END IF
291 res = res / den
292 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
293 $ abs( x( 1, 2 ) ).LT.unfl .AND.
294 $ abs( b( 1, 1 ) ).LE.smlnum*
295 $ abs( ca*a( 1, 1 )-wr*d1 ) )
296 $ res = zero
297 IF( scale.GT.one )
298 $ res = res + one / eps
299 res = res + abs( xnorm-
300 $ abs( x( 1, 1 ) )-
301 $ abs( x( 1, 2 ) ) ) /
302 $ max( smlnum, xnorm ) / eps
303 IF( info.NE.0 .AND. info.NE.1 )
304 $ res = res + one / eps
305 knt = knt + 1
306 IF( res.GT.rmax ) THEN
307 lmax = knt
308 rmax = res
309 END IF
310 40 CONTINUE
311 50 CONTINUE
312 60 CONTINUE
313 70 CONTINUE
314
315 na = 2
316 nw = 1
317 DO 100 ia = 1, 3
318 a( 1, 1 ) = vab( ia )
319 a( 1, 2 ) = -three*vab( ia )
320 a( 2, 1 ) = -seven*vab( ia )
321 a( 2, 2 ) = twnone*vab( ia )
322 DO 90 ib = 1, 3
323 b( 1, 1 ) = vab( ib )
324 b( 2, 1 ) = -two*vab( ib )
325 DO 80 iwr = 1, 4
326 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
327 $ one ) THEN
328 wr = vwr( iwr )*a( 1, 1 )
329 ELSE
330 wr = vwr( iwr )
331 END IF
332 wi = zero
333 CALL slaln2( ltrans( itrans ), na, nw,
334 $ smin, ca, a, 2, d1, d2, b, 2,
335 $ wr, wi, x, 2, scale, xnorm,
336 $ info )
337 IF( info.LT.0 )
338 $ ninfo( 1 ) = ninfo( 1 ) + 1
339 IF( info.GT.0 )
340 $ ninfo( 2 ) = ninfo( 2 ) + 1
341 IF( itrans.EQ.1 ) THEN
342 tmp = a( 1, 2 )
343 a( 1, 2 ) = a( 2, 1 )
344 a( 2, 1 ) = tmp
345 END IF
346 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
347 $ x( 1, 1 )+( ca*a( 1, 2 ) )*
348 $ x( 2, 1 )-scale*b( 1, 1 ) )
349 res = res + abs( ( ca*a( 2, 1 ) )*
350 $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
351 $ x( 2, 1 )-scale*b( 2, 1 ) )
352 IF( info.EQ.0 ) THEN
353 den = max( eps*( max( abs( ca*a( 1,
354 $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
355 $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
356 $ 2 )-wr*d2 ) )*max( abs( x( 1,
357 $ 1 ) ), abs( x( 2, 1 ) ) ) ),
358 $ smlnum )
359 ELSE
360 den = max( eps*( max( smin / eps,
361 $ max( abs( ca*a( 1,
362 $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
363 $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
364 $ 2 )-wr*d2 ) ) )*max( abs( x( 1,
365 $ 1 ) ), abs( x( 2, 1 ) ) ) ),
366 $ smlnum )
367 END IF
368 res = res / den
369 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
370 $ abs( x( 2, 1 ) ).LT.unfl .AND.
371 $ abs( b( 1, 1 ) )+abs( b( 2, 1 ) ).LE.
372 $ smlnum*( abs( ca*a( 1,
373 $ 1 )-wr*d1 )+abs( ca*a( 1,
374 $ 2 ) )+abs( ca*a( 2,
375 $ 1 ) )+abs( ca*a( 2, 2 )-wr*d2 ) ) )
376 $ res = zero
377 IF( scale.GT.one )
378 $ res = res + one / eps
379 res = res + abs( xnorm-
380 $ max( abs( x( 1, 1 ) ), abs( x( 2,
381 $ 1 ) ) ) ) / max( smlnum, xnorm ) /
382 $ eps
383 IF( info.NE.0 .AND. info.NE.1 )
384 $ res = res + one / eps
385 knt = knt + 1
386 IF( res.GT.rmax ) THEN
387 lmax = knt
388 rmax = res
389 END IF
390 80 CONTINUE
391 90 CONTINUE
392 100 CONTINUE
393
394 na = 2
395 nw = 2
396 DO 140 ia = 1, 3
397 a( 1, 1 ) = vab( ia )*two
398 a( 1, 2 ) = -three*vab( ia )
399 a( 2, 1 ) = -seven*vab( ia )
400 a( 2, 2 ) = twnone*vab( ia )
401 DO 130 ib = 1, 3
402 b( 1, 1 ) = vab( ib )
403 b( 2, 1 ) = -two*vab( ib )
404 b( 1, 2 ) = four*vab( ib )
405 b( 2, 2 ) = -seven*vab( ib )
406 DO 120 iwr = 1, 4
407 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
408 $ one ) THEN
409 wr = vwr( iwr )*a( 1, 1 )
410 ELSE
411 wr = vwr( iwr )
412 END IF
413 DO 110 iwi = 1, 4
414 IF( d1.EQ.one .AND. d2.EQ.one .AND.
415 $ ca.EQ.one ) THEN
416 wi = vwi( iwi )*a( 1, 1 )
417 ELSE
418 wi = vwi( iwi )
419 END IF
420 CALL slaln2( ltrans( itrans ), na, nw,
421 $ smin, ca, a, 2, d1, d2, b,
422 $ 2, wr, wi, x, 2, scale,
423 $ xnorm, info )
424 IF( info.LT.0 )
425 $ ninfo( 1 ) = ninfo( 1 ) + 1
426 IF( info.GT.0 )
427 $ ninfo( 2 ) = ninfo( 2 ) + 1
428 IF( itrans.EQ.1 ) THEN
429 tmp = a( 1, 2 )
430 a( 1, 2 ) = a( 2, 1 )
431 a( 2, 1 ) = tmp
432 END IF
433 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
434 $ x( 1, 1 )+( ca*a( 1, 2 ) )*
435 $ x( 2, 1 )+( wi*d1 )*x( 1, 2 )-
436 $ scale*b( 1, 1 ) )
437 res = res + abs( ( ca*a( 1,
438 $ 1 )-wr*d1 )*x( 1, 2 )+
439 $ ( ca*a( 1, 2 ) )*x( 2, 2 )-
440 $ ( wi*d1 )*x( 1, 1 )-scale*
441 $ b( 1, 2 ) )
442 res = res + abs( ( ca*a( 2, 1 ) )*
443 $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
444 $ x( 2, 1 )+( wi*d2 )*x( 2, 2 )-
445 $ scale*b( 2, 1 ) )
446 res = res + abs( ( ca*a( 2, 1 ) )*
447 $ x( 1, 2 )+( ca*a( 2, 2 )-wr*d2 )*
448 $ x( 2, 2 )-( wi*d2 )*x( 2, 1 )-
449 $ scale*b( 2, 2 ) )
450 IF( info.EQ.0 ) THEN
451 den = max( eps*( max( abs( ca*a( 1,
452 $ 1 )-wr*d1 )+abs( ca*a( 1,
453 $ 2 ) )+abs( wi*d1 ),
454 $ abs( ca*a( 2,
455 $ 1 ) )+abs( ca*a( 2,
456 $ 2 )-wr*d2 )+abs( wi*d2 ) )*
457 $ max( abs( x( 1,
458 $ 1 ) )+abs( x( 2, 1 ) ),
459 $ abs( x( 1, 2 ) )+abs( x( 2,
460 $ 2 ) ) ) ), smlnum )
461 ELSE
462 den = max( eps*( max( smin / eps,
463 $ max( abs( ca*a( 1,
464 $ 1 )-wr*d1 )+abs( ca*a( 1,
465 $ 2 ) )+abs( wi*d1 ),
466 $ abs( ca*a( 2,
467 $ 1 ) )+abs( ca*a( 2,
468 $ 2 )-wr*d2 )+abs( wi*d2 ) ) )*
469 $ max( abs( x( 1,
470 $ 1 ) )+abs( x( 2, 1 ) ),
471 $ abs( x( 1, 2 ) )+abs( x( 2,
472 $ 2 ) ) ) ), smlnum )
473 END IF
474 res = res / den
475 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
476 $ abs( x( 2, 1 ) ).LT.unfl .AND.
477 $ abs( x( 1, 2 ) ).LT.unfl .AND.
478 $ abs( x( 2, 2 ) ).LT.unfl .AND.
479 $ abs( b( 1, 1 ) )+
480 $ abs( b( 2, 1 ) ).LE.smlnum*
481 $ ( abs( ca*a( 1, 1 )-wr*d1 )+
482 $ abs( ca*a( 1, 2 ) )+abs( ca*a( 2,
483 $ 1 ) )+abs( ca*a( 2,
484 $ 2 )-wr*d2 )+abs( wi*d2 )+abs( wi*
485 $ d1 ) ) )res = zero
486 IF( scale.GT.one )
487 $ res = res + one / eps
488 res = res + abs( xnorm-
489 $ max( abs( x( 1, 1 ) )+abs( x( 1,
490 $ 2 ) ), abs( x( 2,
491 $ 1 ) )+abs( x( 2, 2 ) ) ) ) /
492 $ max( smlnum, xnorm ) / eps
493 IF( info.NE.0 .AND. info.NE.1 )
494 $ res = res + one / eps
495 knt = knt + 1
496 IF( res.GT.rmax ) THEN
497 lmax = knt
498 rmax = res
499 END IF
500 110 CONTINUE
501 120 CONTINUE
502 130 CONTINUE
503 140 CONTINUE
504 150 CONTINUE
505 160 CONTINUE
506 170 CONTINUE
507 180 CONTINUE
508 190 CONTINUE
509
510 RETURN
511
512
513
subroutine slaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
real function slamch(cmach)
SLAMCH