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