1 /*
2
3 Copyright (C) 2014, The University of Texas at Austin
4
5 This file is part of libflame and is available under the 3-Clause
6 BSD license, which can be found in the LICENSE file at the top-level
7 directory, or at http://opensource.org/licenses/BSD-3-Clause
8
9 */
10
11 #include "FLAME.h"
12
FLA_Bidiag_UT_u_unb_var2(FLA_Obj A,FLA_Obj TU,FLA_Obj TV)13 FLA_Error FLA_Bidiag_UT_u_unb_var2( FLA_Obj A, FLA_Obj TU, FLA_Obj TV )
14 {
15 return FLA_Bidiag_UT_u_step_unb_var2( A, TU, TV );
16 }
17
FLA_Bidiag_UT_u_step_unb_var2(FLA_Obj A,FLA_Obj T,FLA_Obj S)18 FLA_Error FLA_Bidiag_UT_u_step_unb_var2( FLA_Obj A, FLA_Obj T, FLA_Obj S )
19 {
20 FLA_Obj ATL, ATR, A00, a01, A02,
21 ABL, ABR, a10t, alpha11, a12t,
22 A20, a21, A22;
23 FLA_Obj TTL, TTR, T00, t01, T02,
24 TBL, TBR, t10t, tau11, t12t,
25 T20, t21, T22;
26 FLA_Obj STL, STR, S00, s01, S02,
27 SBL, SBR, s10t, sigma11, s12t,
28 S20, s21, S22;
29 FLA_Obj yT, y01,
30 yB, psi11,
31 y21;
32 FLA_Obj zT, z01,
33 zB, zeta11,
34 z21;
35 FLA_Obj vT, v01,
36 vB, nu11,
37 v21;
38 FLA_Obj v, y, z;
39
40 FLA_Obj beta;
41
42 FLA_Obj a12t_l, a12t_r;
43 FLA_Obj v21_t,
44 v21_b;
45
46 FLA_Datatype datatype_A;
47 dim_t m_A, n_A;
48 dim_t b_alg;
49
50
51 b_alg = FLA_Obj_length( T );
52
53 datatype_A = FLA_Obj_datatype( A );
54 m_A = FLA_Obj_length( A );
55 n_A = FLA_Obj_width( A );
56
57 FLA_Obj_create( datatype_A, 1, 1, 0, 0, &beta );
58 FLA_Obj_create( datatype_A, n_A, 1, 0, 0, &v );
59 FLA_Obj_create( datatype_A, n_A, 1, 0, 0, &y );
60 FLA_Obj_create( datatype_A, m_A, 1, 0, 0, &z );
61
62 FLA_Part_2x2( A, &ATL, &ATR,
63 &ABL, &ABR, 0, 0, FLA_TL );
64 FLA_Part_2x2( T, &TTL, &TTR,
65 &TBL, &TBR, 0, 0, FLA_TL );
66 FLA_Part_2x2( S, &STL, &STR,
67 &SBL, &SBR, 0, 0, FLA_TL );
68 FLA_Part_2x1( v, &vT,
69 &vB, 0, FLA_TOP );
70 FLA_Part_2x1( y, &yT,
71 &yB, 0, FLA_TOP );
72 FLA_Part_2x1( z, &zT,
73 &zB, 0, FLA_TOP );
74
75 while ( FLA_Obj_length( ATL ) < b_alg )
76 {
77 FLA_Repart_2x2_to_3x3( ATL, /**/ ATR, &A00, /**/ &a01, &A02,
78 /* ************* */ /* ************************** */
79 &a10t, /**/ &alpha11, &a12t,
80 ABL, /**/ ABR, &A20, /**/ &a21, &A22,
81 1, 1, FLA_BR );
82 FLA_Repart_2x2_to_3x3( TTL, /**/ TTR, &T00, /**/ &t01, &T02,
83 /* ************* */ /* ************************** */
84 &t10t, /**/ &tau11, &t12t,
85 TBL, /**/ TBR, &T20, /**/ &t21, &T22,
86 1, 1, FLA_BR );
87 FLA_Repart_2x2_to_3x3( STL, /**/ STR, &S00, /**/ &s01, &S02,
88 /* ************* */ /* ************************** */
89 &s10t, /**/ &sigma11, &s12t,
90 SBL, /**/ SBR, &S20, /**/ &s21, &S22,
91 1, 1, FLA_BR );
92 FLA_Repart_2x1_to_3x1( vT, &v01,
93 /* ** */ /* ***** */
94 &nu11,
95 vB, &v21, 1, FLA_BOTTOM );
96 FLA_Repart_2x1_to_3x1( yT, &y01,
97 /* ** */ /* ***** */
98 &psi11,
99 yB, &y21, 1, FLA_BOTTOM );
100 FLA_Repart_2x1_to_3x1( zT, &z01,
101 /* ** */ /* ***** */
102 &zeta11,
103 zB, &z21, 1, FLA_BOTTOM );
104
105 /*------------------------------------------------------------*/
106
107 // [ alpha11_new, u21, tau11 ] = House2( alpha11, a21 );
108 FLA_Househ2_UT( FLA_LEFT,
109 alpha11,
110 a21, tau11 );
111
112 if ( FLA_Obj_width( A22 ) > 0 )
113 {
114 // y21' = a12t + u21' * A22;
115 // y21 = conj(a12t) + A22' * u21;
116 FLA_Copyt( FLA_CONJ_TRANSPOSE, a12t, y21 );
117 FLA_Gemvc( FLA_CONJ_TRANSPOSE, FLA_NO_CONJUGATE, FLA_ONE, A22, a21, FLA_ONE, y21 );
118
119 // y21 = y21 / tau11;
120 FLA_Inv_scalc( FLA_NO_CONJUGATE, tau11, y21 );
121
122 // a12t = a12t - conj(y21)^T;
123 FLA_Axpyt( FLA_CONJ_TRANSPOSE, FLA_MINUS_ONE, y21, a12t );
124
125 FLA_Part_1x2( a12t, &a12t_l, &a12t_r, 1, FLA_LEFT );
126 FLA_Part_2x1( v21, &v21_t,
127 &v21_b, 1, FLA_TOP );
128
129 // [ a12t_l, v12t_b, sigma11 ] = House2( a12t_l, a12t_r );
130 FLA_Househ2_UT( FLA_RIGHT, a12t_l, a12t_r, sigma11 );
131
132 // v21_t = 1;
133 // v21_b = a12t_r^T;
134 FLA_Set( FLA_ONE, v21_t );
135 FLA_Copyt( FLA_TRANSPOSE, a12t_r, v21_b );
136
137 // beta = - y21' * v21;
138 FLA_Dotc( FLA_CONJUGATE, y21, v21, beta );
139 FLA_Scal( FLA_MINUS_ONE, beta );
140
141 // z21 = ( A22 - u21 * y21' ) * v21 / sigma11;
142 // = ( A22 * v21 - u21 * y21' * v21 ) / sigma11;
143 // = ( A22 * v21 + beta * u21 ) / sigma11;
144 FLA_Copy( a21, z21 );
145 FLA_Gemvc( FLA_NO_TRANSPOSE, FLA_NO_CONJUGATE, FLA_ONE, A22, v21, beta, z21 );
146 FLA_Inv_scalc( FLA_NO_CONJUGATE, sigma11, z21 );
147
148 // A22 = A22 - u21 * y21' - z21 * v21';
149 FLA_Gerc( FLA_NO_CONJUGATE, FLA_CONJUGATE, FLA_MINUS_ONE, a21, y21, A22 );
150 FLA_Gerc( FLA_NO_CONJUGATE, FLA_CONJUGATE, FLA_MINUS_ONE, z21, v21, A22 );
151
152 // s01 = conj(V02) * v21;
153 FLA_Gemv( FLA_CONJ_NO_TRANSPOSE, FLA_ONE, A02, v21, FLA_ZERO, s01 );
154 }
155
156 // t01 = a10t' + U20' * u21;
157 FLA_Copyt( FLA_CONJ_TRANSPOSE, a10t, t01 );
158 FLA_Gemv( FLA_CONJ_TRANSPOSE, FLA_ONE, A20, a21, FLA_ONE, t01 );
159
160 /*------------------------------------------------------------*/
161
162 FLA_Cont_with_3x3_to_2x2( &ATL, /**/ &ATR, A00, a01, /**/ A02,
163 a10t, alpha11, /**/ a12t,
164 /* ************** */ /* ************************ */
165 &ABL, /**/ &ABR, A20, a21, /**/ A22,
166 FLA_TL );
167 FLA_Cont_with_3x3_to_2x2( &TTL, /**/ &TTR, T00, t01, /**/ T02,
168 t10t, tau11, /**/ t12t,
169 /* ************** */ /* ************************ */
170 &TBL, /**/ &TBR, T20, t21, /**/ T22,
171 FLA_TL );
172 FLA_Cont_with_3x3_to_2x2( &STL, /**/ &STR, S00, s01, /**/ S02,
173 s10t, sigma11, /**/ s12t,
174 /* ************** */ /* ************************ */
175 &SBL, /**/ &SBR, S20, s21, /**/ S22,
176 FLA_TL );
177 FLA_Cont_with_3x1_to_2x1( &vT, v01,
178 nu11,
179 /* ** */ /* ***** */
180 &vB, v21, FLA_TOP );
181 FLA_Cont_with_3x1_to_2x1( &yT, y01,
182 psi11,
183 /* ** */ /* ***** */
184 &yB, y21, FLA_TOP );
185 FLA_Cont_with_3x1_to_2x1( &zT, z01,
186 zeta11,
187 /* ** */ /* ***** */
188 &zB, z21, FLA_TOP );
189 }
190
191 FLA_Obj_free( &beta );
192 FLA_Obj_free( &v );
193 FLA_Obj_free( &y );
194 FLA_Obj_free( &z );
195
196 return FLA_SUCCESS;
197 }
198
199