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_var5(FLA_Obj A,FLA_Obj Y,FLA_Obj B)13 FLA_Error FLA_Eig_gest_iu_unb_var5( 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 yL, yR, y10t, psi11, y12t;
24
25 //FLA_Obj y12t_t,
26 // y12t_b;
27
28 FLA_Obj psi11, y12t,
29 y21, Y22;
30
31 FLA_Part_2x2( A, &ATL, &ATR,
32 &ABL, &ABR, 0, 0, FLA_TL );
33
34 FLA_Part_2x2( B, &BTL, &BTR,
35 &BBL, &BBR, 0, 0, FLA_TL );
36
37 //FLA_Part_1x2( Y, &yL, &yR, 0, FLA_LEFT );
38
39 FLA_Part_2x2( Y, &psi11, &y12t,
40 &y21, &Y22, 1, 1, FLA_TL );
41
42 while ( FLA_Obj_length( ATL ) < FLA_Obj_length( A ) ){
43
44 FLA_Repart_2x2_to_3x3( ATL, /**/ ATR, &A00, /**/ &a01, &A02,
45 /* ************* */ /* ************************** */
46 &a10t, /**/ &alpha11, &a12t,
47 ABL, /**/ ABR, &A20, /**/ &a21, &A22,
48 1, 1, FLA_BR );
49
50 FLA_Repart_2x2_to_3x3( BTL, /**/ BTR, &B00, /**/ &b01, &B02,
51 /* ************* */ /* ************************* */
52 &b10t, /**/ &beta11, &b12t,
53 BBL, /**/ BBR, &B20, /**/ &b21, &B22,
54 1, 1, FLA_BR );
55
56 //FLA_Repart_1x2_to_1x3( yL, /**/ yR, &y10t, /**/ &psi11, &y12t,
57 // 1, FLA_RIGHT );
58
59 /*------------------------------------------------------------*/
60
61 //FLA_Part_2x1( y12t, &y12t_t,
62 // &y12t_b, 1, FLA_TOP );
63
64 // alpha11 = inv(beta11) * alpha11 * inv(conj(beta11));
65 // = inv(beta11) * alpha11 * inv(beta11);
66 FLA_Inv_scal_external( beta11, alpha11 );
67 FLA_Inv_scal_external( beta11, alpha11 );
68
69 //// y12t = alpha11 * b12t;
70 //FLA_Copy_external( b12t, y12t_t );
71 //FLA_Scal_external( alpha11, y12t_t );
72 // psi11 = - 1/2 * alpha11;
73 FLA_Copy_external( alpha11, psi11 );
74 FLA_Scal_external( FLA_MINUS_ONE_HALF, psi11 );
75
76 // a12t = inv(conj(beta11)) * a12t;
77 // = inv(beta11) * a12t;
78 FLA_Inv_scal_external( beta11, a12t );
79
80 //// a12t = a12t - 1/2 * y12t;
81 //FLA_Axpy_external( FLA_MINUS_ONE_HALF, y12t_t, a12t );
82 // a12t = a12t - 1/2 * alpha11 * b12t;
83 FLA_Axpy_external( psi11, b12t, a12t );
84
85 // A22 = A22 - a12t' * b12t - b12t' * a12t;
86 FLA_Her2c_external( FLA_UPPER_TRIANGULAR, FLA_CONJUGATE,
87 FLA_MINUS_ONE, a12t, b12t, A22 );
88
89 //// a12t = a12t - 1/2 * y12t;
90 //FLA_Axpy_external( FLA_MINUS_ONE_HALF, y12t_t, a12t );
91 // a12t = a12t - 1/2 * alpha11 * b12t;
92 FLA_Axpy_external( psi11, b12t, a12t );
93
94 // a12t = a12t * inv( tril( B22 ) );
95 // a12t^T = inv( triu( B22 )^T ) * a12t^T;
96 FLA_Trsv_external( FLA_UPPER_TRIANGULAR, FLA_TRANSPOSE, FLA_NONUNIT_DIAG,
97 B22, a12t );
98
99 /*------------------------------------------------------------*/
100
101 FLA_Cont_with_3x3_to_2x2( &ATL, /**/ &ATR, A00, a01, /**/ A02,
102 a10t, alpha11, /**/ a12t,
103 /* ************** */ /* ************************ */
104 &ABL, /**/ &ABR, A20, a21, /**/ A22,
105 FLA_TL );
106
107 FLA_Cont_with_3x3_to_2x2( &BTL, /**/ &BTR, B00, b01, /**/ B02,
108 b10t, beta11, /**/ b12t,
109 /* ************** */ /* *********************** */
110 &BBL, /**/ &BBR, B20, b21, /**/ B22,
111 FLA_TL );
112
113 //FLA_Cont_with_1x3_to_1x2( &yL, /**/ &yR, y10t, psi11, /**/ y12t,
114 // FLA_LEFT );
115 }
116
117 return FLA_SUCCESS;
118 }
119
120