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_Eig_gest_iu_unb_var1(FLA_Obj A,FLA_Obj Y,FLA_Obj B)13 FLA_Error FLA_Eig_gest_iu_unb_var1( FLA_Obj A, FLA_Obj Y, FLA_Obj B )
14 {
15 FLA_Obj ATL, ATR, A00, a01, A02,
16 ABL, ABR, a10t, alpha11, a12t,
17 A20, a21, A22;
18
19 FLA_Obj BTL, BTR, B00, b01, B02,
20 BBL, BBR, b10t, beta11, b12t,
21 B20, b21, B22;
22
23 FLA_Obj yT, y01,
24 yB, psi11,
25 y21;
26
27 FLA_Obj y01_l, y01_r;
28
29 FLA_Part_2x2( A, &ATL, &ATR,
30 &ABL, &ABR, 0, 0, FLA_TL );
31
32 FLA_Part_2x2( B, &BTL, &BTR,
33 &BBL, &BBR, 0, 0, FLA_TL );
34
35 FLA_Part_2x1( Y, &yT,
36 &yB, 0, FLA_TOP );
37
38 while ( FLA_Obj_length( ATL ) < FLA_Obj_length( A ) ){
39
40 FLA_Repart_2x2_to_3x3( ATL, /**/ ATR, &A00, /**/ &a01, &A02,
41 /* ************* */ /* ************************** */
42 &a10t, /**/ &alpha11, &a12t,
43 ABL, /**/ ABR, &A20, /**/ &a21, &A22,
44 1, 1, FLA_BR );
45
46 FLA_Repart_2x2_to_3x3( BTL, /**/ BTR, &B00, /**/ &b01, &B02,
47 /* ************* */ /* ************************* */
48 &b10t, /**/ &beta11, &b12t,
49 BBL, /**/ BBR, &B20, /**/ &b21, &B22,
50 1, 1, FLA_BR );
51
52 FLA_Repart_2x1_to_3x1( yT, &y01,
53 /* ** */ /* ***** */
54 &psi11,
55 yB, &y21, 1, FLA_BOTTOM );
56
57 /*------------------------------------------------------------*/
58
59 FLA_Part_1x2( y01, &y01_l, &y01_r, 1, FLA_LEFT );
60
61 // y01 = A00 * b01;
62 FLA_Hemvc_external( FLA_UPPER_TRIANGULAR, FLA_NO_CONJUGATE,
63 FLA_ONE, A00, b01, FLA_ZERO, y01_l );
64
65 // a01 = inv( triu( B00 )' ) * a01;
66 FLA_Trsv_external( FLA_UPPER_TRIANGULAR, FLA_CONJ_TRANSPOSE, FLA_NONUNIT_DIAG,
67 B00, a01 );
68
69 // a01 = a01 - 1/2 * y01;
70 FLA_Axpy_external( FLA_MINUS_ONE_HALF, y01_l, a01 );
71
72 // alpha11 = alpha11 - a01' * b01 - b01' * a01;
73 FLA_Dot2cs_external( FLA_CONJUGATE, FLA_MINUS_ONE, a01, b01, FLA_ONE, alpha11 );
74
75 // alpha11 = inv(beta11) * alpha11 * inv(conj(beta11));
76 // = inv(beta11) * alpha11 * inv(beta11);
77 FLA_Inv_scal_external( beta11, alpha11 );
78 FLA_Inv_scal_external( beta11, alpha11 );
79
80 // a01 = a01 - 1/2 * y01;
81 FLA_Axpy_external( FLA_MINUS_ONE_HALF, y01_l, a01 );
82
83 // a01 = a01 * inv(beta11);
84 FLA_Inv_scal_external( beta11, a01 );
85
86 /*------------------------------------------------------------*/
87
88 FLA_Cont_with_3x3_to_2x2( &ATL, /**/ &ATR, A00, a01, /**/ A02,
89 a10t, alpha11, /**/ a12t,
90 /* ************** */ /* ************************ */
91 &ABL, /**/ &ABR, A20, a21, /**/ A22,
92 FLA_TL );
93
94 FLA_Cont_with_3x3_to_2x2( &BTL, /**/ &BTR, B00, b01, /**/ B02,
95 b10t, beta11, /**/ b12t,
96 /* ************** */ /* *********************** */
97 &BBL, /**/ &BBR, B20, b21, /**/ B22,
98 FLA_TL );
99
100 FLA_Cont_with_3x1_to_2x1( &yT, y01,
101 psi11,
102 /* ** */ /* ***** */
103 &yB, y21, FLA_TOP );
104 }
105
106 return FLA_SUCCESS;
107 }
108
109