libflame  revision_anchor
Functions
FLA_Tevd_francis_v_opt_var1.c File Reference

(r)

Functions

FLA_Error FLA_Tevd_francis_v_opt_var1 (FLA_Obj shift, FLA_Obj g, FLA_Obj d, FLA_Obj e)
 
FLA_Error FLA_Tevd_francis_v_ops_var1 (int m_A, float *buff_shift, scomplex *buff_g, int inc_g, float *buff_d, int inc_d, float *buff_e, int inc_e)
 
FLA_Error FLA_Tevd_francis_v_opd_var1 (int m_A, double *buff_shift, dcomplex *buff_g, int inc_g, double *buff_d, int inc_d, double *buff_e, int inc_e)
 

Function Documentation

◆ FLA_Tevd_francis_v_opd_var1()

FLA_Error FLA_Tevd_francis_v_opd_var1 ( int  m_A,
double *  buff_shift,
dcomplex buff_g,
int  inc_g,
double *  buff_d,
int  inc_d,
double *  buff_e,
int  inc_e 
)
86 {
87  double eps, safmin;
88  double temp0, temp1;
89  double bulge;
90  int ij_deflated;
91  int i;
92 
93  // Initialize the deflation index.
94  ij_deflated = FLA_SUCCESS;
95 
96  // Initialize the bulge variable to zero.
97  bulge = 0.0;
98 
99  // Query epsilon and safmin.
100  eps = FLA_Mach_params_opd( FLA_MACH_EPS );
101  safmin = FLA_Mach_params_opd( FLA_MACH_SFMIN );
102 
103  // Apply the rotations in forward order.
104  for ( i = 0; i < m_A - 1; ++i )
105  {
106  double* alpha00 = buff_d + (i-1)*inc_d;
107  double* alpha10 = buff_e + (i-1)*inc_e;
108  double* alpha20 = &bulge;
109 
110  double* alpha11 = buff_d + (i )*inc_d;
111  double* alpha21 = buff_e + (i )*inc_e;
112  double* alpha22 = buff_d + (i+1)*inc_d;
113 
114  double* alpha31 = &bulge;
115  double* alpha32 = buff_e + (i+1)*inc_e;
116 
117  double* gamma1 = &(buff_g[(i )*inc_g].real);
118  double* sigma1 = &(buff_g[(i )*inc_g].imag);
119 
120  double alpha10_new;
121 
122  int m_behind = i;
123  int m_ahead = m_A - i - 2;
124 
125  /*------------------------------------------------------------*/
126 
127  if ( i == 0 )
128  {
129  // Induce an iteration that introduces the bulge by
130  // changing the addresses of alpha10 and alpha20.
131  temp0 = *buff_d - *buff_shift;
132  temp1 = *buff_e;
133  alpha10 = &temp0;
134  alpha20 = &temp1;
135 
136  // Compute a new Givens rotation that introduces the bulge.
137  MAC_Givens2_opd( alpha10,
138  //FLA_Givens2_opd( alpha10,
139  alpha20,
140  gamma1,
141  sigma1,
142  &alpha10_new );
143 
144  // We don't apply the Givens rotation to the 2x1 vector at
145  // alpha10 when introducing the bulge.
146  }
147  else
148  {
149  // Compute a new Givens rotation to push the bulge.
150  MAC_Givens2_opd( alpha10,
151  //FLA_Givens2_opd( alpha10,
152  alpha20,
153  gamma1,
154  sigma1,
155  &alpha10_new );
156 
157  // Apply the Givens rotation to the 2x1 vector from which it
158  // was computed, which annihilates alpha20.
159  *alpha10 = alpha10_new;
160  *alpha20 = 0.0;
161  }
162 
163  // Apply the Givens rotation to the 2x2 submatrix at alpha11.
164  MAC_Apply_GTG_opd( gamma1,
165  sigma1,
166  alpha11,
167  alpha21,
168  alpha22 );
169 
170  if ( m_ahead > 0 )
171  {
172  // Apply the Givens rotation to the 1x2 vector below the 2x2
173  // submatrix. This should move the bulge to alpha31.
174  MAC_Apply_G_1x2_opd( gamma1,
175  sigma1,
176  alpha31,
177  alpha32 );
178 
179  // Check for deflation after applying the rotations, except after
180  // applying the initial bulge-introducing rotations.
181  if ( m_behind > 0 )
182  {
183  // We check for deflation in the previous column now that we
184  // are done modifying it. If deflation occurred, record the
185  // index.
186  if ( MAC_Tevd_eigval_converged_opd( eps, safmin, *alpha00, *alpha10, *alpha11 ) )
187  {
188  ij_deflated = i - 1;
189  }
190  }
191 
192  // Sanity check. If the bulge ever disappears, something may be wrong.
193  if ( *alpha31 == 0.0 )
194  {
195  printf( "FLA_Tevd_francis_v_opt_var1: bulge disappeared!\n" );
196  if ( MAC_Tevd_eigval_converged_opd( eps, safmin, *alpha11, *alpha21, *alpha22 ) )
197  {
198  printf( "FLA_Tevd_francis_v_opt_var1: deflation detected (col %d)\n", i );
199  printf( "FLA_Tevd_francis_v_opt_var1: alpha11 = %23.19e\n", *alpha11 );
200  printf( "FLA_Tevd_francis_v_opt_var1: alpha21 alpha22 = %23.19e %23.19e\n", *alpha21, *alpha22 );
201  return i;
202  }
203  else
204  {
205  printf( "FLA_Tevd_francis_v_opt_var1: but NO deflation detected! (col %d)\n", i );
206  printf( "FLA_Tevd_francis_v_opt_var1: alpha11 = %23.19e\n", *alpha11 );
207  printf( "FLA_Tevd_francis_v_opt_var1: alpha21 alpha22 = %23.19e %23.19e\n", *alpha21, *alpha22 );
208  FLA_Abort();
209  return FLA_FAILURE;
210  }
211  }
212  }
213 
214  /*------------------------------------------------------------*/
215  }
216 
217  // Return the index of column where deflation most recently occurred,
218  // or FLA_SUCCESS if no deflation was detected.
219  return ij_deflated;
220 }
float real
Definition: FLA_f2c.h:30
void FLA_Abort(void)
Definition: FLA_Error.c:248
double FLA_Mach_params_opd(FLA_Machval machval)
Definition: FLA_Mach_params.c:74
int i
Definition: bl1_axmyv2.c:145
double temp1
Definition: bl1_axpyv2b.c:146
rho_c imag
Definition: bl1_axpyv2bdotaxpy.c:483

References FLA_Abort(), FLA_Mach_params_opd(), i, imag, and temp1.

Referenced by FLA_Tevd_eigval_v_opd_var1(), FLA_Tevd_eigval_v_opd_var3(), and FLA_Tevd_francis_v_opt_var1().

◆ FLA_Tevd_francis_v_ops_var1()

FLA_Error FLA_Tevd_francis_v_ops_var1 ( int  m_A,
float *  buff_shift,
scomplex buff_g,
int  inc_g,
float *  buff_d,
int  inc_d,
float *  buff_e,
int  inc_e 
)
75 {
76  return FLA_SUCCESS;
77 }

Referenced by FLA_Tevd_francis_v_opt_var1().

◆ FLA_Tevd_francis_v_opt_var1()

FLA_Error FLA_Tevd_francis_v_opt_var1 ( FLA_Obj  shift,
FLA_Obj  g,
FLA_Obj  d,
FLA_Obj  e 
)
14 {
15  FLA_Datatype datatype;
16  int m_A;
17  int inc_g;
18  int inc_d;
19  int inc_e;
20 
21  datatype = FLA_Obj_datatype( d );
22 
23  m_A = FLA_Obj_vector_dim( d );
24 
25  inc_g = FLA_Obj_vector_inc( g );
26  inc_d = FLA_Obj_vector_inc( d );
27  inc_e = FLA_Obj_vector_inc( e );
28 
29 
30  switch ( datatype )
31  {
32  case FLA_FLOAT:
33  {
34  float* buff_shift = FLA_FLOAT_PTR( shift );
35  scomplex* buff_g = FLA_COMPLEX_PTR( g );
36  float* buff_d = FLA_FLOAT_PTR( d );
37  float* buff_e = FLA_FLOAT_PTR( e );
38 
40  buff_shift,
41  buff_g, inc_g,
42  buff_d, inc_d,
43  buff_e, inc_e );
44 
45  break;
46  }
47 
48  case FLA_DOUBLE:
49  {
50  double* buff_shift = FLA_DOUBLE_PTR( shift );
51  dcomplex* buff_g = FLA_DOUBLE_COMPLEX_PTR( g );
52  double* buff_d = FLA_DOUBLE_PTR( d );
53  double* buff_e = FLA_DOUBLE_PTR( e );
54 
56  buff_shift,
57  buff_g, inc_g,
58  buff_d, inc_d,
59  buff_e, inc_e );
60 
61  break;
62  }
63  }
64 
65  return FLA_SUCCESS;
66 }
FLA_Error FLA_Tevd_francis_v_opd_var1(int m_A, double *buff_shift, dcomplex *buff_g, int inc_g, double *buff_d, int inc_d, double *buff_e, int inc_e)
Definition: FLA_Tevd_francis_v_opt_var1.c:81
FLA_Error FLA_Tevd_francis_v_ops_var1(int m_A, float *buff_shift, scomplex *buff_g, int inc_g, float *buff_d, int inc_d, float *buff_e, int inc_e)
Definition: FLA_Tevd_francis_v_opt_var1.c:70
dim_t FLA_Obj_vector_inc(FLA_Obj obj)
Definition: FLA_Query.c:145
dim_t FLA_Obj_vector_dim(FLA_Obj obj)
Definition: FLA_Query.c:137
FLA_Datatype FLA_Obj_datatype(FLA_Obj obj)
Definition: FLA_Query.c:13
int FLA_Datatype
Definition: FLA_type_defs.h:49
Definition: blis_type_defs.h:138
Definition: blis_type_defs.h:133

References FLA_Obj_datatype(), FLA_Obj_vector_dim(), FLA_Obj_vector_inc(), FLA_Tevd_francis_v_opd_var1(), and FLA_Tevd_francis_v_ops_var1().