165
166
167
168
169
170
171 LOGICAL LREAL, LTRAN
172 INTEGER INFO, LDT, N
173 REAL SCALE, W
174
175
176 REAL B( * ), T( LDT, * ), WORK( * ), X( * )
177
178
179
180
181
182 REAL ZERO, ONE
183 parameter( zero = 0.0e+0, one = 1.0e+0 )
184
185
186 LOGICAL NOTRAN
187 INTEGER I, IERR, J, J1, J2, JNEXT, K, N1, N2
188 REAL BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
189 $ SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
190
191
192 REAL D( 2, 2 ), V( 2, 2 )
193
194
195 INTEGER ISAMAX
196 REAL SASUM, SDOT, SLAMCH, SLANGE
198
199
201
202
203 INTRINSIC abs, max
204
205
206
207
208
209 notran = .NOT.ltran
210 info = 0
211
212
213
214 IF( n.EQ.0 )
215 $ RETURN
216
217
218
220 smlnum =
slamch(
'S' ) / eps
221 bignum = one / smlnum
222
223 xnorm =
slange(
'M', n, n, t, ldt, d )
224 IF( .NOT.lreal )
225 $ xnorm = max( xnorm, abs( w ),
slange(
'M', n, 1, b, n, d ) )
226 smin = max( smlnum, eps*xnorm )
227
228
229
230
231 work( 1 ) = zero
232 DO 10 j = 2, n
233 work( j ) =
sasum( j-1, t( 1, j ), 1 )
234 10 CONTINUE
235
236 IF( .NOT.lreal ) THEN
237 DO 20 i = 2, n
238 work( i ) = work( i ) + abs( b( i ) )
239 20 CONTINUE
240 END IF
241
242 n2 = 2*n
243 n1 = n
244 IF( .NOT.lreal )
245 $ n1 = n2
247 xmax = abs( x( k ) )
248 scale = one
249
250 IF( xmax.GT.bignum ) THEN
251 scale = bignum / xmax
252 CALL sscal( n1, scale, x, 1 )
253 xmax = bignum
254 END IF
255
256 IF( lreal ) THEN
257
258 IF( notran ) THEN
259
260
261
262 jnext = n
263 DO 30 j = n, 1, -1
264 IF( j.GT.jnext )
265 $ GO TO 30
266 j1 = j
267 j2 = j
268 jnext = j - 1
269 IF( j.GT.1 ) THEN
270 IF( t( j, j-1 ).NE.zero ) THEN
271 j1 = j - 1
272 jnext = j - 2
273 END IF
274 END IF
275
276 IF( j1.EQ.j2 ) THEN
277
278
279
280
281
282
283 xj = abs( x( j1 ) )
284 tjj = abs( t( j1, j1 ) )
285 tmp = t( j1, j1 )
286 IF( tjj.LT.smin ) THEN
287 tmp = smin
288 tjj = smin
289 info = 1
290 END IF
291
292 IF( xj.EQ.zero )
293 $ GO TO 30
294
295 IF( tjj.LT.one ) THEN
296 IF( xj.GT.bignum*tjj ) THEN
297 rec = one / xj
298 CALL sscal( n, rec, x, 1 )
299 scale = scale*rec
300 xmax = xmax*rec
301 END IF
302 END IF
303 x( j1 ) = x( j1 ) / tmp
304 xj = abs( x( j1 ) )
305
306
307
308
309 IF( xj.GT.one ) THEN
310 rec = one / xj
311 IF( work( j1 ).GT.( bignum-xmax )*rec ) THEN
312 CALL sscal( n, rec, x, 1 )
313 scale = scale*rec
314 END IF
315 END IF
316 IF( j1.GT.1 ) THEN
317 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
319 xmax = abs( x( k ) )
320 END IF
321
322 ELSE
323
324
325
326
327
328
329 d( 1, 1 ) = x( j1 )
330 d( 2, 1 ) = x( j2 )
331 CALL slaln2( .false., 2, 1, smin, one, t( j1, j1 ),
332 $ ldt, one, one, d, 2, zero, zero, v, 2,
333 $ scaloc, xnorm, ierr )
334 IF( ierr.NE.0 )
335 $ info = 2
336
337 IF( scaloc.NE.one ) THEN
338 CALL sscal( n, scaloc, x, 1 )
339 scale = scale*scaloc
340 END IF
341 x( j1 ) = v( 1, 1 )
342 x( j2 ) = v( 2, 1 )
343
344
345
346
347 xj = max( abs( v( 1, 1 ) ), abs( v( 2, 1 ) ) )
348 IF( xj.GT.one ) THEN
349 rec = one / xj
350 IF( max( work( j1 ), work( j2 ) ).GT.
351 $ ( bignum-xmax )*rec ) THEN
352 CALL sscal( n, rec, x, 1 )
353 scale = scale*rec
354 END IF
355 END IF
356
357
358
359 IF( j1.GT.1 ) THEN
360 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
361 CALL saxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x, 1 )
363 xmax = abs( x( k ) )
364 END IF
365
366 END IF
367
368 30 CONTINUE
369
370 ELSE
371
372
373
374 jnext = 1
375 DO 40 j = 1, n
376 IF( j.LT.jnext )
377 $ GO TO 40
378 j1 = j
379 j2 = j
380 jnext = j + 1
381 IF( j.LT.n ) THEN
382 IF( t( j+1, j ).NE.zero ) THEN
383 j2 = j + 1
384 jnext = j + 2
385 END IF
386 END IF
387
388 IF( j1.EQ.j2 ) THEN
389
390
391
392
393
394
395 xj = abs( x( j1 ) )
396 IF( xmax.GT.one ) THEN
397 rec = one / xmax
398 IF( work( j1 ).GT.( bignum-xj )*rec ) THEN
399 CALL sscal( n, rec, x, 1 )
400 scale = scale*rec
401 xmax = xmax*rec
402 END IF
403 END IF
404
405 x( j1 ) = x( j1 ) -
sdot( j1-1, t( 1, j1 ), 1, x, 1 )
406
407 xj = abs( x( j1 ) )
408 tjj = abs( t( j1, j1 ) )
409 tmp = t( j1, j1 )
410 IF( tjj.LT.smin ) THEN
411 tmp = smin
412 tjj = smin
413 info = 1
414 END IF
415
416 IF( tjj.LT.one ) THEN
417 IF( xj.GT.bignum*tjj ) THEN
418 rec = one / xj
419 CALL sscal( n, rec, x, 1 )
420 scale = scale*rec
421 xmax = xmax*rec
422 END IF
423 END IF
424 x( j1 ) = x( j1 ) / tmp
425 xmax = max( xmax, abs( x( j1 ) ) )
426
427 ELSE
428
429
430
431
432
433
434 xj = max( abs( x( j1 ) ), abs( x( j2 ) ) )
435 IF( xmax.GT.one ) THEN
436 rec = one / xmax
437 IF( max( work( j2 ), work( j1 ) ).GT.( bignum-xj )*
438 $ rec ) THEN
439 CALL sscal( n, rec, x, 1 )
440 scale = scale*rec
441 xmax = xmax*rec
442 END IF
443 END IF
444
445 d( 1, 1 ) = x( j1 ) -
sdot( j1-1, t( 1, j1 ), 1, x,
446 $ 1 )
447 d( 2, 1 ) = x( j2 ) -
sdot( j1-1, t( 1, j2 ), 1, x,
448 $ 1 )
449
450 CALL slaln2( .true., 2, 1, smin, one, t( j1, j1 ),
451 $ ldt, one, one, d, 2, zero, zero, v, 2,
452 $ scaloc, xnorm, ierr )
453 IF( ierr.NE.0 )
454 $ info = 2
455
456 IF( scaloc.NE.one ) THEN
457 CALL sscal( n, scaloc, x, 1 )
458 scale = scale*scaloc
459 END IF
460 x( j1 ) = v( 1, 1 )
461 x( j2 ) = v( 2, 1 )
462 xmax = max( abs( x( j1 ) ), abs( x( j2 ) ), xmax )
463
464 END IF
465 40 CONTINUE
466 END IF
467
468 ELSE
469
470 sminw = max( eps*abs( w ), smin )
471 IF( notran ) THEN
472
473
474
475 jnext = n
476 DO 70 j = n, 1, -1
477 IF( j.GT.jnext )
478 $ GO TO 70
479 j1 = j
480 j2 = j
481 jnext = j - 1
482 IF( j.GT.1 ) THEN
483 IF( t( j, j-1 ).NE.zero ) THEN
484 j1 = j - 1
485 jnext = j - 2
486 END IF
487 END IF
488
489 IF( j1.EQ.j2 ) THEN
490
491
492
493
494
495 z = w
496 IF( j1.EQ.1 )
497 $ z = b( 1 )
498 xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
499 tjj = abs( t( j1, j1 ) ) + abs( z )
500 tmp = t( j1, j1 )
501 IF( tjj.LT.sminw ) THEN
502 tmp = sminw
503 tjj = sminw
504 info = 1
505 END IF
506
507 IF( xj.EQ.zero )
508 $ GO TO 70
509
510 IF( tjj.LT.one ) THEN
511 IF( xj.GT.bignum*tjj ) THEN
512 rec = one / xj
513 CALL sscal( n2, rec, x, 1 )
514 scale = scale*rec
515 xmax = xmax*rec
516 END IF
517 END IF
518 CALL sladiv( x( j1 ), x( n+j1 ), tmp, z, sr, si )
519 x( j1 ) = sr
520 x( n+j1 ) = si
521 xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
522
523
524
525
526 IF( xj.GT.one ) THEN
527 rec = one / xj
528 IF( work( j1 ).GT.( bignum-xmax )*rec ) THEN
529 CALL sscal( n2, rec, x, 1 )
530 scale = scale*rec
531 END IF
532 END IF
533
534 IF( j1.GT.1 ) THEN
535 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
536 CALL saxpy( j1-1, -x( n+j1 ), t( 1, j1 ), 1,
537 $ x( n+1 ), 1 )
538
539 x( 1 ) = x( 1 ) + b( j1 )*x( n+j1 )
540 x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 )
541
542 xmax = zero
543 DO 50 k = 1, j1 - 1
544 xmax = max( xmax, abs( x( k ) )+
545 $ abs( x( k+n ) ) )
546 50 CONTINUE
547 END IF
548
549 ELSE
550
551
552
553 d( 1, 1 ) = x( j1 )
554 d( 2, 1 ) = x( j2 )
555 d( 1, 2 ) = x( n+j1 )
556 d( 2, 2 ) = x( n+j2 )
557 CALL slaln2( .false., 2, 2, sminw, one, t( j1, j1 ),
558 $ ldt, one, one, d, 2, zero, -w, v, 2,
559 $ scaloc, xnorm, ierr )
560 IF( ierr.NE.0 )
561 $ info = 2
562
563 IF( scaloc.NE.one ) THEN
564 CALL sscal( 2*n, scaloc, x, 1 )
565 scale = scaloc*scale
566 END IF
567 x( j1 ) = v( 1, 1 )
568 x( j2 ) = v( 2, 1 )
569 x( n+j1 ) = v( 1, 2 )
570 x( n+j2 ) = v( 2, 2 )
571
572
573
574
575 xj = max( abs( v( 1, 1 ) )+abs( v( 1, 2 ) ),
576 $ abs( v( 2, 1 ) )+abs( v( 2, 2 ) ) )
577 IF( xj.GT.one ) THEN
578 rec = one / xj
579 IF( max( work( j1 ), work( j2 ) ).GT.
580 $ ( bignum-xmax )*rec ) THEN
581 CALL sscal( n2, rec, x, 1 )
582 scale = scale*rec
583 END IF
584 END IF
585
586
587
588 IF( j1.GT.1 ) THEN
589 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
590 CALL saxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x, 1 )
591
592 CALL saxpy( j1-1, -x( n+j1 ), t( 1, j1 ), 1,
593 $ x( n+1 ), 1 )
594 CALL saxpy( j1-1, -x( n+j2 ), t( 1, j2 ), 1,
595 $ x( n+1 ), 1 )
596
597 x( 1 ) = x( 1 ) + b( j1 )*x( n+j1 ) +
598 $ b( j2 )*x( n+j2 )
599 x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -
600 $ b( j2 )*x( j2 )
601
602 xmax = zero
603 DO 60 k = 1, j1 - 1
604 xmax = max( abs( x( k ) )+abs( x( k+n ) ),
605 $ xmax )
606 60 CONTINUE
607 END IF
608
609 END IF
610 70 CONTINUE
611
612 ELSE
613
614
615
616 jnext = 1
617 DO 80 j = 1, n
618 IF( j.LT.jnext )
619 $ GO TO 80
620 j1 = j
621 j2 = j
622 jnext = j + 1
623 IF( j.LT.n ) THEN
624 IF( t( j+1, j ).NE.zero ) THEN
625 j2 = j + 1
626 jnext = j + 2
627 END IF
628 END IF
629
630 IF( j1.EQ.j2 ) THEN
631
632
633
634
635
636
637 xj = abs( x( j1 ) ) + abs( x( j1+n ) )
638 IF( xmax.GT.one ) THEN
639 rec = one / xmax
640 IF( work( j1 ).GT.( bignum-xj )*rec ) THEN
641 CALL sscal( n2, rec, x, 1 )
642 scale = scale*rec
643 xmax = xmax*rec
644 END IF
645 END IF
646
647 x( j1 ) = x( j1 ) -
sdot( j1-1, t( 1, j1 ), 1, x, 1 )
648 x( n+j1 ) = x( n+j1 ) -
sdot( j1-1, t( 1, j1 ), 1,
649 $ x( n+1 ), 1 )
650 IF( j1.GT.1 ) THEN
651 x( j1 ) = x( j1 ) - b( j1 )*x( n+1 )
652 x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1 )
653 END IF
654 xj = abs( x( j1 ) ) + abs( x( j1+n ) )
655
656 z = w
657 IF( j1.EQ.1 )
658 $ z = b( 1 )
659
660
661
662
663 tjj = abs( t( j1, j1 ) ) + abs( z )
664 tmp = t( j1, j1 )
665 IF( tjj.LT.sminw ) THEN
666 tmp = sminw
667 tjj = sminw
668 info = 1
669 END IF
670
671 IF( tjj.LT.one ) THEN
672 IF( xj.GT.bignum*tjj ) THEN
673 rec = one / xj
674 CALL sscal( n2, rec, x, 1 )
675 scale = scale*rec
676 xmax = xmax*rec
677 END IF
678 END IF
679 CALL sladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si )
680 x( j1 ) = sr
681 x( j1+n ) = si
682 xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax )
683
684 ELSE
685
686
687
688
689
690
691 xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),
692 $ abs( x( j2 ) )+abs( x( n+j2 ) ) )
693 IF( xmax.GT.one ) THEN
694 rec = one / xmax
695 IF( max( work( j1 ), work( j2 ) ).GT.
696 $ ( bignum-xj ) / xmax ) THEN
697 CALL sscal( n2, rec, x, 1 )
698 scale = scale*rec
699 xmax = xmax*rec
700 END IF
701 END IF
702
703 d( 1, 1 ) = x( j1 ) -
sdot( j1-1, t( 1, j1 ), 1, x,
704 $ 1 )
705 d( 2, 1 ) = x( j2 ) -
sdot( j1-1, t( 1, j2 ), 1, x,
706 $ 1 )
707 d( 1, 2 ) = x( n+j1 ) -
sdot( j1-1, t( 1, j1 ), 1,
708 $ x( n+1 ), 1 )
709 d( 2, 2 ) = x( n+j2 ) -
sdot( j1-1, t( 1, j2 ), 1,
710 $ x( n+1 ), 1 )
711 d( 1, 1 ) = d( 1, 1 ) - b( j1 )*x( n+1 )
712 d( 2, 1 ) = d( 2, 1 ) - b( j2 )*x( n+1 )
713 d( 1, 2 ) = d( 1, 2 ) + b( j1 )*x( 1 )
714 d( 2, 2 ) = d( 2, 2 ) + b( j2 )*x( 1 )
715
716 CALL slaln2( .true., 2, 2, sminw, one, t( j1, j1 ),
717 $ ldt, one, one, d, 2, zero, w, v, 2,
718 $ scaloc, xnorm, ierr )
719 IF( ierr.NE.0 )
720 $ info = 2
721
722 IF( scaloc.NE.one ) THEN
723 CALL sscal( n2, scaloc, x, 1 )
724 scale = scaloc*scale
725 END IF
726 x( j1 ) = v( 1, 1 )
727 x( j2 ) = v( 2, 1 )
728 x( n+j1 ) = v( 1, 2 )
729 x( n+j2 ) = v( 2, 2 )
730 xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),
731 $ abs( x( j2 ) )+abs( x( n+j2 ) ), xmax )
732
733 END IF
734
735 80 CONTINUE
736
737 END IF
738
739 END IF
740
741 RETURN
742
743
744
real function sasum(n, sx, incx)
SASUM
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
real function sdot(n, sx, incx, sy, incy)
SDOT
integer function isamax(n, sx, incx)
ISAMAX
subroutine sladiv(a, b, c, d, p, q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
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
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine sscal(n, sa, sx, incx)
SSCAL