libflame  revision_anchor
Functions
FLA_Tevd_find_submatrix.c File Reference

(r)

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)
 

Function Documentation

◆ 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 
)

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().

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

◆ 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 }