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_nl_unb_var1(FLA_Obj A,FLA_Obj Y,FLA_Obj B)13 FLA_Error FLA_Eig_gest_nl_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 y21_l, y21_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( y21, &y21_l, &y21_r, 1, FLA_LEFT );
60
61 // y21 = A22 * b21;
62 FLA_Hemv_external( FLA_LOWER_TRIANGULAR,
63 FLA_ONE, A22, b21, FLA_ZERO, y21_l );
64
65 // a21 = a21 * beta11;
66 FLA_Scal_external( beta11, a21 );
67
68 // a21 = a21 + 1/2 * y21;
69 FLA_Axpy_external( FLA_ONE_HALF, y21_l, a21 );
70
71 // alpha11 = conj(beta11) * alpha11 * beta11;
72 // = beta11 * alpha11 * beta11;
73 FLA_Scal_external( beta11, alpha11 );
74 FLA_Scal_external( beta11, alpha11 );
75
76 // alpha11 = alpha11 + a21' * b21 + b21' * a21;
77 FLA_Dot2cs_external( FLA_CONJUGATE, FLA_ONE, a21, b21, FLA_ONE, alpha11 );
78
79 // a21 = a21 + 1/2 * y21;
80 FLA_Axpy_external( FLA_ONE_HALF, y21_l, a21 );
81
82 // a21 = tril( B22 )' * a21;
83 FLA_Trmv_external( FLA_LOWER_TRIANGULAR, FLA_CONJ_TRANSPOSE, FLA_NONUNIT_DIAG,
84 B22, a21 );
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