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