libflame  revision_anchor
Functions
FLA_Wilkshift_tridiag.c File Reference

(r)

Functions

FLA_Error FLA_Wilkshift_tridiag (FLA_Obj delta1, FLA_Obj epsilon, FLA_Obj delta2, FLA_Obj kappa)
 
FLA_Error FLA_Wilkshift_tridiag_ops (float delta1, float epsilon, float delta2, float *kappa)
 
double fla_dlapy2 (double x, double y)
 
FLA_Error FLA_Wilkshift_tridiag_opd (double delta1, double epsilon, double delta2, double *kappa)
 

Function Documentation

◆ fla_dlapy2()

double fla_dlapy2 ( double  x,
double  y 
)

◆ FLA_Wilkshift_tridiag()

FLA_Error FLA_Wilkshift_tridiag ( FLA_Obj  delta1,
FLA_Obj  epsilon,
FLA_Obj  delta2,
FLA_Obj  kappa 
)
59 {
60  FLA_Datatype datatype;
61 
62  datatype = FLA_Obj_datatype( delta1 );
63 
64  if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
65  FLA_Wilkshift_tridiag_check( delta1, epsilon, delta2, kappa );
66 
67  switch ( datatype )
68  {
69  case FLA_FLOAT:
70  {
71  float* delta1_p = ( float* ) FLA_FLOAT_PTR( delta1 );
72  float* epsilon_p = ( float* ) FLA_FLOAT_PTR( epsilon );
73  float* delta2_p = ( float* ) FLA_FLOAT_PTR( delta2 );
74  float* kappa_p = ( float* ) FLA_FLOAT_PTR( kappa );
75 
76  FLA_Wilkshift_tridiag_ops( *delta1_p,
77  *epsilon_p,
78  *delta2_p,
79  kappa_p );
80 
81  break;
82  }
83 
84  case FLA_DOUBLE:
85  {
86  double* delta1_p = ( double* ) FLA_DOUBLE_PTR( delta1 );
87  double* epsilon_p = ( double* ) FLA_DOUBLE_PTR( epsilon );
88  double* delta2_p = ( double* ) FLA_DOUBLE_PTR( delta2 );
89  double* kappa_p = ( double* ) FLA_DOUBLE_PTR( kappa );
90 
91  FLA_Wilkshift_tridiag_opd( *delta1_p,
92  *epsilon_p,
93  *delta2_p,
94  kappa_p );
95 
96  break;
97  }
98  }
99 
100  return FLA_SUCCESS;
101 }
FLA_Error FLA_Wilkshift_tridiag_ops(float delta1, float epsilon, float delta2, float *kappa)
Definition: FLA_Wilkshift_tridiag.c:105
FLA_Error FLA_Wilkshift_tridiag_opd(double delta1, double epsilon, double delta2, double *kappa)
Definition: FLA_Wilkshift_tridiag.c:155
FLA_Error FLA_Wilkshift_tridiag_check(FLA_Obj delta1, FLA_Obj epsilon, FLA_Obj delta2, FLA_Obj kappa)
Definition: FLA_Wilkshift_tridiag_check.c:13
unsigned int FLA_Check_error_level(void)
Definition: FLA_Check.c:18
FLA_Datatype FLA_Obj_datatype(FLA_Obj obj)
Definition: FLA_Query.c:13
int FLA_Datatype
Definition: FLA_type_defs.h:49

References FLA_Check_error_level(), FLA_Obj_datatype(), FLA_Wilkshift_tridiag_check(), FLA_Wilkshift_tridiag_opd(), and FLA_Wilkshift_tridiag_ops().

◆ FLA_Wilkshift_tridiag_opd()

FLA_Error FLA_Wilkshift_tridiag_opd ( double  delta1,
double  epsilon,
double  delta2,
double *  kappa 
)
159 {
160  double a = delta1;
161  double c = epsilon;
162  double d = delta2;
163  double p, q, r, s, k;
164 
165  // Begin with kappa equal to d.
166  k = d;
167 
168  // Compute a scaling factor to promote numerical stability.
169  s = fabs( a ) + 2.0 * fabs( c ) + fabs( d );
170 
171  if ( s == 0.0 ) return FLA_SUCCESS;
172 
173  // Compute q with scaling applied.
174  q = ( c / s ) * ( c / s );
175 
176  if ( q != 0.0 )
177  {
178 
179  // Compute p = 0.5*( a - d ) with scaling applied.
180  p = 0.5 * ( ( a / s ) - ( d / s ) );
181 
182  // Compute r = sqrt( p*p - q ).
183  r = sqrt( p * p + q );
184 
185  // If p*r is negative, then we need to negate r to ensure we use the
186  // larger of the two eigenvalues.
187  if ( p * r < 0.0 ) r = -r;
188 
189  // Compute the Wilkinson shift with scaling removed:
190  // k = lambda_min + d
191  // = d + lambda_min
192  // = d + (-q / lambda_max)
193  // = d - q / lambda_max
194  // = d - q / (p + r)
195  k = k - s * ( q / ( p + r ) );
196 
197 /*
198  // LAPACK method:
199 
200  p = 0.5 * ( ( a ) - ( d ) ) / c ;
201  //r = sqrt( p * p + 1.0 );
202  r = fla_dlapy2( p, 1.0 );
203  if ( p < 0.0 ) r = -r;
204  k = k - ( c / ( p + r ) );
205 */
206  }
207 
208  // Save the result.
209  *kappa = k;
210 
211  return FLA_SUCCESS;
212 }

Referenced by FLA_Tevd_eigval_n_opd_var1(), FLA_Tevd_eigval_v_opd_var1(), FLA_Tevd_eigval_v_opd_var3(), FLA_Tevd_find_perfshift_opd(), and FLA_Wilkshift_tridiag().

◆ FLA_Wilkshift_tridiag_ops()

FLA_Error FLA_Wilkshift_tridiag_ops ( float  delta1,
float  epsilon,
float  delta2,
float *  kappa 
)
109 {
110  float a = delta1;
111  float c = epsilon;
112  float d = delta2;
113  float p, q, r, s, k;
114 
115  // Begin with kappa equal to d.
116  k = d;
117 
118  // Compute a scaling factor to promote numerical stability.
119  s = fabs( a ) + 2.0F * fabs( c ) + fabs( d );
120 
121  if ( s == 0.0F ) return FLA_SUCCESS;
122 
123  // Compute q with scaling applied.
124  q = ( c / s ) * ( c / s );
125 
126  if ( q != 0.0F )
127  {
128  // Compute p = 0.5*( a - d ) with scaling applied.
129  p = 0.5F * ( ( a / s ) - ( d / s ) );
130 
131  // Compute r = sqrt( p*p - q ).
132  r = sqrt( p * p + q );
133 
134  // If p*r is negative, then we need to negate r to ensure we use the
135  // larger of the two eigenvalues.
136  if ( p * r < 0.0F ) r = -r;
137 
138  // Compute the Wilkinson shift with scaling removed:
139  // k = lambda_min + d
140  // = d + lambda_min
141  // = d + (-q / lambda_max)
142  // = d - q / lambda_max
143  // = d - q / (p + r)
144  k = k - s * ( q / ( p + r ) );
145  }
146 
147  // Save the result.
148  *kappa = k;
149 
150  return FLA_SUCCESS;
151 }

Referenced by FLA_Wilkshift_tridiag().