711{
712 double* buff_1 = FLA_DOUBLE_PTR(
FLA_ONE );
713 double* buff_0 = FLA_DOUBLE_PTR(
FLA_ZERO );
715
716 double alpha12;
717 double minus_conj_alpha12;
718 double psi11_minus_alpha12;
719 double minus_inv_tau11;
720 double beta;
721 double last_elem;
723
724
725 int b_alg = m_TS;
726
727
728
729
730
731
732
733
734
735
736
737 double* buff_w = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
738 double* buff_al = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
739 double* buff_ap = (
double* )
FLA_malloc( n_A *
sizeof( *buff_A ) );
740 double* buff_u = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
741 double* buff_up = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
742 double* buff_v = (
double* )
FLA_malloc( n_A *
sizeof( *buff_A ) );
743 double* buff_d = (
double* )
FLA_malloc( n_A *
sizeof( *buff_A ) );
744 double* buff_e = (
double* )
FLA_malloc( n_A *
sizeof( *buff_A ) );
745 double* buff_f = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
746 double* buff_g = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
747 int inc_w = 1;
748 int inc_al = 1;
749 int inc_ap = 1;
750 int inc_u = 1;
751 int inc_up = 1;
752 int inc_v = 1;
753 int inc_d = 1;
754 int inc_e = 1;
755 int inc_f = 1;
756 int inc_g = 1;
757
758
759
761 b_alg,
762 buff_0,
763 buff_Y, rs_Y, cs_Y );
765 b_alg,
766 buff_0,
767 buff_Z, rs_Z, cs_Z );
768
769 for (
i = 0;
i < b_alg; ++
i )
770 {
771 double* a10t = buff_A + (0 )*cs_A + (
i )*rs_A;
772 double* A20 = buff_A + (0 )*cs_A + (
i+1)*rs_A;
773 double* a01 = buff_A + (
i )*cs_A + (0 )*rs_A;
774 double* alpha11 = buff_A + (
i )*cs_A + (
i )*rs_A;
775 double* a21 = buff_A + (
i )*cs_A + (
i+1)*rs_A;
776 double* A02 = buff_A + (
i+1)*cs_A + (0 )*rs_A;
777 double* a12t = buff_A + (
i+1)*cs_A + (
i )*rs_A;
778 double* A22 = buff_A + (
i+1)*cs_A + (
i+1)*rs_A;
779
780 double* y10t = buff_Y + (0 )*cs_Y + (
i )*rs_Y;
781 double* Y20 = buff_Y + (0 )*cs_Y + (
i+1)*rs_Y;
782 double* y21 = buff_Y + (
i )*cs_Y + (
i+1)*rs_Y;
783
784 double* z10t = buff_Z + (0 )*cs_Z + (
i )*rs_Z;
785 double* Z20 = buff_Z + (0 )*cs_Z + (
i+1)*rs_Z;
786 double* z21 = buff_Z + (
i )*cs_Z + (
i+1)*rs_Z;
787
788 double* t01 = buff_T + (
i )*cs_T + (0 )*rs_T;
789 double* tau11 = buff_T + (
i )*cs_T + (
i )*rs_T;
790
791 double* s01 = buff_S + (
i )*cs_S + (0 )*rs_S;
792 double* sigma11 = buff_S + (
i )*cs_S + (
i )*rs_S;
793
794 double* w21 = buff_w + (
i+1)*inc_w;
795
796 double* a22l = buff_al + (
i+1)*inc_al;
797
798 double* a12p = buff_ap + (
i+1)*inc_ap;
799
800 double* u21 = buff_u + (
i+1)*inc_u;
801
802 double* u21p = buff_up + (
i+1)*inc_up;
803
804 double* v21 = buff_v + (
i+1)*inc_v;
805
806 double* d0 = buff_d + (0 )*inc_d;
807
808 double* e0 = buff_e + (0 )*inc_e;
809
810 double* f0 = buff_f + (0 )*inc_f;
811
812 double* g0 = buff_g + (0 )*inc_g;
813
814 double* a12p_t = a12p + (0 )*inc_ap;
815 double* a12p_b = a12p + (1 )*inc_ap;
816
817 double* v21_t = v21 + (0 )*inc_v;
818 double* v21_b = v21 + (1 )*inc_v;
819
820 double* a01_b = a01 + (0 )*cs_A + (
i-1)*rs_A;
821
822 double* a12t_l = a12t + (0 )*cs_A + (0 )*rs_A;
823 double* a12t_r = a12t + (1 )*cs_A + (0 )*rs_A;
824
825 double* A02_l = A02 + (0 )*cs_A + (0 )*rs_A;
826
827 double* A22_l = A22 + (0 )*cs_A + (0 )*rs_A;
828
829 double* Y20_t = Y20 + (0 )*cs_Y + (0 )*rs_Y;
830
831 double* ABL = a10t;
832 double* ZBL = z10t;
833
834 double* a2 = alpha11;
835
836 int m_ahead = m_A -
i - 1;
837 int n_ahead = n_A -
i - 1;
840
841
842
843 if ( m_behind > 0 )
844 {
845
846
847 last_elem = *a01_b;
848 *a01_b = *buff_1;
849 }
850
851
852
855 m_ahead + 1,
856 n_behind,
857 buff_m1,
858 ABL, rs_A, cs_A,
859 y10t, cs_Y,
860 buff_1,
861 a2, rs_A );
864 m_ahead + 1,
865 n_behind,
866 buff_m1,
867 ZBL, rs_Z, cs_Z,
868 a01, rs_A,
869 buff_1,
870 a2, rs_A );
871
872
873
876 n_ahead,
877 n_behind,
878 buff_m1,
879 Y20, rs_Y, cs_Y,
880 a10t, cs_A,
881 buff_1,
882 a12t, cs_A );
885 m_behind,
886 n_ahead,
887 buff_m1,
888 A02, rs_A, cs_A,
889 z10t, cs_Z,
890 buff_1,
891 a12t, cs_A );
892
893 if ( m_behind > 0 )
894 {
895
896 *a01_b = last_elem;
897 }
898
899
900
901
902
904 alpha11,
905 a21, rs_A,
906 tau11 );
908 m_ahead,
909 a21, rs_A,
910 u21p, inc_up );
911
912 if ( n_ahead > 0 )
913 {
914
915
916 bl1_ddiv3( buff_m1, tau11, &minus_inv_tau11 );
917
918
919
921 n_ahead,
922 a12t, cs_A,
923 a12p, inc_ap );
925 n_ahead,
926 &minus_inv_tau11,
927 a12t, cs_A,
928 a12p, inc_ap );
929
930
931
934 m_ahead,
935 n_behind,
936 buff_1,
937 A20, rs_A, cs_A,
938 u21p, inc_up,
939 buff_0,
940 d0, inc_d );
943 m_ahead,
944 n_behind,
945 buff_1,
946 Z20, rs_Z, cs_Z,
947 u21p, inc_up,
948 buff_0,
949 e0, inc_e );
950
951
952
954 n_behind,
955 a10t, cs_A,
956 t01, rs_T );
958 n_behind,
959 buff_1,
960 d0, inc_d,
961 t01, rs_T );
962
963
964
965
967 buff_0,
968 y21, rs_Y );
971 n_ahead,
972 n_behind,
973 buff_m1,
974 Y20, rs_Y, cs_Y,
975 d0, inc_d,
976 buff_1,
977 y21, rs_Y );
980 m_behind,
981 n_ahead,
982 buff_m1,
983 A02, rs_A, cs_A,
984 e0, inc_e,
985 buff_1,
986 y21, rs_Y );
987
988
991 m_ahead,
992 n_ahead,
993 buff_1,
994 A22, rs_A, cs_A,
995 u21p, inc_up,
996 buff_1,
997 y21, rs_Y );
998
999
1001 n_ahead,
1002 &minus_inv_tau11,
1003 y21, rs_Y,
1004 a12p, inc_ap );
1005
1006
1009 m_ahead,
1010 n_ahead,
1011 buff_1,
1012 A22, rs_A, cs_A,
1013 a12p, inc_ap,
1014 buff_0,
1015 w21, inc_w );
1016
1017
1018
1021 n_ahead,
1022 n_behind,
1023 buff_1,
1024 Y20, rs_Y, cs_Y,
1025 a12p, inc_ap,
1026 buff_0,
1027 f0, inc_f );
1030 m_behind,
1031 n_ahead,
1032 buff_1,
1033 A02, rs_A, cs_A,
1034 a12p, inc_ap,
1035 buff_0,
1036 g0, inc_g );
1037
1038
1039
1042 m_ahead,
1043 n_behind,
1044 buff_m1,
1045 A20, rs_A, cs_A,
1046 f0, inc_f,
1047 buff_1,
1048 w21, inc_w );
1051 m_ahead,
1052 n_behind,
1053 buff_m1,
1054 Z20, rs_Z, cs_Z,
1055 g0, inc_g,
1056 buff_1,
1057 w21, inc_w );
1058
1059
1060
1061
1063 m_ahead,
1064 A22_l, rs_A,
1065 a22l, inc_al );
1068 m_ahead,
1069 n_behind,
1070 buff_m1,
1071 A20, rs_A, cs_A,
1072 Y20_t, cs_Y,
1073 buff_1,
1074 a22l, inc_al );
1077 m_ahead,
1078 n_behind,
1079 buff_m1,
1080 Z20, rs_Z, cs_Z,
1081 A02_l, rs_A,
1082 buff_1,
1083 a22l, inc_al );
1084
1085
1087 n_ahead,
1088 buff_1,
1089 a12t, cs_A,
1090 y21, rs_Y );
1091
1092
1093
1094
1095
1097 a12p_t,
1098 a12p_b, inc_ap,
1099 &alpha12,
1100 &psi11_minus_alpha12,
1101 sigma11 );
1102
1103
1104
1105
1106
1108 n_ahead,
1109 a12p, inc_ap,
1110 v21, inc_v );
1111 bl1_dmult4( buff_m1, &alpha12, v21_t, v21_t );
1113 n_ahead,
1114 &psi11_minus_alpha12,
1115 v21, inc_v );
1117 v21_b, inc_v );
1118
1119
1120
1121 bl1_dcopyconj( &alpha12, &minus_conj_alpha12 );
1122 bl1_dneg1( &minus_conj_alpha12 );
1123
1124
1125
1126
1128 n_behind,
1129 g0, inc_g,
1130 s01, rs_S );
1132 n_behind,
1133 &minus_conj_alpha12,
1134 A02_l, rs_A,
1135 s01, rs_S );
1137 n_behind,
1138 &psi11_minus_alpha12,
1139 s01, rs_S );
1140
1141
1142
1143 *a12t_l = alpha12;
1145 n_ahead - 1,
1146 v21_b, inc_v,
1147 a12t_r, cs_A );
1148 }
1149
1150
1152 m_ahead,
1153 u21p, inc_up,
1154 u21, inc_u );
1155
1156 if ( n_ahead > 0 )
1157 {
1158
1159
1160
1162 n_ahead,
1163 y21, rs_Y,
1164 v21, inc_v,
1165 &beta );
1166 bl1_dscals( &minus_inv_tau11, &beta );
1167
1168
1169
1170
1171
1173 m_ahead,
1174 w21, inc_w,
1175 z21, rs_Z );
1177 m_ahead,
1178 &minus_conj_alpha12,
1179 a22l, inc_al,
1180 z21, rs_Z );
1182 m_ahead,
1183 &psi11_minus_alpha12,
1184 z21, rs_Z );
1186 m_ahead,
1187 &beta,
1188 u21, inc_u,
1189 z21, rs_Z );
1190
1191
1192
1194 n_ahead,
1195 tau11,
1196 y21, rs_Y );
1198 m_ahead,
1199 sigma11,
1200 z21, rs_Z );
1201 }
1202 else
1203 {
1204
1205
1207 n_behind,
1208 a10t, cs_A,
1209 t01, rs_T );
1212 m_ahead,
1213 n_behind,
1214 buff_1,
1215 A20, rs_A, cs_A,
1216 u21, inc_u,
1217 buff_1,
1218 t01, rs_T );
1219 }
1220
1221
1222
1223 }
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1245
1246 return FLA_SUCCESS;
1247}
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_dinvscalv(conj1_t conj, int n, double *alpha, double *x, int incx)
Definition bl1_invscalv.c:26
void bl1_dsetv(int m, double *sigma, double *x, int incx)
Definition bl1_setv.c:39
void bl1_dsetm(int m, int n, double *sigma, double *a, int a_rs, int a_cs)
Definition bl1_setm.c:45