libflame  revision_anchor
Functions
FLA_Svv_2x2.c File Reference

(r)

Functions

FLA_Error FLA_Svv_2x2 (FLA_Obj alpha11, FLA_Obj alpha12, FLA_Obj alpha22, FLA_Obj sigma1, FLA_Obj sigma2, FLA_Obj gammaL, FLA_Obj sigmaL, FLA_Obj gammaR, FLA_Obj sigmaR)
 
FLA_Error FLA_Svv_2x2_ops (float *alpha11, float *alpha12, float *alpha22, float *sigma1, float *sigma2, float *gammaL, float *sigmaL, float *gammaR, float *sigmaR)
 
FLA_Error FLA_Svv_2x2_opd (double *alpha11, double *alpha12, double *alpha22, double *sigma1, double *sigma2, double *gammaL, double *sigmaL, double *gammaR, double *sigmaR)
 

Function Documentation

◆ FLA_Svv_2x2()

FLA_Error FLA_Svv_2x2 ( FLA_Obj  alpha11,
FLA_Obj  alpha12,
FLA_Obj  alpha22,
FLA_Obj  sigma1,
FLA_Obj  sigma2,
FLA_Obj  gammaL,
FLA_Obj  sigmaL,
FLA_Obj  gammaR,
FLA_Obj  sigmaR 
)
39 {
40  FLA_Datatype datatype;
41 
42  datatype = FLA_Obj_datatype( alpha11 );
43 
44  switch ( datatype )
45  {
46  case FLA_FLOAT:
47  {
48  float* buff_alpha11 = FLA_FLOAT_PTR( alpha11 );
49  float* buff_alpha12 = FLA_FLOAT_PTR( alpha12 );
50  float* buff_alpha22 = FLA_FLOAT_PTR( alpha22 );
51  float* buff_sigma1 = FLA_FLOAT_PTR( sigma1 );
52  float* buff_sigma2 = FLA_FLOAT_PTR( sigma2 );
53  float* buff_gammaL = FLA_FLOAT_PTR( gammaL );
54  float* buff_sigmaL = FLA_FLOAT_PTR( sigmaL );
55  float* buff_gammaR = FLA_FLOAT_PTR( gammaR );
56  float* buff_sigmaR = FLA_FLOAT_PTR( sigmaR );
57 
58  FLA_Svv_2x2_ops( buff_alpha11,
59  buff_alpha12,
60  buff_alpha22,
61  buff_sigma1,
62  buff_sigma2,
63  buff_gammaL,
64  buff_sigmaL,
65  buff_gammaR,
66  buff_sigmaR );
67 
68  break;
69  }
70 
71  case FLA_DOUBLE:
72  {
73  double* buff_alpha11 = FLA_DOUBLE_PTR( alpha11 );
74  double* buff_alpha12 = FLA_DOUBLE_PTR( alpha12 );
75  double* buff_alpha22 = FLA_DOUBLE_PTR( alpha22 );
76  double* buff_sigma1 = FLA_DOUBLE_PTR( sigma1 );
77  double* buff_sigma2 = FLA_DOUBLE_PTR( sigma2 );
78  double* buff_gammaL = FLA_DOUBLE_PTR( gammaL );
79  double* buff_sigmaL = FLA_DOUBLE_PTR( sigmaL );
80  double* buff_gammaR = FLA_DOUBLE_PTR( gammaR );
81  double* buff_sigmaR = FLA_DOUBLE_PTR( sigmaR );
82 
83  FLA_Svv_2x2_opd( buff_alpha11,
84  buff_alpha12,
85  buff_alpha22,
86  buff_sigma1,
87  buff_sigma2,
88  buff_gammaL,
89  buff_sigmaL,
90  buff_gammaR,
91  buff_sigmaR );
92 
93  break;
94  }
95  }
96 
97  return FLA_SUCCESS;
98 }
FLA_Error FLA_Svv_2x2_ops(float *alpha11, float *alpha12, float *alpha22, float *sigma1, float *sigma2, float *gammaL, float *sigmaL, float *gammaR, float *sigmaR)
Definition: FLA_Svv_2x2.c:102
FLA_Error FLA_Svv_2x2_opd(double *alpha11, double *alpha12, double *alpha22, double *sigma1, double *sigma2, double *gammaL, double *sigmaL, double *gammaR, double *sigmaR)
Definition: FLA_Svv_2x2.c:290
FLA_Datatype FLA_Obj_datatype(FLA_Obj obj)
Definition: FLA_Query.c:13
int FLA_Datatype
Definition: FLA_type_defs.h:49

References FLA_Obj_datatype(), FLA_Svv_2x2_opd(), and FLA_Svv_2x2_ops().

◆ FLA_Svv_2x2_opd()

FLA_Error FLA_Svv_2x2_opd ( double *  alpha11,
double *  alpha12,
double *  alpha22,
double *  sigma1,
double *  sigma2,
double *  gammaL,
double *  sigmaL,
double *  gammaR,
double *  sigmaR 
)
299 {
300  double zero = 0.0;
301  double half = 0.5;
302  double one = 1.0;
303  double two = 2.0;
304  double four = 4.0;
305 
306  double eps;
307 
308  double f, g, h;
309  double clt, crt, slt, srt;
310  double a, d, fa, ft, ga, gt, ha, ht, l;
311  double m, mm, r, s, t, temp, tsign, tt;
312  double ssmin, ssmax;
313  double csl, snl;
314  double csr, snr;
315 
316  int gasmal, swap;
317  int pmax;
318 
319  f = *alpha11;
320  g = *alpha12;
321  h = *alpha22;
322 
323  eps = FLA_Mach_params_opd( FLA_MACH_EPS );
324 
325  ft = f;
326  fa = fabs( f );
327  ht = h;
328  ha = fabs( h );
329 
330  // pmax points to the maximum absolute element of matrix.
331  // pmax = 1 if f largest in absolute values.
332  // pmax = 2 if g largest in absolute values.
333  // pmax = 3 if h largest in absolute values.
334 
335  pmax = 1;
336 
337  swap = ( ha > fa );
338  if ( swap )
339  {
340  pmax = 3;
341 
342  temp = ft;
343  ft = ht;
344  ht = temp;
345 
346  temp = fa;
347  fa = ha;
348  ha = temp;
349  }
350 
351  gt = g;
352  ga = fabs( g );
353 
354  if ( ga == zero )
355  {
356  // Diagonal matrix case.
357 
358  ssmin = ha;
359  ssmax = fa;
360  clt = one;
361  slt = zero;
362  crt = one;
363  srt = zero;
364  }
365  else
366  {
367  gasmal = TRUE;
368 
369  if ( ga > fa )
370  {
371  pmax = 2;
372 
373  if ( ( fa / ga ) < eps )
374  {
375  // Case of very large ga.
376 
377  gasmal = FALSE;
378 
379  ssmax = ga;
380 
381  if ( ha > one ) ssmin = fa / ( ga / ha );
382  else ssmin = ( fa / ga ) * ha;
383 
384  clt = one;
385  slt = ht / gt;
386  crt = ft / gt;
387  srt = one;
388  }
389  }
390 
391  if ( gasmal )
392  {
393  // Normal case.
394 
395  d = fa - ha;
396 
397  if ( d == fa ) l = one;
398  else l = d / fa;
399 
400  m = gt / ft;
401 
402  t = two - l;
403 
404  mm = m * m;
405  tt = t * t;
406  s = sqrt( tt + mm );
407 
408  if ( l == zero ) r = fabs( m );
409  else r = sqrt( l * l + mm );
410 
411  a = half * ( s + r );
412 
413  ssmin = ha / a;
414  ssmax = fa * a;
415 
416  if ( mm == zero )
417  {
418  // Here, m is tiny.
419 
420  if ( l == zero ) t = signof( two, ft ) * signof( one, gt );
421  else t = gt / signof( d, ft ) + m / t;
422  }
423  else
424  {
425  t = ( m / ( s + t ) + m / ( r + l ) ) * ( one + a );
426  }
427 
428  l = sqrt( t*t + four );
429  crt = two / l;
430  srt = t / l;
431  clt = ( crt + srt * m ) / a;
432  slt = ( ht / ft ) * srt / a;
433  }
434  }
435 
436  if ( swap )
437  {
438  csl = srt;
439  snl = crt;
440  csr = slt;
441  snr = clt;
442  }
443  else
444  {
445  csl = clt;
446  snl = slt;
447  csr = crt;
448  snr = srt;
449  }
450 
451 
452  // Correct the signs of ssmax and ssmin.
453 
454  if ( pmax == 1 )
455  tsign = signof( one, csr ) * signof( one, csl ) * signof( one, f );
456  else if ( pmax == 2 )
457  tsign = signof( one, snr ) * signof( one, csl ) * signof( one, g );
458  else // if ( pmax == 3 )
459  tsign = signof( one, snr ) * signof( one, snl ) * signof( one, h );
460 
461  ssmax = signof( ssmax, tsign );
462  ssmin = signof( ssmin, tsign * signof( one, f ) * signof( one, h ) );
463 
464  // Save the output values.
465 
466  *sigma1 = ssmin;
467  *sigma2 = ssmax;
468  *gammaL = csl;
469  *sigmaL = snl;
470  *gammaR = csr;
471  *sigmaR = snr;
472 
473  return FLA_SUCCESS;
474 }
double FLA_Mach_params_opd(FLA_Machval machval)
Definition: FLA_Mach_params.c:74
dcomplex temp
Definition: bl1_axpyv2b.c:301

References FLA_Mach_params_opd(), and temp.

Referenced by FLA_Bsvd_iteracc_v_opd_var1(), and FLA_Svv_2x2().

◆ FLA_Svv_2x2_ops()

FLA_Error FLA_Svv_2x2_ops ( float *  alpha11,
float *  alpha12,
float *  alpha22,
float *  sigma1,
float *  sigma2,
float *  gammaL,
float *  sigmaL,
float *  gammaR,
float *  sigmaR 
)
111 {
112  float zero = 0.0F;
113  float half = 0.5F;
114  float one = 1.0F;
115  float two = 2.0F;
116  float four = 4.0F;
117 
118  float eps;
119 
120  float f, g, h;
121  float clt, crt, slt, srt;
122  float a, d, fa, ft, ga, gt, ha, ht, l;
123  float m, mm, r, s, t, temp, tsign, tt;
124  float ssmin, ssmax;
125  float csl, snl;
126  float csr, snr;
127 
128  int gasmal, swap;
129  int pmax;
130 
131  f = *alpha11;
132  g = *alpha12;
133  h = *alpha22;
134 
135  eps = FLA_Mach_params_ops( FLA_MACH_EPS );
136 
137  ft = f;
138  fa = fabsf( f );
139  ht = h;
140  ha = fabsf( h );
141 
142  // pmax points to the maximum absolute element of matrix.
143  // pmax = 1 if f largest in absolute values.
144  // pmax = 2 if g largest in absolute values.
145  // pmax = 3 if h largest in absolute values.
146 
147  pmax = 1;
148 
149  swap = ( ha > fa );
150  if ( swap )
151  {
152  pmax = 3;
153 
154  temp = ft;
155  ft = ht;
156  ht = temp;
157 
158  temp = fa;
159  fa = ha;
160  ha = temp;
161  }
162 
163  gt = g;
164  ga = fabsf( g );
165 
166  if ( ga == zero )
167  {
168  // Diagonal matrix case.
169 
170  ssmin = ha;
171  ssmax = fa;
172  clt = one;
173  slt = zero;
174  crt = one;
175  srt = zero;
176  }
177  else
178  {
179  gasmal = TRUE;
180 
181  if ( ga > fa )
182  {
183  pmax = 2;
184 
185  if ( ( fa / ga ) < eps )
186  {
187  // Case of very large ga.
188 
189  gasmal = FALSE;
190 
191  ssmax = ga;
192 
193  if ( ha > one ) ssmin = fa / ( ga / ha );
194  else ssmin = ( fa / ga ) * ha;
195 
196  clt = one;
197  slt = ht / gt;
198  crt = ft / gt;
199  srt = one;
200  }
201  }
202 
203  if ( gasmal )
204  {
205  // Normal case.
206 
207  d = fa - ha;
208 
209  if ( d == fa ) l = one;
210  else l = d / fa;
211 
212  m = gt / ft;
213 
214  t = two - l;
215 
216  mm = m * m;
217  tt = t * t;
218  s = sqrtf( tt + mm );
219 
220  if ( l == zero ) r = fabsf( m );
221  else r = sqrtf( l * l + mm );
222 
223  a = half * ( s + r );
224 
225  ssmin = ha / a;
226  ssmax = fa * a;
227 
228  if ( mm == zero )
229  {
230  // Here, m is tiny.
231 
232  if ( l == zero ) t = signof( two, ft ) * signof( one, gt );
233  else t = gt / signof( d, ft ) + m / t;
234  }
235  else
236  {
237  t = ( m / ( s + t ) + m / ( r + l ) ) * ( one + a );
238  }
239 
240  l = sqrtf( t*t + four );
241  crt = two / l;
242  srt = t / l;
243  clt = ( crt + srt * m ) / a;
244  slt = ( ht / ft ) * srt / a;
245  }
246  }
247 
248  if ( swap )
249  {
250  csl = srt;
251  snl = crt;
252  csr = slt;
253  snr = clt;
254  }
255  else
256  {
257  csl = clt;
258  snl = slt;
259  csr = crt;
260  snr = srt;
261  }
262 
263 
264  // Correct the signs of ssmax and ssmin.
265 
266  if ( pmax == 1 )
267  tsign = signof( one, csr ) * signof( one, csl ) * signof( one, f );
268  else if ( pmax == 2 )
269  tsign = signof( one, snr ) * signof( one, csl ) * signof( one, g );
270  else // if ( pmax == 3 )
271  tsign = signof( one, snr ) * signof( one, snl ) * signof( one, h );
272 
273  ssmax = signof( ssmax, tsign );
274  ssmin = signof( ssmin, tsign * signof( one, f ) * signof( one, h ) );
275 
276  // Save the output values.
277 
278  *sigma1 = ssmin;
279  *sigma2 = ssmax;
280  *gammaL = csl;
281  *sigmaL = snl;
282  *gammaR = csr;
283  *sigmaR = snr;
284 
285  return FLA_SUCCESS;
286 }
float FLA_Mach_params_ops(FLA_Machval machval)
Definition: FLA_Mach_params.c:47

References FLA_Mach_params_ops(), and temp.

Referenced by FLA_Bsvd_iteracc_v_ops_var1(), and FLA_Svv_2x2().