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