182
183
184
185
186
187
188 CHARACTER UPLO
189 INTEGER INFO, KB, LDA, LDW, N, NB
190
191
192 INTEGER IPIV( * )
193 REAL A( LDA, * ), W( LDW, * )
194
195
196
197
198
199 REAL ZERO, ONE
200 parameter( zero = 0.0e+0, one = 1.0e+0 )
201 REAL EIGHT, SEVTEN
202 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
203
204
205 LOGICAL DONE
206 INTEGER IMAX, ITEMP, J, JJ, JMAX, JP1, JP2, K, KK,
207 $ KW, KKW, KP, KSTEP, P, II
208
209 REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
210 $ STEMP, R1, ROWMAX, T, SFMIN
211
212
213 LOGICAL LSAME
214 INTEGER ISAMAX
215 REAL SLAMCH
217
218
220
221
222 INTRINSIC abs, max, min, sqrt
223
224
225
226 info = 0
227
228
229
230 alpha = ( one+sqrt( sevten ) ) / eight
231
232
233
235
236 IF(
lsame( uplo,
'U' ) )
THEN
237
238
239
240
241
242
243
244 k = n
245 10 CONTINUE
246
247
248
249 kw = nb + k - n
250
251
252
253 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
254 $ GO TO 30
255
256 kstep = 1
257 p = k
258
259
260
261 CALL scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
262 IF( k.LT.n )
263 $
CALL sgemv(
'No transpose', k, n-k, -one, a( 1, k+1 ),
264 $ lda, w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
265
266
267
268
269 absakk = abs( w( k, kw ) )
270
271
272
273
274
275 IF( k.GT.1 ) THEN
276 imax =
isamax( k-1, w( 1, kw ), 1 )
277 colmax = abs( w( imax, kw ) )
278 ELSE
279 colmax = zero
280 END IF
281
282 IF( max( absakk, colmax ).EQ.zero ) THEN
283
284
285
286 IF( info.EQ.0 )
287 $ info = k
288 kp = k
289 CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
290 ELSE
291
292
293
294
295
296
297
298
299 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
300
301
302
303 kp = k
304
305 ELSE
306
307 done = .false.
308
309
310
311 12 CONTINUE
312
313
314
315
316
317
318 CALL scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ),
319 $ 1 )
320 CALL scopy( k-imax, a( imax, imax+1 ), lda,
321 $ w( imax+1, kw-1 ), 1 )
322
323 IF( k.LT.n )
324 $
CALL sgemv(
'No transpose', k, n-k, -one,
325 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
326 $ one, w( 1, kw-1 ), 1 )
327
328
329
330
331
332 IF( imax.NE.k ) THEN
333 jmax = imax +
isamax( k-imax, w( imax+1, kw-1 ),
334 $ 1 )
335 rowmax = abs( w( jmax, kw-1 ) )
336 ELSE
337 rowmax = zero
338 END IF
339
340 IF( imax.GT.1 ) THEN
341 itemp =
isamax( imax-1, w( 1, kw-1 ), 1 )
342 stemp = abs( w( itemp, kw-1 ) )
343 IF( stemp.GT.rowmax ) THEN
344 rowmax = stemp
345 jmax = itemp
346 END IF
347 END IF
348
349
350
351
352
353 IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
354 $ THEN
355
356
357
358
359 kp = imax
360
361
362
363 CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
364
365 done = .true.
366
367
368
369
370 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
371 $ THEN
372
373
374
375
376 kp = imax
377 kstep = 2
378 done = .true.
379 ELSE
380
381
382
383 p = imax
384 colmax = rowmax
385 imax = jmax
386
387
388
389 CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
390
391 END IF
392
393
394
395 IF( .NOT. done ) GOTO 12
396
397 END IF
398
399
400
401 kk = k - kstep + 1
402
403
404
405 kkw = nb + kk - n
406
407 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
408
409
410
411 CALL scopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
412 CALL scopy( p, a( 1, k ), 1, a( 1, p ), 1 )
413
414
415
416
417 CALL sswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
418 CALL sswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
419 $ ldw )
420 END IF
421
422
423
424 IF( kp.NE.kk ) THEN
425
426
427
428 a( kp, k ) = a( kk, k )
429 CALL scopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
430 $ lda )
431 CALL scopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
432
433
434
435
436 CALL sswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ),
437 $ lda )
438 CALL sswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
439 $ ldw )
440 END IF
441
442 IF( kstep.EQ.1 ) THEN
443
444
445
446
447
448
449
450
451
452 CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
453 IF( k.GT.1 ) THEN
454 IF( abs( a( k, k ) ).GE.sfmin ) THEN
455 r1 = one / a( k, k )
456 CALL sscal( k-1, r1, a( 1, k ), 1 )
457 ELSE IF( a( k, k ).NE.zero ) THEN
458 DO 14 ii = 1, k - 1
459 a( ii, k ) = a( ii, k ) / a( k, k )
460 14 CONTINUE
461 END IF
462 END IF
463
464 ELSE
465
466
467
468
469
470
471
472
473
474 IF( k.GT.2 ) THEN
475
476
477
478 d12 = w( k-1, kw )
479 d11 = w( k, kw ) / d12
480 d22 = w( k-1, kw-1 ) / d12
481 t = one / ( d11*d22-one )
482 DO 20 j = 1, k - 2
483 a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
484 $ d12 )
485 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
486 $ d12 )
487 20 CONTINUE
488 END IF
489
490
491
492 a( k-1, k-1 ) = w( k-1, kw-1 )
493 a( k-1, k ) = w( k-1, kw )
494 a( k, k ) = w( k, kw )
495 END IF
496 END IF
497
498
499
500 IF( kstep.EQ.1 ) THEN
501 ipiv( k ) = kp
502 ELSE
503 ipiv( k ) = -p
504 ipiv( k-1 ) = -kp
505 END IF
506
507
508
509 k = k - kstep
510 GO TO 10
511
512 30 CONTINUE
513
514
515
516
517
518 CALL sgemmtr(
'Upper',
'No transpose',
'Transpose', k, n-k,
519 $ -one, a( 1, k+1 ), lda, w( 1, kw+1 ), ldw,
520 $ one, a( 1, 1 ), lda )
521
522
523
524
525 j = k + 1
526 60 CONTINUE
527
528 kstep = 1
529 jp1 = 1
530 jj = j
531 jp2 = ipiv( j )
532 IF( jp2.LT.0 ) THEN
533 jp2 = -jp2
534 j = j + 1
535 jp1 = -ipiv( j )
536 kstep = 2
537 END IF
538
539 j = j + 1
540 IF( jp2.NE.jj .AND. j.LE.n )
541 $
CALL sswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
542 jj = j - 1
543 IF( jp1.NE.jj .AND. kstep.EQ.2 )
544 $
CALL sswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
545 IF( j.LE.n )
546 $ GO TO 60
547
548
549
550 kb = n - k
551
552 ELSE
553
554
555
556
557
558
559
560 k = 1
561 70 CONTINUE
562
563
564
565 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
566 $ GO TO 90
567
568 kstep = 1
569 p = k
570
571
572
573 CALL scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
574 IF( k.GT.1 )
575 $
CALL sgemv(
'No transpose', n-k+1, k-1, -one, a( k, 1 ),
576 $ lda, w( k, 1 ), ldw, one, w( k, k ), 1 )
577
578
579
580
581 absakk = abs( w( k, k ) )
582
583
584
585
586
587 IF( k.LT.n ) THEN
588 imax = k +
isamax( n-k, w( k+1, k ), 1 )
589 colmax = abs( w( imax, k ) )
590 ELSE
591 colmax = zero
592 END IF
593
594 IF( max( absakk, colmax ).EQ.zero ) THEN
595
596
597
598 IF( info.EQ.0 )
599 $ info = k
600 kp = k
601 CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
602 ELSE
603
604
605
606
607
608
609
610
611 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
612
613
614
615 kp = k
616
617 ELSE
618
619 done = .false.
620
621
622
623 72 CONTINUE
624
625
626
627
628
629
630 CALL scopy( imax-k, a( imax, k ), lda, w( k, k+1 ),
631 $ 1)
632 CALL scopy( n-imax+1, a( imax, imax ), 1,
633 $ w( imax, k+1 ), 1 )
634 IF( k.GT.1 )
635 $
CALL sgemv(
'No transpose', n-k+1, k-1, -one,
636 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
637 $ one, w( k, k+1 ), 1 )
638
639
640
641
642
643 IF( imax.NE.k ) THEN
644 jmax = k - 1 +
isamax( imax-k, w( k, k+1 ), 1 )
645 rowmax = abs( w( jmax, k+1 ) )
646 ELSE
647 rowmax = zero
648 END IF
649
650 IF( imax.LT.n ) THEN
651 itemp = imax +
isamax( n-imax, w( imax+1, k+1 ),
652 $ 1)
653 stemp = abs( w( itemp, k+1 ) )
654 IF( stemp.GT.rowmax ) THEN
655 rowmax = stemp
656 jmax = itemp
657 END IF
658 END IF
659
660
661
662
663
664 IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
665 $ THEN
666
667
668
669
670 kp = imax
671
672
673
674 CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ),
675 $ 1 )
676
677 done = .true.
678
679
680
681
682 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
683 $ THEN
684
685
686
687
688 kp = imax
689 kstep = 2
690 done = .true.
691 ELSE
692
693
694
695 p = imax
696 colmax = rowmax
697 imax = jmax
698
699
700
701 CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ),
702 $ 1 )
703
704 END IF
705
706
707
708 IF( .NOT. done ) GOTO 72
709
710 END IF
711
712
713
714 kk = k + kstep - 1
715
716 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
717
718
719
720 CALL scopy( p-k, a( k, k ), 1, a( p, k ), lda )
721 CALL scopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
722
723
724
725
726 CALL sswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
727 CALL sswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
728 END IF
729
730
731
732 IF( kp.NE.kk ) THEN
733
734
735
736 a( kp, k ) = a( kk, k )
737 CALL scopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ),
738 $ lda )
739 CALL scopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
740
741
742
743 CALL sswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
744 CALL sswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
745 END IF
746
747 IF( kstep.EQ.1 ) THEN
748
749
750
751
752
753
754
755
756
757 CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
758 IF( k.LT.n ) THEN
759 IF( abs( a( k, k ) ).GE.sfmin ) THEN
760 r1 = one / a( k, k )
761 CALL sscal( n-k, r1, a( k+1, k ), 1 )
762 ELSE IF( a( k, k ).NE.zero ) THEN
763 DO 74 ii = k + 1, n
764 a( ii, k ) = a( ii, k ) / a( k, k )
765 74 CONTINUE
766 END IF
767 END IF
768
769 ELSE
770
771
772
773
774
775
776
777
778 IF( k.LT.n-1 ) THEN
779
780
781
782 d21 = w( k+1, k )
783 d11 = w( k+1, k+1 ) / d21
784 d22 = w( k, k ) / d21
785 t = one / ( d11*d22-one )
786 DO 80 j = k + 2, n
787 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
788 $ d21 )
789 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
790 $ d21 )
791 80 CONTINUE
792 END IF
793
794
795
796 a( k, k ) = w( k, k )
797 a( k+1, k ) = w( k+1, k )
798 a( k+1, k+1 ) = w( k+1, k+1 )
799 END IF
800 END IF
801
802
803
804 IF( kstep.EQ.1 ) THEN
805 ipiv( k ) = kp
806 ELSE
807 ipiv( k ) = -p
808 ipiv( k+1 ) = -kp
809 END IF
810
811
812
813 k = k + kstep
814 GO TO 70
815
816 90 CONTINUE
817
818
819
820
821
822 CALL sgemmtr(
'Lower',
'No transpose',
'Transpose', n-k+1,
823 $ k-1, -one, a( k, 1 ), lda, w( k, 1 ), ldw,
824 $ one, a( k, k ), lda )
825
826
827
828
829 j = k - 1
830 120 CONTINUE
831
832 kstep = 1
833 jp1 = 1
834 jj = j
835 jp2 = ipiv( j )
836 IF( jp2.LT.0 ) THEN
837 jp2 = -jp2
838 j = j - 1
839 jp1 = -ipiv( j )
840 kstep = 2
841 END IF
842
843 j = j - 1
844 IF( jp2.NE.jj .AND. j.GE.1 )
845 $
CALL sswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
846 jj = j + 1
847 IF( jp1.NE.jj .AND. kstep.EQ.2 )
848 $
CALL sswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
849 IF( j.GE.1 )
850 $ GO TO 120
851
852
853
854 kb = k - 1
855
856 END IF
857 RETURN
858
859
860
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
subroutine sgemmtr(uplo, transa, transb, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMMTR
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
integer function isamax(n, sx, incx)
ISAMAX
real function slamch(cmach)
SLAMCH
logical function lsame(ca, cb)
LSAME
subroutine sscal(n, sa, sx, incx)
SSCAL
subroutine sswap(n, sx, incx, sy, incy)
SSWAP