251
252
253
254
255
256
257 CHARACTER DIRECT, SIDE, STOREV, TRANS
258 INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
259
260
261 COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
262 $ V( LDV, * ), WORK( LDWORK, * )
263
264
265
266
267
268 COMPLEX*16 ONE, ZERO
269 parameter( one = (1.0,0.0), zero = (0.0,0.0) )
270
271
272 INTEGER I, J, MP, NP, KP
273 LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
274
275
276 LOGICAL LSAME
278
279
281
282
283 INTRINSIC conjg
284
285
286
287
288
289 IF( m.LE.0 .OR. n.LE.0 .OR. k.LE.0 .OR. l.LT.0 ) RETURN
290
291 IF(
lsame( storev,
'C' ) )
THEN
292 column = .true.
293 row = .false.
294 ELSE IF (
lsame( storev,
'R' ) )
THEN
295 column = .false.
296 row = .true.
297 ELSE
298 column = .false.
299 row = .false.
300 END IF
301
302 IF(
lsame( side,
'L' ) )
THEN
303 left = .true.
304 right = .false.
305 ELSE IF(
lsame( side,
'R' ) )
THEN
306 left = .false.
307 right = .true.
308 ELSE
309 left = .false.
310 right = .false.
311 END IF
312
313 IF(
lsame( direct,
'F' ) )
THEN
314 forward = .true.
315 backward = .false.
316 ELSE IF(
lsame( direct,
'B' ) )
THEN
317 forward = .false.
318 backward = .true.
319 ELSE
320 forward = .false.
321 backward = .false.
322 END IF
323
324
325
326 IF( column .AND. forward .AND. left ) THEN
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343 mp = min( m-l+1, m )
344 kp = min( l+1, k )
345
346 DO j = 1, n
347 DO i = 1, l
348 work( i, j ) = b( m-l+i, j )
349 END DO
350 END DO
351 CALL ztrmm(
'L',
'U',
'C',
'N', l, n, one, v( mp, 1 ), ldv,
352 $ work, ldwork )
353 CALL zgemm(
'C',
'N', l, n, m-l, one, v, ldv, b, ldb,
354 $ one, work, ldwork )
355 CALL zgemm(
'C',
'N', k-l, n, m, one, v( 1, kp ), ldv,
356 $ b, ldb, zero, work( kp, 1 ), ldwork )
357
358 DO j = 1, n
359 DO i = 1, k
360 work( i, j ) = work( i, j ) + a( i, j )
361 END DO
362 END DO
363
364 CALL ztrmm(
'L',
'U', trans,
'N', k, n, one, t, ldt,
365 $ work, ldwork )
366
367 DO j = 1, n
368 DO i = 1, k
369 a( i, j ) = a( i, j ) - work( i, j )
370 END DO
371 END DO
372
373 CALL zgemm(
'N',
'N', m-l, n, k, -one, v, ldv, work, ldwork,
374 $ one, b, ldb )
375 CALL zgemm(
'N',
'N', l, n, k-l, -one, v( mp, kp ), ldv,
376 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
377 CALL ztrmm(
'L',
'U',
'N',
'N', l, n, one, v( mp, 1 ), ldv,
378 $ work, ldwork )
379 DO j = 1, n
380 DO i = 1, l
381 b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
382 END DO
383 END DO
384
385
386
387 ELSE IF( column .AND. forward .AND. right ) THEN
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403 np = min( n-l+1, n )
404 kp = min( l+1, k )
405
406 DO j = 1, l
407 DO i = 1, m
408 work( i, j ) = b( i, n-l+j )
409 END DO
410 END DO
411 CALL ztrmm(
'R',
'U',
'N',
'N', m, l, one, v( np, 1 ), ldv,
412 $ work, ldwork )
413 CALL zgemm(
'N',
'N', m, l, n-l, one, b, ldb,
414 $ v, ldv, one, work, ldwork )
415 CALL zgemm(
'N',
'N', m, k-l, n, one, b, ldb,
416 $ v( 1, kp ), ldv, zero, work( 1, kp ), ldwork )
417
418 DO j = 1, k
419 DO i = 1, m
420 work( i, j ) = work( i, j ) + a( i, j )
421 END DO
422 END DO
423
424 CALL ztrmm(
'R',
'U', trans,
'N', m, k, one, t, ldt,
425 $ work, ldwork )
426
427 DO j = 1, k
428 DO i = 1, m
429 a( i, j ) = a( i, j ) - work( i, j )
430 END DO
431 END DO
432
433 CALL zgemm(
'N',
'C', m, n-l, k, -one, work, ldwork,
434 $ v, ldv, one, b, ldb )
435 CALL zgemm(
'N',
'C', m, l, k-l, -one, work( 1, kp ), ldwork,
436 $ v( np, kp ), ldv, one, b( 1, np ), ldb )
437 CALL ztrmm(
'R',
'U',
'C',
'N', m, l, one, v( np, 1 ), ldv,
438 $ work, ldwork )
439 DO j = 1, l
440 DO i = 1, m
441 b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
442 END DO
443 END DO
444
445
446
447 ELSE IF( column .AND. backward .AND. left ) THEN
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464 mp = min( l+1, m )
465 kp = min( k-l+1, k )
466
467 DO j = 1, n
468 DO i = 1, l
469 work( k-l+i, j ) = b( i, j )
470 END DO
471 END DO
472
473 CALL ztrmm(
'L',
'L',
'C',
'N', l, n, one, v( 1, kp ), ldv,
474 $ work( kp, 1 ), ldwork )
475 CALL zgemm(
'C',
'N', l, n, m-l, one, v( mp, kp ), ldv,
476 $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
477 CALL zgemm(
'C',
'N', k-l, n, m, one, v, ldv,
478 $ b, ldb, zero, work, ldwork )
479
480 DO j = 1, n
481 DO i = 1, k
482 work( i, j ) = work( i, j ) + a( i, j )
483 END DO
484 END DO
485
486 CALL ztrmm(
'L',
'L', trans,
'N', k, n, one, t, ldt,
487 $ work, ldwork )
488
489 DO j = 1, n
490 DO i = 1, k
491 a( i, j ) = a( i, j ) - work( i, j )
492 END DO
493 END DO
494
495 CALL zgemm(
'N',
'N', m-l, n, k, -one, v( mp, 1 ), ldv,
496 $ work, ldwork, one, b( mp, 1 ), ldb )
497 CALL zgemm(
'N',
'N', l, n, k-l, -one, v, ldv,
498 $ work, ldwork, one, b, ldb )
499 CALL ztrmm(
'L',
'L',
'N',
'N', l, n, one, v( 1, kp ), ldv,
500 $ work( kp, 1 ), ldwork )
501 DO j = 1, n
502 DO i = 1, l
503 b( i, j ) = b( i, j ) - work( k-l+i, j )
504 END DO
505 END DO
506
507
508
509 ELSE IF( column .AND. backward .AND. right ) THEN
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525 np = min( l+1, n )
526 kp = min( k-l+1, k )
527
528 DO j = 1, l
529 DO i = 1, m
530 work( i, k-l+j ) = b( i, j )
531 END DO
532 END DO
533 CALL ztrmm(
'R',
'L',
'N',
'N', m, l, one, v( 1, kp ), ldv,
534 $ work( 1, kp ), ldwork )
535 CALL zgemm(
'N',
'N', m, l, n-l, one, b( 1, np ), ldb,
536 $ v( np, kp ), ldv, one, work( 1, kp ), ldwork )
537 CALL zgemm(
'N',
'N', m, k-l, n, one, b, ldb,
538 $ v, ldv, zero, work, ldwork )
539
540 DO j = 1, k
541 DO i = 1, m
542 work( i, j ) = work( i, j ) + a( i, j )
543 END DO
544 END DO
545
546 CALL ztrmm(
'R',
'L', trans,
'N', m, k, one, t, ldt,
547 $ work, ldwork )
548
549 DO j = 1, k
550 DO i = 1, m
551 a( i, j ) = a( i, j ) - work( i, j )
552 END DO
553 END DO
554
555 CALL zgemm(
'N',
'C', m, n-l, k, -one, work, ldwork,
556 $ v( np, 1 ), ldv, one, b( 1, np ), ldb )
557 CALL zgemm(
'N',
'C', m, l, k-l, -one, work, ldwork,
558 $ v, ldv, one, b, ldb )
559 CALL ztrmm(
'R',
'L',
'C',
'N', m, l, one, v( 1, kp ), ldv,
560 $ work( 1, kp ), ldwork )
561 DO j = 1, l
562 DO i = 1, m
563 b( i, j ) = b( i, j ) - work( i, k-l+j )
564 END DO
565 END DO
566
567
568
569 ELSE IF( row .AND. forward .AND. left ) THEN
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585 mp = min( m-l+1, m )
586 kp = min( l+1, k )
587
588 DO j = 1, n
589 DO i = 1, l
590 work( i, j ) = b( m-l+i, j )
591 END DO
592 END DO
593 CALL ztrmm(
'L',
'L',
'N',
'N', l, n, one, v( 1, mp ), ldv,
594 $ work, ldb )
595 CALL zgemm(
'N',
'N', l, n, m-l, one, v, ldv,b, ldb,
596 $ one, work, ldwork )
597 CALL zgemm(
'N',
'N', k-l, n, m, one, v( kp, 1 ), ldv,
598 $ b, ldb, zero, work( kp, 1 ), ldwork )
599
600 DO j = 1, n
601 DO i = 1, k
602 work( i, j ) = work( i, j ) + a( i, j )
603 END DO
604 END DO
605
606 CALL ztrmm(
'L',
'U', trans,
'N', k, n, one, t, ldt,
607 $ work, ldwork )
608
609 DO j = 1, n
610 DO i = 1, k
611 a( i, j ) = a( i, j ) - work( i, j )
612 END DO
613 END DO
614
615 CALL zgemm(
'C',
'N', m-l, n, k, -one, v, ldv, work, ldwork,
616 $ one, b, ldb )
617 CALL zgemm(
'C',
'N', l, n, k-l, -one, v( kp, mp ), ldv,
618 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
619 CALL ztrmm(
'L',
'L',
'C',
'N', l, n, one, v( 1, mp ), ldv,
620 $ work, ldwork )
621 DO j = 1, n
622 DO i = 1, l
623 b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
624 END DO
625 END DO
626
627
628
629 ELSE IF( row .AND. forward .AND. right ) THEN
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644 np = min( n-l+1, n )
645 kp = min( l+1, k )
646
647 DO j = 1, l
648 DO i = 1, m
649 work( i, j ) = b( i, n-l+j )
650 END DO
651 END DO
652 CALL ztrmm(
'R',
'L',
'C',
'N', m, l, one, v( 1, np ), ldv,
653 $ work, ldwork )
654 CALL zgemm(
'N',
'C', m, l, n-l, one, b, ldb, v, ldv,
655 $ one, work, ldwork )
656 CALL zgemm(
'N',
'C', m, k-l, n, one, b, ldb,
657 $ v( kp, 1 ), ldv, zero, work( 1, kp ), ldwork )
658
659 DO j = 1, k
660 DO i = 1, m
661 work( i, j ) = work( i, j ) + a( i, j )
662 END DO
663 END DO
664
665 CALL ztrmm(
'R',
'U', trans,
'N', m, k, one, t, ldt,
666 $ work, ldwork )
667
668 DO j = 1, k
669 DO i = 1, m
670 a( i, j ) = a( i, j ) - work( i, j )
671 END DO
672 END DO
673
674 CALL zgemm(
'N',
'N', m, n-l, k, -one, work, ldwork,
675 $ v, ldv, one, b, ldb )
676 CALL zgemm(
'N',
'N', m, l, k-l, -one, work( 1, kp ), ldwork,
677 $ v( kp, np ), ldv, one, b( 1, np ), ldb )
678 CALL ztrmm(
'R',
'L',
'N',
'N', m, l, one, v( 1, np ), ldv,
679 $ work, ldwork )
680 DO j = 1, l
681 DO i = 1, m
682 b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
683 END DO
684 END DO
685
686
687
688 ELSE IF( row .AND. backward .AND. left ) THEN
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704 mp = min( l+1, m )
705 kp = min( k-l+1, k )
706
707 DO j = 1, n
708 DO i = 1, l
709 work( k-l+i, j ) = b( i, j )
710 END DO
711 END DO
712 CALL ztrmm(
'L',
'U',
'N',
'N', l, n, one, v( kp, 1 ), ldv,
713 $ work( kp, 1 ), ldwork )
714 CALL zgemm(
'N',
'N', l, n, m-l, one, v( kp, mp ), ldv,
715 $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
716 CALL zgemm(
'N',
'N', k-l, n, m, one, v, ldv, b, ldb,
717 $ zero, work, ldwork )
718
719 DO j = 1, n
720 DO i = 1, k
721 work( i, j ) = work( i, j ) + a( i, j )
722 END DO
723 END DO
724
725 CALL ztrmm(
'L',
'L ', trans,
'N', k, n, one, t, ldt,
726 $ work, ldwork )
727
728 DO j = 1, n
729 DO i = 1, k
730 a( i, j ) = a( i, j ) - work( i, j )
731 END DO
732 END DO
733
734 CALL zgemm(
'C',
'N', m-l, n, k, -one, v( 1, mp ), ldv,
735 $ work, ldwork, one, b( mp, 1 ), ldb )
736 CALL zgemm(
'C',
'N', l, n, k-l, -one, v, ldv,
737 $ work, ldwork, one, b, ldb )
738 CALL ztrmm(
'L',
'U',
'C',
'N', l, n, one, v( kp, 1 ), ldv,
739 $ work( kp, 1 ), ldwork )
740 DO j = 1, n
741 DO i = 1, l
742 b( i, j ) = b( i, j ) - work( k-l+i, j )
743 END DO
744 END DO
745
746
747
748 ELSE IF( row .AND. backward .AND. right ) THEN
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763 np = min( l+1, n )
764 kp = min( k-l+1, k )
765
766 DO j = 1, l
767 DO i = 1, m
768 work( i, k-l+j ) = b( i, j )
769 END DO
770 END DO
771 CALL ztrmm(
'R',
'U',
'C',
'N', m, l, one, v( kp, 1 ), ldv,
772 $ work( 1, kp ), ldwork )
773 CALL zgemm(
'N',
'C', m, l, n-l, one, b( 1, np ), ldb,
774 $ v( kp, np ), ldv, one, work( 1, kp ), ldwork )
775 CALL zgemm(
'N',
'C', m, k-l, n, one, b, ldb, v, ldv,
776 $ zero, work, ldwork )
777
778 DO j = 1, k
779 DO i = 1, m
780 work( i, j ) = work( i, j ) + a( i, j )
781 END DO
782 END DO
783
784 CALL ztrmm(
'R',
'L', trans,
'N', m, k, one, t, ldt,
785 $ work, ldwork )
786
787 DO j = 1, k
788 DO i = 1, m
789 a( i, j ) = a( i, j ) - work( i, j )
790 END DO
791 END DO
792
793 CALL zgemm(
'N',
'N', m, n-l, k, -one, work, ldwork,
794 $ v( 1, np ), ldv, one, b( 1, np ), ldb )
795 CALL zgemm(
'N',
'N', m, l, k-l , -one, work, ldwork,
796 $ v, ldv, one, b, ldb )
797 CALL ztrmm(
'R',
'U',
'N',
'N', m, l, one, v( kp, 1 ), ldv,
798 $ work( 1, kp ), ldwork )
799 DO j = 1, l
800 DO i = 1, m
801 b( i, j ) = b( i, j ) - work( i, k-l+j )
802 END DO
803 END DO
804
805 END IF
806
807 RETURN
808
809
810
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
logical function lsame(ca, cb)
LSAME
subroutine ztrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRMM