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 A( LDA, * ), B( LDB, * ), T( LDT, * ),
262 $ V( LDV, * ), WORK( LDWORK, * )
263
264
265
266
267
268 COMPLEX 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 ctrmm(
'L',
'U',
'C',
'N', l, n, one, v( mp, 1 ), ldv,
352 $ work, ldwork )
353 CALL cgemm(
'C',
'N', l, n, m-l, one, v, ldv, b, ldb,
354 $ one, work, ldwork )
355 CALL cgemm(
'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 ctrmm(
'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 cgemm(
'N',
'N', m-l, n, k, -one, v, ldv, work, ldwork,
374 $ one, b, ldb )
375 CALL cgemm(
'N',
'N', l, n, k-l, -one, v( mp, kp ), ldv,
376 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
377 CALL ctrmm(
'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 ctrmm(
'R',
'U',
'N',
'N', m, l, one, v( np, 1 ), ldv,
412 $ work, ldwork )
413 CALL cgemm(
'N',
'N', m, l, n-l, one, b, ldb,
414 $ v, ldv, one, work, ldwork )
415 CALL cgemm(
'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 ctrmm(
'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 cgemm(
'N',
'C', m, n-l, k, -one, work, ldwork,
434 $ v, ldv, one, b, ldb )
435 CALL cgemm(
'N',
'C', m, l, k-l, -one, work( 1, kp ), ldwork,
436 $ v( np, kp ), ldv, one, b( 1, np ), ldb )
437 CALL ctrmm(
'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 ctrmm(
'L',
'L',
'C',
'N', l, n, one, v( 1, kp ), ldv,
474 $ work( kp, 1 ), ldwork )
475 CALL cgemm(
'C',
'N', l, n, m-l, one, v( mp, kp ), ldv,
476 $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
477 CALL cgemm(
'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 ctrmm(
'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 cgemm(
'N',
'N', m-l, n, k, -one, v( mp, 1 ), ldv,
496 $ work, ldwork, one, b( mp, 1 ), ldb )
497 CALL cgemm(
'N',
'N', l, n, k-l, -one, v, ldv,
498 $ work, ldwork, one, b, ldb )
499 CALL ctrmm(
'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 ctrmm(
'R',
'L',
'N',
'N', m, l, one, v( 1, kp ), ldv,
534 $ work( 1, kp ), ldwork )
535 CALL cgemm(
'N',
'N', m, l, n-l, one, b( 1, np ), ldb,
536 $ v( np, kp ), ldv, one, work( 1, kp ), ldwork )
537 CALL cgemm(
'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 ctrmm(
'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 cgemm(
'N',
'C', m, n-l, k, -one, work, ldwork,
556 $ v( np, 1 ), ldv, one, b( 1, np ), ldb )
557 CALL cgemm(
'N',
'C', m, l, k-l, -one, work, ldwork,
558 $ v, ldv, one, b, ldb )
559 CALL ctrmm(
'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 ctrmm(
'L',
'L',
'N',
'N', l, n, one, v( 1, mp ), ldv,
594 $ work, ldb )
595 CALL cgemm(
'N',
'N', l, n, m-l, one, v, ldv,b, ldb,
596 $ one, work, ldwork )
597 CALL cgemm(
'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 ctrmm(
'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 cgemm(
'C',
'N', m-l, n, k, -one, v, ldv, work, ldwork,
616 $ one, b, ldb )
617 CALL cgemm(
'C',
'N', l, n, k-l, -one, v( kp, mp ), ldv,
618 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
619 CALL ctrmm(
'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 ctrmm(
'R',
'L',
'C',
'N', m, l, one, v( 1, np ), ldv,
653 $ work, ldwork )
654 CALL cgemm(
'N',
'C', m, l, n-l, one, b, ldb, v, ldv,
655 $ one, work, ldwork )
656 CALL cgemm(
'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 ctrmm(
'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 cgemm(
'N',
'N', m, n-l, k, -one, work, ldwork,
675 $ v, ldv, one, b, ldb )
676 CALL cgemm(
'N',
'N', m, l, k-l, -one, work( 1, kp ), ldwork,
677 $ v( kp, np ), ldv, one, b( 1, np ), ldb )
678 CALL ctrmm(
'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 ctrmm(
'L',
'U',
'N',
'N', l, n, one, v( kp, 1 ), ldv,
713 $ work( kp, 1 ), ldwork )
714 CALL cgemm(
'N',
'N', l, n, m-l, one, v( kp, mp ), ldv,
715 $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
716 CALL cgemm(
'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 ctrmm(
'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 cgemm(
'C',
'N', m-l, n, k, -one, v( 1, mp ), ldv,
735 $ work, ldwork, one, b( mp, 1 ), ldb )
736 CALL cgemm(
'C',
'N', l, n, k-l, -one, v, ldv,
737 $ work, ldwork, one, b, ldb )
738 CALL ctrmm(
'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 ctrmm(
'R',
'U',
'C',
'N', m, l, one, v( kp, 1 ), ldv,
772 $ work( 1, kp ), ldwork )
773 CALL cgemm(
'N',
'C', m, l, n-l, one, b( 1, np ), ldb,
774 $ v( kp, np ), ldv, one, work( 1, kp ), ldwork )
775 CALL cgemm(
'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 ctrmm(
'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 cgemm(
'N',
'N', m, n-l, k, -one, work, ldwork,
794 $ v( 1, np ), ldv, one, b( 1, np ), ldb )
795 CALL cgemm(
'N',
'N', m, l, k-l , -one, work, ldwork,
796 $ v, ldv, one, b, ldb )
797 CALL ctrmm(
'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 cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
logical function lsame(ca, cb)
LSAME
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM