libflame  revision_anchor
Functions
FLA_Tevd_n.h File Reference

(r)

Go to the source code of this file.

Functions

FLA_Error FLA_Tevd_find_submatrix_ops (int m_A, int ij_begin, float *buff_d, int inc_d, float *buff_e, int inc_e, int *ijTL, int *ijBR)
 
FLA_Error FLA_Tevd_find_submatrix_opd (int m_A, int ij_begin, double *buff_d, int inc_d, double *buff_e, int inc_e, int *ijTL, int *ijBR)
 
FLA_Error FLA_Norm1_tridiag (FLA_Obj d, FLA_Obj e, FLA_Obj norm)
 
FLA_Error FLA_Norm1_tridiag_ops (int m_A, float *buff_d, int inc_d, float *buff_e, int inc_e, float *norm)
 
FLA_Error FLA_Norm1_tridiag_opd (int m_A, double *buff_d, int inc_d, double *buff_e, int inc_e, double *norm)
 
FLA_Error FLA_Tevd_n_opt_var1 (dim_t n_iter_max, FLA_Obj d, FLA_Obj e, FLA_Obj G, FLA_Obj U)
 
FLA_Error FLA_Tevd_n_ops_var1 (int m_A, int m_U, int n_G, int n_iter_max, float *buff_d, int inc_d, float *buff_e, int inc_e, scomplex *buff_G, int rs_G, int cs_G)
 
FLA_Error FLA_Tevd_n_opd_var1 (int m_A, int m_U, int n_G, int n_iter_max, double *buff_d, int inc_d, double *buff_e, int inc_e, dcomplex *buff_G, int rs_G, int cs_G)
 
FLA_Error FLA_Tevd_n_opc_var1 (int m_A, int m_U, int n_G, int n_iter_max, float *buff_d, int inc_d, float *buff_e, int inc_e, scomplex *buff_G, int rs_G, int cs_G)
 
FLA_Error FLA_Tevd_n_opz_var1 (int m_A, int m_U, int n_G, int n_iter_max, double *buff_d, int inc_d, double *buff_e, int inc_e, dcomplex *buff_G, int rs_G, int cs_G)
 

Function Documentation

◆ FLA_Norm1_tridiag()

FLA_Error FLA_Norm1_tridiag ( FLA_Obj  d,
FLA_Obj  e,
FLA_Obj  norm 
)
14 {
15  FLA_Datatype datatype;
16  int m_A;
17  int inc_d;
18  int inc_e;
19 
20  datatype = FLA_Obj_datatype( d );
21 
22  m_A = FLA_Obj_vector_dim( d );
23 
24  inc_d = FLA_Obj_vector_inc( d );
25  inc_e = FLA_Obj_vector_inc( e );
26 
27 
28  switch ( datatype )
29  {
30  case FLA_FLOAT:
31  {
32  float* buff_d = FLA_FLOAT_PTR( d );
33  float* buff_e = FLA_FLOAT_PTR( e );
34  float* buff_norm = FLA_FLOAT_PTR( norm );
35 
37  buff_d, inc_d,
38  buff_e, inc_e,
39  buff_norm );
40 
41  break;
42  }
43 
44  case FLA_DOUBLE:
45  {
46  double* buff_d = FLA_DOUBLE_PTR( d );
47  double* buff_e = FLA_DOUBLE_PTR( e );
48  double* buff_norm = FLA_DOUBLE_PTR( norm );
49 
51  buff_d, inc_d,
52  buff_e, inc_e,
53  buff_norm );
54 
55  break;
56  }
57  }
58 
59  return FLA_SUCCESS;
60 }
FLA_Error FLA_Norm1_tridiag_opd(int m_A, double *buff_d, int inc_d, double *buff_e, int inc_e, double *norm)
Definition: FLA_Norm1_tridiag.c:111
FLA_Error FLA_Norm1_tridiag_ops(int m_A, float *buff_d, int inc_d, float *buff_e, int inc_e, float *norm)
Definition: FLA_Norm1_tridiag.c:64
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

◆ FLA_Norm1_tridiag_opd()

FLA_Error FLA_Norm1_tridiag_opd ( int  m_A,
double *  buff_d,
int  inc_d,
double *  buff_e,
int  inc_e,
double *  norm 
)
115 {
116  double* d = buff_d;
117  double* e = buff_e;
118  double nm;
119  int i;
120 
121  if ( m_A == 1 )
122  {
123  nm = fabs( *d );
124  }
125  else
126  {
127  double d_first = d[ (0 )*inc_d ];
128  double e_first = e[ (0 )*inc_e ];
129  double e_last = e[ (m_A-2)*inc_e ];
130  double d_last = d[ (m_A-1)*inc_d ];
131 
132  // Record the maximum of the absolute row/column sums for the
133  // first and last row/columns.
134  nm = max( fabs( d_first ) + fabs( e_first ),
135  fabs( e_last ) + fabs( d_last ) );
136 
137  for ( i = 1; i < m_A - 2; ++i )
138  {
139  double e0 = e[ (i-1)*inc_e ];
140  double e1 = e[ (i )*inc_e ];
141  double d1 = d[ (i )*inc_d ];
142 
143  // Update nm with the absolute row/column sum for the ith
144  // row/column.
145  nm = max( nm, fabs( e0 ) +
146  fabs( d1 ) +
147  fabs( e1 ) );
148  }
149  }
150 
151  *norm = nm;
152 
153  return FLA_SUCCESS;
154 }
int i
Definition: bl1_axmyv2.c:145

◆ FLA_Norm1_tridiag_ops()

FLA_Error FLA_Norm1_tridiag_ops ( int  m_A,
float *  buff_d,
int  inc_d,
float *  buff_e,
int  inc_e,
float *  norm 
)
68 {
69  float* d = buff_d;
70  float* e = buff_e;
71  float nm;
72  int i;
73 
74  if ( m_A == 1 )
75  {
76  nm = fabs( *d );
77  }
78  else
79  {
80  float d_first = d[ (0 )*inc_d ];
81  float e_first = e[ (0 )*inc_e ];
82  float e_last = e[ (m_A-2)*inc_e ];
83  float d_last = d[ (m_A-1)*inc_d ];
84 
85  // Record the maximum of the absolute row/column sums for the
86  // first and last row/columns.
87  nm = max( fabs( d_first ) + fabs( e_first ),
88  fabs( e_last ) + fabs( d_last ) );
89 
90  for ( i = 1; i < m_A - 2; ++i )
91  {
92  float e0 = e[ (i-1)*inc_e ];
93  float e1 = e[ (i )*inc_e ];
94  float d1 = d[ (i )*inc_d ];
95 
96  // Update nm with the absolute row/column sum for the ith
97  // row/column.
98  nm = max( nm, fabs( e0 ) +
99  fabs( d1 ) +
100  fabs( e1 ) );
101  }
102  }
103 
104  *norm = nm;
105 
106  return FLA_SUCCESS;
107 }

◆ FLA_Tevd_find_submatrix_opd()

FLA_Error FLA_Tevd_find_submatrix_opd ( int  m_A,
int  ij_begin,
double *  buff_d,
int  inc_d,
double *  buff_e,
int  inc_e,
int *  ijTL,
int *  ijBR 
)
34 {
35  double rzero = bl1_d0();
36  double eps;
37  int ij_tl;
38  int ij_br;
39 
40  // Initialize some numerical constants.
41  eps = FLA_Mach_params_opd( FLA_MACH_EPS );
42 
43  // Search for the first non-zero subdiagonal element starting at
44  // the index specified by ij_begin.
45  for ( ij_tl = ij_begin; ij_tl < m_A - 1; ++ij_tl )
46  {
47  double* d1 = buff_d + (ij_tl )*inc_d;
48  double* d2 = buff_d + (ij_tl+1)*inc_d;
49  double* e1 = buff_e + (ij_tl )*inc_e;
50  double abs_e1 = fabs( *e1 );
51 
52  // If we encounter a non-zero subdiagonal element that is close
53  // enough to zero, set it to zero.
54  if ( abs_e1 != rzero )
55  {
56  if ( abs_e1 <= eps * sqrt( fabs( *d1 ) ) *
57  sqrt( fabs( *d2 ) ) )
58  {
59 #ifdef PRINTF
60 printf( "FLA_Tevd_find_submatrix_opd: nudging non-zero subdiagonal element (e1) to zero.\n" );
61 printf( " d[%3d] = %22.19e\n", ij_tl, *d1 );
62 printf( " e[%3d] d[%3d] = %22.19e %22.19e\n", ij_tl, ij_tl+1, *e1, *d2 );
63 #endif
64  *e1 = rzero;
65  }
66  }
67 
68  // If we find a non-zero element, record it and break out of this
69  // loop.
70  if ( *e1 != rzero )
71  {
72 #ifdef PRINTF
73 printf( "FLA_Tevd_find_submatrix_opd: found non-zero subdiagonal element\n" );
74 printf( " e[%3d] = %22.19e\n", ij_tl, *e1 );
75 #endif
76  *ijTL = ij_tl;
77  break;
78  }
79  }
80 
81  // If ij_tl was incremented all the way up to m_A - 1, then we didn't
82  // find any non-zeros.
83  if ( ij_tl == m_A - 1 )
84  {
85 #ifdef PRINTF
86 printf( "FLA_Tevd_find_submatrix_opd: no submatrices found.\n" );
87 #endif
88  return FLA_FAILURE;
89  }
90 
91  // If we've gotten this far, then a non-zero subdiagonal element was
92  // found. Now we must walk the remaining portion of the subdiagonal
93  // to find the first zero element, or if one is not found, we simply
94  // use the last element of the subdiagonal.
95  for ( ij_br = ij_tl; ij_br < m_A - 1; ++ij_br )
96  {
97  double* d1 = buff_d + (ij_br )*inc_d;
98  double* d2 = buff_d + (ij_br+1)*inc_d;
99  double* e1 = buff_e + (ij_br )*inc_e;
100  double abs_e1 = fabs( *e1 );
101 
102  // If we encounter a non-zero subdiagonal element that is close
103  // enough to zero, set it to zero.
104  if ( abs_e1 != rzero )
105  {
106  if ( abs_e1 <= eps * sqrt( fabs( *d1 ) ) *
107  sqrt( fabs( *d2 ) ) )
108  {
109 #ifdef PRINTF
110 printf( "FLA_Tevd_find_submatrix_opd: nudging non-zero subdiagonal element (e1) to zero.\n" );
111 printf( " d[%3d] = %22.19e\n", ij_br, *d1 );
112 printf( " e[%3d] d[%3d] = %22.19e %22.19e\n", ij_br, ij_br+1, *e1, *d2 );
113 #endif
114  *e1 = rzero;
115  }
116  }
117 
118  // If we find a zero element, record it and break out of this
119  // loop.
120  if ( *e1 == rzero )
121  {
122 #ifdef PRINTF
123 printf( "FLA_Tevd_find_submatrix_opd: found zero subdiagonal element\n" );
124 printf( " e[%3d] = %22.19e\n", ij_br, *e1 );
125 #endif
126  break;
127  }
128  }
129 
130  // If a zero element was found, then ij_br should hold the index of
131  // that element. If a zero element was not found, then ij_br should
132  // hold m_A - 1. Either way, we save the value and return success.
133  *ijBR = ij_br;
134 
135  return FLA_SUCCESS;
136 }
double FLA_Mach_params_opd(FLA_Machval machval)
Definition: FLA_Mach_params.c:74
double bl1_d0(void)
Definition: bl1_constants.c:118

References bl1_d0(), and FLA_Mach_params_opd().

Referenced by FLA_Tevd_n_opz_var1(), FLA_Tevd_v_opd_var1(), FLA_Tevd_v_opd_var2(), FLA_Tevd_v_opz_var1(), and FLA_Tevd_v_opz_var2().

◆ FLA_Tevd_find_submatrix_ops()

FLA_Error FLA_Tevd_find_submatrix_ops ( int  m_A,
int  ij_begin,
float *  buff_d,
int  inc_d,
float *  buff_e,
int  inc_e,
int *  ijTL,
int *  ijBR 
)
20 {
21  FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );
22 
23  return FLA_SUCCESS;
24 }

◆ FLA_Tevd_n_opc_var1()

FLA_Error FLA_Tevd_n_opc_var1 ( int  m_A,
int  m_U,
int  n_G,
int  n_iter_max,
float *  buff_d,
int  inc_d,
float *  buff_e,
int  inc_e,
scomplex buff_G,
int  rs_G,
int  cs_G 
)
162 {
163  return FLA_SUCCESS;
164 }

Referenced by FLA_Tevd_n_opt_var1().

◆ FLA_Tevd_n_opd_var1()

FLA_Error FLA_Tevd_n_opd_var1 ( int  m_A,
int  m_U,
int  n_G,
int  n_iter_max,
double *  buff_d,
int  inc_d,
double *  buff_e,
int  inc_e,
dcomplex buff_G,
int  rs_G,
int  cs_G 
)
151 {
152  return FLA_SUCCESS;
153 }

Referenced by FLA_Tevd_n_opt_var1().

◆ FLA_Tevd_n_ops_var1()

FLA_Error FLA_Tevd_n_ops_var1 ( int  m_A,
int  m_U,
int  n_G,
int  n_iter_max,
float *  buff_d,
int  inc_d,
float *  buff_e,
int  inc_e,
scomplex buff_G,
int  rs_G,
int  cs_G 
)
138 {
139  return FLA_SUCCESS;
140 }

Referenced by FLA_Tevd_n_opt_var1().

◆ FLA_Tevd_n_opt_var1()

FLA_Error FLA_Tevd_n_opt_var1 ( dim_t  n_iter_max,
FLA_Obj  d,
FLA_Obj  e,
FLA_Obj  G,
FLA_Obj  U 
)
14 {
15  FLA_Error r_val = FLA_SUCCESS;
16  FLA_Datatype datatype;
17  int m_A, m_U, n_G;
18  int inc_d;
19  int inc_e;
20  int rs_G, cs_G;
21 
22  datatype = FLA_Obj_datatype( U );
23 
24  m_A = FLA_Obj_vector_dim( d );
25  m_U = FLA_Obj_vector_dim( d );
26  n_G = FLA_Obj_width( G );
27 
28  inc_d = FLA_Obj_vector_inc( d );
29  inc_e = FLA_Obj_vector_inc( e );
30 
31  rs_G = FLA_Obj_row_stride( G );
32  cs_G = FLA_Obj_col_stride( G );
33 
34 /*
35 FLA_Obj de, deL, deR, deLT, deLB;
36 FLA_Obj_create( FLA_DOUBLE, m_A, 2, 0, 0, &de );
37 FLA_Part_1x2( de, &deL, &deR, 1, FLA_LEFT );
38 FLA_Part_2x1( deL, &deLT,
39  &deLB, 1, FLA_BOTTOM );
40 FLA_Copy( d, deR );
41 FLA_Copy( e, deLT );
42 FLA_Set( FLA_ZERO, deLB );
43 //FLA_Obj_show( "de", de, "%21.17e", "" );
44 */
45 
46  switch ( datatype )
47  {
48  case FLA_FLOAT:
49  {
50  float* buff_d = FLA_FLOAT_PTR( d );
51  float* buff_e = FLA_FLOAT_PTR( e );
52  scomplex* buff_G = FLA_COMPLEX_PTR( G );
53 
54  r_val = FLA_Tevd_n_ops_var1( m_A,
55  m_U,
56  n_G,
57  n_iter_max,
58  buff_d, inc_d,
59  buff_e, inc_e,
60  buff_G, rs_G, cs_G );
61 
62  break;
63  }
64 
65  case FLA_DOUBLE:
66  {
67  double* buff_d = FLA_DOUBLE_PTR( d );
68  double* buff_e = FLA_DOUBLE_PTR( e );
69  dcomplex* buff_G = FLA_DOUBLE_COMPLEX_PTR( G );
70 
71  r_val = FLA_Tevd_n_opd_var1( m_A,
72  m_U,
73  n_G,
74  n_iter_max,
75  buff_d, inc_d,
76  buff_e, inc_e,
77  buff_G, rs_G, cs_G );
78 
79  break;
80  }
81 
82  case FLA_COMPLEX:
83  {
84  float* buff_d = FLA_FLOAT_PTR( d );
85  float* buff_e = FLA_FLOAT_PTR( e );
86  scomplex* buff_G = FLA_COMPLEX_PTR( G );
87 
88  r_val = FLA_Tevd_n_opc_var1( m_A,
89  m_U,
90  n_G,
91  n_iter_max,
92  buff_d, inc_d,
93  buff_e, inc_e,
94  buff_G, rs_G, cs_G );
95 
96  break;
97  }
98 
99  case FLA_DOUBLE_COMPLEX:
100  {
101  double* buff_d = FLA_DOUBLE_PTR( d );
102  double* buff_e = FLA_DOUBLE_PTR( e );
103  dcomplex* buff_G = FLA_DOUBLE_COMPLEX_PTR( G );
104 
105  r_val = FLA_Tevd_n_opz_var1( m_A,
106  m_U,
107  n_G,
108  n_iter_max,
109  buff_d, inc_d,
110  buff_e, inc_e,
111  buff_G, rs_G, cs_G );
112 
113  break;
114  }
115  }
116 /*
117 FLA_Copy( d, deR );
118 FLA_Copy( e, deLT );
119 FLA_Set( FLA_ZERO, deLB );
120 FLA_Sort( FLA_FORWARD, deR );
121 FLA_Obj_show( "de after", de, "%21.17e", "" );
122 double eps = FLA_Mach_params_opd( FLA_MACH_EPS );
123 printf( "epsilon = %21.17e\n", eps );
124 FLA_Obj_free( &de );
125 */
126  return r_val;
127 }
FLA_Error FLA_Tevd_n_opd_var1(int m_A, int m_U, int n_G, int n_iter_max, double *buff_d, int inc_d, double *buff_e, int inc_e, dcomplex *buff_G, int rs_G, int cs_G)
Definition: FLA_Tevd_n_opt_var1.c:144
FLA_Error FLA_Tevd_n_opc_var1(int m_A, int m_U, int n_G, int n_iter_max, float *buff_d, int inc_d, float *buff_e, int inc_e, scomplex *buff_G, int rs_G, int cs_G)
Definition: FLA_Tevd_n_opt_var1.c:155
FLA_Error FLA_Tevd_n_opz_var1(int m_A, int m_U, int n_G, int n_iter_max, double *buff_d, int inc_d, double *buff_e, int inc_e, dcomplex *buff_G, int rs_G, int cs_G)
Definition: FLA_Tevd_n_opt_var1.c:168
FLA_Error FLA_Tevd_n_ops_var1(int m_A, int m_U, int n_G, int n_iter_max, float *buff_d, int inc_d, float *buff_e, int inc_e, scomplex *buff_G, int rs_G, int cs_G)
Definition: FLA_Tevd_n_opt_var1.c:131
dim_t FLA_Obj_width(FLA_Obj obj)
Definition: FLA_Query.c:123
dim_t FLA_Obj_row_stride(FLA_Obj obj)
Definition: FLA_Query.c:167
dim_t FLA_Obj_col_stride(FLA_Obj obj)
Definition: FLA_Query.c:174
int FLA_Error
Definition: FLA_type_defs.h:47
Definition: blis_type_defs.h:138
Definition: blis_type_defs.h:133

References FLA_Obj_col_stride(), FLA_Obj_datatype(), FLA_Obj_row_stride(), FLA_Obj_vector_dim(), FLA_Obj_vector_inc(), FLA_Obj_width(), FLA_Tevd_n_opc_var1(), FLA_Tevd_n_opd_var1(), FLA_Tevd_n_ops_var1(), and FLA_Tevd_n_opz_var1().

◆ FLA_Tevd_n_opz_var1()

FLA_Error FLA_Tevd_n_opz_var1 ( int  m_A,
int  m_U,
int  n_G,
int  n_iter_max,
double *  buff_d,
int  inc_d,
double *  buff_e,
int  inc_e,
dcomplex buff_G,
int  rs_G,
int  cs_G 
)
175 {
176  dcomplex one = bl1_z1();
177  double rone = bl1_d1();
178 
179  double eps;
180  double eps2;
181  double safmin;
182  double ssfmin;
183  double safmax;
184  double ssfmax;
185 
186  dcomplex* G;
187  double* d1;
188  double* e1;
189  int r_val;
190  int done;
191  int m_G_sweep_max;
192  int ij_begin;
193  int ijTL, ijBR;
194  int m_A11;
195  int n_iter_perf;
196  int total_deflations;
197  int n_deflations;
198  int n_iter_prev;
199  int n_iter_perf_sweep_max;
200 
201  // Initialize some numerical constants.
202  eps = FLA_Mach_params_opd( FLA_MACH_EPS );
203  eps2 = FLA_Mach_params_opd( FLA_MACH_EPS2 );
204  safmin = FLA_Mach_params_opd( FLA_MACH_SFMIN );
205  safmax = rone / safmin;
206  ssfmax = sqrt( safmax ) / 3.0;
207  ssfmin = sqrt( safmin ) / eps2;
208 
209  // Initialize our completion flag.
210  done = FALSE;
211 
212  // Initialize a counter that holds the maximum number of rows of G
213  // that we would need to initialize for the next sweep.
214  m_G_sweep_max = m_A - 1;
215 
216  // Initialize a counter for the total number of iterations performed.
217  n_iter_prev = 0;
218 
219  // Iterate until the matrix has completely deflated.
220  for ( total_deflations = 0; done != TRUE; )
221  {
222 
223  // Initialize G to contain only identity rotations.
224  bl1_zsetm( m_G_sweep_max,
225  n_G,
226  &one,
227  buff_G, rs_G, cs_G );
228 
229  // Keep track of the maximum number of iterations performed in the
230  // current sweep. This is used when applying the sweep's Givens
231  // rotations.
232  n_iter_perf_sweep_max = 0;
233 
234  // Perform a sweep: Move through the matrix and perform a tridiagonal
235  // EVD on each non-zero submatrix that is encountered. During the
236  // first time through, ijTL will be 0 and ijBR will be m_A - 1.
237  for ( ij_begin = 0; ij_begin < m_A; )
238  {
239 
240 #ifdef PRINTF
241 if ( ij_begin == 0 )
242 printf( "FLA_Tevd_n_opz_var1: beginning new sweep (ij_begin = %d)\n", ij_begin );
243 #endif
244 
245  // Search for the first submatrix along the diagonal that is
246  // bounded by zeroes (or endpoints of the matrix). If no
247  // submatrix is found (ie: if the entire subdiagonal is zero
248  // then FLA_FAILURE is returned. This function also inspects
249  // subdiagonal elements for proximity to zero. If a given
250  // element is close enough to zero, then it is deemed
251  // converged and manually set to zero.
252  r_val = FLA_Tevd_find_submatrix_opd( m_A,
253  ij_begin,
254  buff_d, inc_d,
255  buff_e, inc_e,
256  &ijTL,
257  &ijBR );
258 
259  // Verify that a submatrix was found. If one was not found,
260  // then we are done with the current sweep. Furthermore, if
261  // a submatrix was not found AND we began our search at the
262  // beginning of the matrix (ie: ij_begin == 0), then the
263  // matrix has completely deflated and so we are done with
264  // Francis step iteration.
265  if ( r_val == FLA_FAILURE )
266  {
267  if ( ij_begin == 0 )
268  {
269 #ifdef PRINTF
270 printf( "FLA_Tevd_n_opz_var1: subdiagonal is completely zero.\n" );
271 printf( "FLA_Tevd_n_opz_var1: Francis iteration is done!\n" );
272 #endif
273  done = TRUE;
274  }
275 
276  // Break out of the current sweep so we can apply the last
277  // remaining Givens rotations.
278  break;
279  }
280 
281  // If we got this far, then:
282  // (a) ijTL refers to the index of the first non-zero
283  // subdiagonal along the diagonal, and
284  // (b) ijBR refers to either:
285  // - the first zero element that occurs after ijTL, or
286  // - the the last diagonal element.
287  // Note that ijTL and ijBR also correspond to the first and
288  // last diagonal elements of the submatrix of interest. Thus,
289  // we may compute the dimension of this submatrix as:
290  m_A11 = ijBR - ijTL + 1;
291 
292 #ifdef PRINTF
293 printf( "FLA_Tevd_n_opz_var1: ij_begin = %d\n", ij_begin );
294 printf( "FLA_Tevd_n_opz_var1: ijTL = %d\n", ijTL );
295 printf( "FLA_Tevd_n_opz_var1: ijBR = %d\n", ijBR );
296 printf( "FLA_Tevd_n_opz_var1: m_A11 = %d\n", m_A11 );
297 #endif
298 
299  // Adjust ij_begin, which gets us ready for the next submatrix
300  // search in the current sweep.
301  ij_begin = ijBR + 1;
302 
303  // Index to the submatrices upon which we will operate.
304  d1 = buff_d + ijTL * inc_d;
305  e1 = buff_e + ijTL * inc_e;
306  G = buff_G + ijTL * rs_G;
307 
308  // Compute the 1-norm (which equals the infinity norm since the
309  // matrix is tridiagonal and symmetric) and perform appropriate
310  // scaling if necessary.
311 /*
312  FLA_Norm1_tridiag( m_A11,
313  d1, inc_d,
314  e1, inc_e,
315  &norm );
316 */
317 
318  // Search for a batch of eigenvalues, recursing on deflated
319  // subproblems whenever a split occurs. Iteration continues
320  // as long as:
321  // (a) there is still matrix left to operate on, and
322  // (b) the number of iterations performed in this batch is
323  // less than n_G.
324  // If/when either of the two above conditions fails to hold,
325  // the function returns.
326  n_deflations = FLA_Tevd_iteracc_n_opd_var1( m_A11,
327  n_G,
328  ijTL,
329  d1, inc_d,
330  e1, inc_e,
331  &n_iter_perf );
332 
333  // Record the number of deflations that were observed.
334  total_deflations += n_deflations;
335 
336  // Update the maximum number of iterations performed in the
337  // current sweep.
338  n_iter_perf_sweep_max = max( n_iter_perf_sweep_max, n_iter_perf );
339 
340 #ifdef PRINTF
341 printf( "FLA_Tevd_n_opz_var1: deflations observed = %d\n", n_deflations );
342 printf( "FLA_Tevd_n_opz_var1: total deflations observed = %d\n", total_deflations );
343 printf( "FLA_Tevd_n_opz_var1: num iterations performed = %d\n", n_iter_perf );
344 #endif
345 
346  // Store the most recent value of ijBR in m_G_sweep_max.
347  // When the sweep is done, this value will contain the minimum
348  // number of rows of G we can apply and safely include all
349  // non-identity rotations that were computed during the
350  // eigenvalue searches.
351  m_G_sweep_max = ijBR;
352 
353  // Make sure we haven't exceeded our maximum iteration count.
354  if ( n_iter_prev >= m_A * n_iter_max )
355  {
356 #ifdef PRINTF
357 printf( "FLA_Tevd_n_opz_var1: reached maximum total number of iterations: %d\n", n_iter_prev );
358 #endif
359  FLA_Abort();
360  //return FLA_FAILURE;
361  }
362  }
363 
364  // The sweep is complete.
365 
366  // Increment the total number of iterations previously performed.
367  n_iter_prev += n_iter_perf_sweep_max;
368 
369 #ifdef PRINTF
370 printf( "FLA_Tevd_n_opz_var1: total number of iterations performed: %d\n", n_iter_prev );
371 #endif
372  }
373 
374  //return FLA_SUCCESS;
375  return n_iter_prev;
376 }
FLA_Error FLA_Tevd_iteracc_n_opd_var1(int m_A, int n_G, int ijTL, double *buff_d, int inc_d, double *buff_e, int inc_e, int *n_iter_perf)
Definition: FLA_Tevd_iteracc_n_opt_var1.c:25
FLA_Error FLA_Tevd_find_submatrix_opd(int m_A, int ij_begin, double *buff_d, int inc_d, double *buff_e, int inc_e, int *ijTL, int *ijBR)
Definition: FLA_Tevd_find_submatrix.c:28
void FLA_Abort(void)
Definition: FLA_Error.c:248
dcomplex bl1_z1(void)
Definition: bl1_constants.c:69
double bl1_d1(void)
Definition: bl1_constants.c:54
void bl1_zsetm(int m, int n, dcomplex *sigma, dcomplex *a, int a_rs, int a_cs)
Definition: bl1_setm.c:78

References bl1_d1(), bl1_z1(), bl1_zsetm(), FLA_Abort(), FLA_Mach_params_opd(), FLA_Tevd_find_submatrix_opd(), and FLA_Tevd_iteracc_n_opd_var1().

Referenced by FLA_Tevd_n_opt_var1().