538{
539 double* buff_1 = FLA_DOUBLE_PTR(
FLA_ONE );
540 double* buff_0 = FLA_DOUBLE_PTR(
FLA_ZERO );
542
543 double alpha12;
544 double minus_conj_alpha12;
545 double psi11_minus_alpha12;
546 double minus_inv_tau11;
547 double minus_upsilon11;
548 double minus_conj_nu11;
549 double minus_conj_psi11;
550 double minus_zeta11;
551 double beta;
553
554
555 int b_alg = m_TS;
556
557
558
559
560
561
562
563
564 double* buff_w = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
565 double* buff_ap = (
double* )
FLA_malloc( n_A *
sizeof( *buff_A ) );
566 double* buff_u = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
567 double* buff_up = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
568 double* buff_v = (
double* )
FLA_malloc( n_A *
sizeof( *buff_A ) );
569 double* buff_y = (
double* )
FLA_malloc( n_A *
sizeof( *buff_A ) );
570 double* buff_z = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
571 int inc_w = 1;
572 int inc_ap = 1;
573 int inc_u = 1;
574 int inc_up = 1;
575 int inc_v = 1;
576 int inc_y = 1;
577 int inc_z = 1;
578
579 for (
i = 0;
i < b_alg; ++
i )
580 {
581 double* a10t = buff_A + (0 )*cs_A + (
i )*rs_A;
582 double* A20 = buff_A + (0 )*cs_A + (
i+1)*rs_A;
583 double* alpha11 = buff_A + (
i )*cs_A + (
i )*rs_A;
584 double* a21 = buff_A + (
i )*cs_A + (
i+1)*rs_A;
585 double* A02 = buff_A + (
i+1)*cs_A + (0 )*rs_A;
586 double* a12t = buff_A + (
i+1)*cs_A + (
i )*rs_A;
587 double* A22 = buff_A + (
i+1)*cs_A + (
i+1)*rs_A;
588
589 double* t01 = buff_T + (
i )*cs_T + (0 )*rs_T;
590 double* tau11 = buff_T + (
i )*cs_T + (
i )*rs_T;
591
592 double* s01 = buff_S + (
i )*cs_S + (0 )*rs_S;
593 double* sigma11 = buff_S + (
i )*cs_S + (
i )*rs_S;
594
595 double* w21 = buff_w + (
i+1)*inc_w;
596
597 double* a12p = buff_ap + (
i+1)*inc_ap;
598
599 double* upsilon11 = buff_u + (
i )*inc_u;
600 double* u21 = buff_u + (
i+1)*inc_u;
601
602 double* u21p = buff_up + (
i+1)*inc_up;
603
604 double* nu11 = buff_v + (
i )*inc_v;
605 double* v21 = buff_v + (
i+1)*inc_v;
606
607 double* psi11 = buff_y + (
i )*inc_y;
608 double* y21 = buff_y + (
i+1)*inc_y;
609
610 double* zeta11 = buff_z + (
i )*inc_z;
611 double* z21 = buff_z + (
i+1)*inc_z;
612
613 double* a12p_t = a12p + (0 )*inc_ap;
614 double* a12p_b = a12p + (1 )*inc_ap;
615
616 double* v21_t = v21 + (0 )*inc_v;
617 double* v21_b = v21 + (1 )*inc_v;
618
619 double* a12t_l = a12t + (0 )*cs_A + (0 )*rs_A;
620 double* a12t_r = a12t + (1 )*cs_A + (0 )*rs_A;
621
622 double* A22_l = A22 + (0 )*cs_A + (0 )*rs_A;
623
624 int m_ahead = m_A -
i - 1;
625 int n_ahead = n_A -
i - 1;
628
629
630
631 if ( m_behind > 0 )
632 {
633
634
635 bl1_dmult3( buff_m1, upsilon11, &minus_upsilon11 );
636
637
638
639 bl1_dmult3( buff_m1, zeta11, &minus_zeta11 );
640
641
642
643 bl1_dcopyconj( psi11, &minus_conj_psi11 );
644 bl1_dscals( buff_m1, &minus_conj_psi11 );
645
646
647
648 bl1_dcopyconj( nu11, &minus_conj_nu11 );
649 bl1_dscals( buff_m1, &minus_conj_nu11 );
650
651
652
654 1,
655 &minus_conj_psi11,
656 upsilon11, 1,
657 alpha11, 1 );
659 1,
660 &minus_conj_nu11,
661 zeta11, 1,
662 alpha11, 1 );
663
664
665
667 m_ahead,
668 &minus_conj_psi11,
669 u21, inc_u,
670 a21, rs_A );
672 m_ahead,
673 &minus_conj_nu11,
674 z21, inc_z,
675 a21, rs_A );
676
677
678
680 n_ahead,
681 &minus_upsilon11,
682 y21, inc_y,
683 a12t, cs_A );
685 n_ahead,
686 &minus_zeta11,
687 v21, inc_v,
688 a12t, cs_A );
689 }
690
691
692
693
694
696 alpha11,
697 a21, rs_A,
698 tau11 );
700 m_ahead,
701 a21, rs_A,
702 u21p, inc_up );
703
704 if ( n_ahead > 0 )
705 {
706
707
708 bl1_ddiv3( buff_m1, tau11, &minus_inv_tau11 );
709
710
711
713 n_ahead,
714 a12t, cs_A,
715 a12p, inc_ap );
717 n_ahead,
718 &minus_inv_tau11,
719 a12t, cs_A,
720 a12p, inc_ap );
721 }
722
723 if ( m_behind > 0 )
724 {
725
726
729 m_ahead,
730 n_ahead,
731 buff_m1,
732 u21, inc_u,
733 y21, inc_y,
734 A22, rs_A, cs_A );
737 m_ahead,
738 n_ahead,
739 buff_m1,
740 z21, inc_z,
741 v21, inc_v,
742 A22, rs_A, cs_A );
743 }
744
745 if ( n_ahead > 0 )
746 {
747
750 m_ahead,
751 n_ahead,
752 buff_1,
753 A22, rs_A, cs_A,
754 u21p, inc_up,
755 buff_0,
756 y21, inc_y );
757
758
760 n_ahead,
761 &minus_inv_tau11,
762 y21, inc_y,
763 a12p, inc_ap );
764
765
768 m_ahead,
769 n_ahead,
770 buff_1,
771 A22, rs_A, cs_A,
772 a12p, inc_ap,
773 buff_0,
774 w21, inc_w );
775
776
778 n_ahead,
779 buff_1,
780 a12t, cs_A,
781 y21, inc_y );
782
783
784
785
786
788 a12p_t,
789 a12p_b, inc_ap,
790 &alpha12,
791 &psi11_minus_alpha12,
792 sigma11 );
793
794
795
796
797
799 n_ahead,
800 a12p, inc_ap,
801 v21, inc_v );
802 bl1_dmult4( buff_m1, &alpha12, v21_t, v21_t );
804 n_ahead,
805 &psi11_minus_alpha12,
806 v21, inc_v );
808 v21_b, inc_v );
809
810
811
812 *a12t_l = alpha12;
814 n_ahead - 1,
815 v21_b, inc_v,
816 a12t_r, cs_A );
817 }
818
819
821 m_ahead,
822 u21p, inc_up,
823 u21, inc_u );
824
825 if ( n_ahead > 0 )
826 {
827
828
829
831 n_ahead,
832 y21, inc_y,
833 v21, inc_v,
834 &beta );
835 bl1_dscals( &minus_inv_tau11, &beta );
836
837
838
839 bl1_dcopyconj( &alpha12, &minus_conj_alpha12 );
840 bl1_dneg1( &minus_conj_alpha12 );
841
842
843
844
845
847 m_ahead,
848 w21, inc_w,
849 z21, inc_z );
851 m_ahead,
852 &minus_conj_alpha12,
853 A22_l, rs_A,
854 z21, inc_z );
856 m_ahead,
857 &psi11_minus_alpha12,
858 z21, inc_z );
860 m_ahead,
861 &beta,
862 u21, inc_u,
863 z21, inc_z );
864
865
866
868 n_ahead,
869 tau11,
870 y21, inc_y );
872 m_ahead,
873 sigma11,
874 z21, inc_z );
875
876
879 m_behind,
880 n_ahead,
881 buff_1,
882 A02, rs_A, cs_A,
883 v21, inc_v,
884 buff_0,
885 s01, rs_S );
886 }
887
888
889
891 n_behind,
892 a10t, cs_A,
893 t01, rs_T );
896 m_ahead,
897 n_behind,
898 buff_1,
899 A20, rs_A, cs_A,
900 u21, inc_u,
901 buff_1,
902 t01, rs_T );
903
904 if ( m_behind + 1 == b_alg && n_ahead > 0 )
905 {
906
907
910 m_ahead,
911 n_ahead,
912 buff_m1,
913 u21, inc_u,
914 y21, inc_y,
915 A22, rs_A, cs_A );
918 m_ahead,
919 n_ahead,
920 buff_m1,
921 z21, inc_z,
922 v21, inc_v,
923 A22, rs_A, cs_A );
924 }
925
926
927
928 }
929
930
931
932
933
934
935
936
944
945 return FLA_SUCCESS;
946}
FLA_Error FLA_Househ2_UT_l_opd(int m_x2, double *chi_1, double *x2, int inc_x2, double *tau)
Definition FLA_Househ2_UT.c:274
FLA_Error FLA_Househ2s_UT_r_opd(int m_x2, double *chi_1, double *x2, int inc_x2, double *alpha, double *chi_1_minus_alpha, double *tau)
Definition FLA_Househ2s_UT.c:572
void bl1_daxpyv(conj1_t conj, int n, double *alpha, double *x, int incx, double *y, int incy)
Definition bl1_axpyv.c:21
void bl1_dconjv(int m, double *x, int incx)
Definition bl1_conjv.c:18
void bl1_dcopyv(conj1_t conj, int m, double *x, int incx, double *y, int incy)
Definition bl1_copyv.c:42
void bl1_ddot(conj1_t conj, int n, double *x, int incx, double *y, int incy, double *rho)
Definition bl1_dot.c:26
void bl1_dgemv(trans1_t transa, conj1_t conjx, int m, int n, double *alpha, double *a, int a_rs, int a_cs, double *x, int incx, double *beta, double *y, int incy)
Definition bl1_gemv.c:69
void bl1_dger(conj1_t conjx, conj1_t conjy, int m, int n, double *alpha, double *x, int incx, double *y, int incy, double *a, int a_rs, int a_cs)
Definition bl1_ger.c:62
void bl1_dinvscalv(conj1_t conj, int n, double *alpha, double *x, int incx)
Definition bl1_invscalv.c:26