QUDA  1.0.0
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
svd_quda.h
Go to the documentation of this file.
1 #include <float.h>
2 
3 #ifndef _SVD_QUDA_H_
4 #define _SVD_QUDA_H_
5 
6 #define DEVICEHOST __device__ __host__
7 #define SVDPREC 1e-11
8 #define LOG2 0.69314718055994530942
9 #define INVALID_DOUBLE (-DBL_MAX)
10 
11 
12 //namespace quda{
13 
14 // FIXME replace this with hypot
15 template<class Cmplx>
16 inline DEVICEHOST
17 typename std::remove_reference<decltype(Cmplx::x)>::type cabs(const Cmplx & z)
18 {
19  typedef typename std::remove_reference<decltype(Cmplx::x)>::type real;
20  real max, ratio, square;
21  if(fabs(z.x) > fabs(z.y)){ max = z.x; ratio = z.y/max; }else{ max=z.y; ratio = z.x/max; }
22  square = (max != 0.0) ? max*max*(1.0 + ratio*ratio) : 0.0;
23  return sqrt(square);
24 }
25 
26 
27 template<class T, class U>
28 inline DEVICEHOST typename PromoteTypeId<T,U>::Type quadSum(const T & a, const U & b){
29  typename PromoteTypeId<T,U>::Type ratio, square, max;
30  if (fabs(a) > fabs(b)) { max = a; ratio = b/a; } else { max=b; ratio = a/b; }
31  square = (max != 0.0) ? max*max*(1.0 + ratio*ratio) : 0.0;
32  return sqrt(square);
33 }
34 
35 
36 // In future iterations of the code, I would like to use templates to compute the norm
38 inline float getNorm(const Array<complex<float>,3>& a){
39  float temp1, temp2, temp3;
40  temp1 = cabs(a[0]);
41  temp2 = cabs(a[1]);
42  temp3 = quadSum(temp1,temp2);
43  temp1 = cabs(a[2]);
44  return quadSum(temp1,temp3);
45 }
46 
47 
49 inline double getNorm(const Array<complex<double>,3>& a){
50  double temp1, temp2, temp3;
51  temp1 = cabs(a[0]);
52  temp2 = cabs(a[1]);
53  temp3 = quadSum(temp1,temp2);
54  temp1 = cabs(a[2]);
55  return quadSum(temp1,temp3);
56 }
57 
58 
59 template<class T>
60 inline DEVICEHOST
61 void constructHHMat(const T & tau, const Array<T,3> & v, Matrix<T,3> & hh)
62 {
63  Matrix<T,3> temp1, temp2;
64  outerProd(v,v,&temp1);
65 
66  temp2 = conj(tau)*temp1;
67 
68  setIdentity(&temp1);
69 
70  hh = temp1 - temp2;
71  return;
72 }
73 
74 
75 template<class Real>
76 inline DEVICEHOST
77 void getLambdaMax(const Matrix<Real,3> & b, Real & lambda_max){
78 
79  Real m11 = b(1,1)*b(1,1) + b(0,1)*b(0,1);
80  Real m22 = b(2,2)*b(2,2) + b(1,2)*b(1,2);
81  Real m12 = b(1,1)*b(1,2);
82  Real dm = (m11 - m22) * 0.5;
83 
84  Real norm1 = quadSum(dm, m12);
85  if( dm >= 0.0 ){
86  lambda_max = m22 - (m12*m12)/(norm1 + dm);
87  }else{
88  lambda_max = m22 + (m12*m12)/(norm1 - dm);
89  }
90  return;
91 }
92 
93 
94 template<class Real>
95 inline DEVICEHOST
96 void getGivensRotation(const Real & alpha, const Real & beta, Real & c, Real & s)
97 {
98  Real ratio;
99  if( beta == 0.0 ){
100  c = 1.0; s = 0.0;
101  }else if( fabs(beta) > fabs(alpha) ){
102  ratio = -alpha/beta;
103  s = rsqrt(1.0 + ratio*ratio);
104  c = ratio*s;
105  }else{
106  ratio = -beta/alpha;
107  c = rsqrt(1.0 + ratio*ratio);
108  s = ratio*c;
109  }
110  return;
111 }
112 
113 template<class Real>
114 inline DEVICEHOST
115 void accumGivensRotation(int index, const Real & c, const Real & s, Matrix<Real,3> & m){
116  int index_p1 = index+1;
117  Real alpha, beta;
118  for(int i=0; i<3; ++i){
119  alpha = m(i,index);
120  beta = m(i,index_p1);
121  m(i,index) = alpha*c - beta*s;
122  m(i,index_p1) = alpha*s + beta*c;
123  }
124  return;
125 }
126 
127 template<class Real>
128 inline DEVICEHOST
129 void assignGivensRotation(const Real & c, const Real & s, Matrix<Real,2> & m){
130  m(0,0) = c;
131  m(0,1) = s;
132  m(1,0) = -s;
133  m(1,1) = c;
134  return;
135 }
136 
137 template<class Real>
138 inline DEVICEHOST
139 void swap(Real & a, Real & b){
140  Real temp = a;
141  a = b;
142  b = temp;
143  return;
144 }
145 
146 template<class Real>
147 inline DEVICEHOST
149  // set u and v to the 2x2 identity matrix
150  setIdentity(&u);
151  setIdentity(&v);
152 
153  Real c, s;
154 
155  if( m(0,0) == 0.0 ){
156 
157  getGivensRotation(m(0,1), m(1,1), c, s);
158 
159  m(0,0) = m(0,1)*c - m(1,1)*s;
160  m(0,1) = 0.0;
161  m(1,1) = 0.0;
162 
163  // exchange columns in v
164  v(0,0) = 0.0;
165  v(0,1) = 1.0;
166  v(1,0) = 1.0;
167  v(1,1) = 0.0;
168 
169  // u is a Givens rotation
170  assignGivensRotation(c, s, u);
171 
172  }else if( m(1,1) == 0.0 ){
173 
174  getGivensRotation(m(0,0), m(0,1), c, s);
175 
176  m(0,0) = m(0,0)*c - m(0,1)*s;
177  m(0,1) = 0.0;
178  m(1,1) = 0.0;
179 
180  assignGivensRotation(c,s,v);
181 
182  }else if( m(0,1) != 0.0 ){
183 
184  // need to calculate (m(1,1)**2 + m(0,1)**2 - m(0,0)**2)/(2*m(0,0)*m(0,1))
185  Real abs01 = fabs(m(0,1));
186  Real abs11 = fabs(m(1,1));
187  Real min, max;
188  if( abs01 > abs11 ){ min = abs11; max = abs01; }
189  else { min = abs01; max = abs11; }
190 
191 
192  Real ratio = min/max;
193  Real alpha = 2.0*log(max) + log(1.0 + ratio*ratio);
194 
195  Real abs00 = fabs(m(0,0));
196  Real beta = 2.0*log(abs00);
197 
198  int sign;
199  Real temp;
200 
201  if( alpha > beta ){
202  sign = 1;
203  temp = alpha + log(1.0 - exp(beta-alpha));
204  }else{
205  sign = -1;
206  temp = beta + log(1.0 - exp(alpha-beta));
207  }
208  temp -= LOG2 + log(abs00) + log(abs01);
209  temp = sign*exp(temp);
210 
211  if( m(0,0) < 0.0 ){ temp *= -1.0; }
212  if( m(0,1) < 0.0 ){ temp *= -1.0; }
213 
214  beta = (temp >= 0.0) ? quadSum(1.0, temp) : -quadSum(1.0, temp);
215  temp = 1.0/(temp + beta);
216 
217  // Calculate beta = sqrt(1 + temp**2)
218  beta = quadSum(1.0, temp);
219 
220  c = 1.0/beta;
221  s = temp*c;
222 
223  Matrix<Real,2> p;
224 
225  p(0,0) = c*m(0,0) - s*m(0,1);
226  p(1,0) = - s*m(1,1);
227  p(0,1) = s*m(0,0) + c*m(0,1);
228  p(1,1) = c*m(1,1);
229 
230  assignGivensRotation(c, s, v);
231 
232  // Make the column with the largest norm the first column
233  alpha = quadSum(p(0,0),p(1,0));
234  beta = quadSum(p(0,1),p(1,1));
235 
236  if(alpha < beta){
237  swap(p(0,0),p(0,1));
238  swap(p(1,0),p(1,1));
239 
240  swap(v(0,0),v(0,1));
241  swap(v(1,0),v(1,1));
242  }
243 
244  getGivensRotation(p(0,0), p(1,0), c, s);
245 
246  m(0,0) = p(0,0)*c - s*p(1,0);
247  m(1,1) = p(0,1)*s + c*p(1,1);
248  m(0,1) = 0.0;
249 
250  assignGivensRotation(c,s,u);
251  }
252 
253  return;
254 }
255 
256 // Change this so that the bidiagonal terms are not returned
257 // as a complex matrix
258 template<class Float>
259  DEVICEHOST
260 void getRealBidiagMatrix(const Matrix<complex<Float>,3> &mat, Matrix<complex<Float>,3> &u, Matrix<complex<Float>,3> &v)
261 {
262  typedef complex<Float> Complex;
263 
265  Array<Complex,3> vec;
266  Matrix<Complex,3> temp;
267 
268 
269  const Complex COMPLEX_UNITY(1.0,0.0);
270  const Complex COMPLEX_ZERO = 0.0;
271  // Step 1: build the first left reflector v,
272  // calculate the first left rotation
273  // apply to the original matrix
274  Float x = cabs(mat(1,0));
275  Float y = cabs(mat(2,0));
276  Float norm1 = quadSum(x,y);
277  Float beta;
278  Complex w, tau, z;
279 
280  // Set them to values that would never be set. If at the end of the
281  // below algorithmic flow, these values have been preserved then we
282  // know that the input matrix was the unit matrix
283  u(0,0).x = INVALID_DOUBLE;
284  v(0,0).x = INVALID_DOUBLE;
285 
286  if(norm1 == 0 && mat(0,0).y == 0){
287  p = mat;
288  }else{
289  Array<Complex,3> temp_vec;
290  copyColumn(mat, 0, &temp_vec);
291 
292  beta = (mat(0,0).x > 0.0) ? -getNorm(temp_vec) : getNorm(temp_vec);
293 
294  w.x = mat(0,0).x - beta; // work around for LLVM
295  w.y = mat(0,0).y;
296  Float norm1_inv = 1.0/cabs(w);
297  w = conj(w)*norm1_inv;
298 
299  // Assign the vector elements
300  vec[0] = COMPLEX_UNITY;
301  vec[1] = mat(1,0)*w*norm1_inv;
302  vec[2] = mat(2,0)*w*norm1_inv;
303 
304  // Now compute tau
305  tau.x = (beta - mat(0,0).x)/beta;
306  tau.y = mat(0,0).y/beta;
307  // construct the Householder matrix
308  constructHHMat(tau, vec, temp);
309  p = conj(temp)*mat;
310  u = temp;
311  }
312 
313  // Step 2: build the first right reflector
314  Float norm2 = cabs(p(0,2));
315  if(norm2 != 0.0 || p(0,1).y != 0.0){
316  norm1 = cabs(p(0,1));
317  beta = (p(0,1).x > 0.0) ? -quadSum(norm1,norm2) : quadSum(norm1,norm2);
318  vec[0] = COMPLEX_ZERO;
319  vec[1] = COMPLEX_UNITY;
320 
321  w.x = p(0,1).x-beta; // work around for LLVM
322  w.y = p(0,1).y;
323  Float norm1_inv = 1.0/cabs(w);
324  w = conj(w)*norm1_inv;
325  z = conj(p(0,2))*norm1_inv;
326  vec[2] = z*conj(w);
327 
328  tau.x = (beta - p(0,1).x)/beta; // work around for LLVM
329  tau.y = (- p(0,1).y)/beta;
330  // construct the Householder matrix
331  constructHHMat(tau, vec, temp);
332  p = p*temp;
333  v = temp;
334  }
335 
336  // Step 3: build the second left reflector
337  norm2 = cabs(p(2,1));
338  if(norm2 != 0.0 || p(1,1).y != 0.0){
339  norm1 = cabs(p(1,1));
340  beta = (p(1,1).x > 0) ? -quadSum(norm1,norm2) : quadSum(norm1,norm2);
341 
342  // Set the first two elements of the vector
343  vec[0] = COMPLEX_ZERO;
344  vec[1] = COMPLEX_UNITY;
345 
346  w.x = p(1,1).x - beta; // work around for LLVM
347  w.y = p(1,1).y;
348  Float norm1_inv = 1.0/cabs(w);
349  w = conj(w)*norm1_inv;
350  z = p(2,1)*norm1_inv;
351  vec[2] = z*w;
352 
353  tau.x = (beta - p(1,1).x)/beta;
354  tau.y = p(1,1).y/beta;
355  // I could very easily change the definition of tau
356  // so that we wouldn't need to call conj(temp) below.
357  // I would have to be careful to make sure this change
358  // is consistent with the other parts of the code.
359  constructHHMat(tau, vec, temp);
360  p = conj(temp)*p;
361  u = u*temp;
362  }
363 
364  // Step 4: build the second right reflector
365  setIdentity(&temp);
366  if( p(1,2).y != 0.0 ){
367  beta = p(1,2).x > 0.0 ? -cabs(p(1,2)) : cabs(p(1,2));
368  temp(2,2) = conj(p(1,2))/beta;
369  p(2,2) = p(2,2)*temp(2,2); // This is the only element of p needed below
370  v = v*temp;
371  }
372 
373  // Step 5: build the third left reflector
374  if( p(2,2).y != 0.0 ){
375  beta = p(2,2).x > 0.0 ? -cabs(p(2,2)) : cabs(p(2,2));
376  temp(2,2) = p(2,2)/beta;
377  u = u*temp;
378  }
379 
380  // unit matrix
381  if (u(0,0).x == INVALID_DOUBLE && v(0,0).x == INVALID_DOUBLE) {
382  u = mat;
383  v = mat;
384  }
385 
386  return;
387 }
388 
389 
390 template<class Real>
391  DEVICEHOST
392 void bdSVD(Matrix<Real,3>& u, Matrix<Real,3>& v, Matrix<Real,3>& b, int max_it)
393 {
394 
395  Real c,s;
396 
397  // set u and v matrices equal to the identity
398  setIdentity(&u);
399  setIdentity(&v);
400 
401  Real alpha, beta;
402 
403  int it=0;
404  do{
405 
406  if( fabs(b(0,1)) < SVDPREC*( fabs(b(0,0)) + fabs(b(1,1)) ) ){ b(0,1) = 0.0; }
407  if( fabs(b(1,2)) < SVDPREC*( fabs(b(0,0)) + fabs(b(2,2)) ) ){ b(1,2) = 0.0; }
408 
409  if( b(0,1) != 0.0 && b(1,2) != 0.0 ){
410  if( b(0,0) == 0.0 ){
411 
412  getGivensRotation(-b(0,1), b(1,1), s, c);
413 
414  for(int i=0; i<3; ++i){
415  alpha = u(i,0);
416  beta = u(i,1);
417  u(i,0) = alpha*c - beta*s;
418  u(i,1) = alpha*s + beta*c;
419  }
420 
421  b(1,1) = b(0,1)*s + b(1,1)*c;
422  b(0,1) = 0.0;
423 
424  b(0,2) = -b(1,2)*s;
425  b(1,2) *= c;
426 
427  getGivensRotation(-b(0,2), b(2,2), s, c);
428 
429  for(int i=0; i<3; ++i){
430  alpha = u(i,0);
431  beta = u(i,2);
432  u(i,0) = alpha*c - beta*s;
433  u(i,2) = alpha*s + beta*c;
434  }
435  b(2,2) = b(0,2)*s + b(2,2)*c;
436  b(0,2) = 0.0;
437 
438  }else if( b(1,1) == 0.0 ){
439  // same block
440  getGivensRotation(-b(1,2), b(2,2), s, c);
441  for(int i=0; i<3; ++i){
442  alpha = u(i,1);
443  beta = u(i,2);
444  u(i,1) = alpha*c - beta*s;
445  u(i,2) = alpha*s + beta*c;
446  }
447  b(2,2) = b(1,2)*s + b(2,2)*c;
448  b(1,2) = 0.0;
449  // end block
450  }else if( b(2,2) == 0.0 ){
451 
452  getGivensRotation(b(1,1), -b(1,2), c, s);
453  for(int i=0; i<3; ++i){
454  alpha = v(i,1);
455  beta = v(i,2);
456  v(i,1) = alpha*c + beta*s;
457  v(i,2) = -alpha*s + beta*c;
458  }
459  b(1,1) = b(1,1)*c + b(1,2)*s;
460  b(1,2) = 0.0;
461 
462  b(0,2) = -b(0,1)*s;
463  b(0,1) *= c;
464 
465  // apply second rotation, to remove b02
466  getGivensRotation(b(0,0), -b(0,2), c, s);
467  for(int i=0; i<3; ++i){
468  alpha = v(i,0);
469  beta = v(i,2);
470  v(i,0) = alpha*c + beta*s;
471  v(i,2) = -alpha*s + beta*c;
472  }
473  b(0,0) = b(0,0)*c + b(0,2)*s;
474  b(0,2) = 0.0;
475 
476  }else{ // Else entering normal QR iteration
477 
478  Real lambda_max;
479  getLambdaMax(b, lambda_max); // defined above
480 
481  alpha = b(0,0)*b(0,0) - lambda_max;
482  beta = b(0,0)*b(0,1);
483 
484  // c*beta + s*alpha = 0
485  getGivensRotation(alpha, beta, c, s);
486  // Multiply right v matrix
487  accumGivensRotation(0, c, s, v);
488  // apply to bidiagonal matrix, this generates b(1,0)
489  alpha = b(0,0);
490  beta = b(0,1);
491 
492  b(0,0) = alpha*c - beta*s;
493  b(0,1) = alpha*s + beta*c;
494  b(1,0) = -b(1,1)*s;
495  b(1,1) = b(1,1)*c;
496 
497  // Calculate the second Givens rotation (this time on the left)
498  getGivensRotation(b(0,0), b(1,0), c, s);
499 
500  // Multiply left u matrix
501  accumGivensRotation(0, c, s, u);
502 
503  b(0,0) = b(0,0)*c - b(1,0)*s;
504  alpha = b(0,1);
505  beta = b(1,1);
506  b(0,1) = alpha*c - beta*s;
507  b(1,1) = alpha*s + beta*c;
508  b(0,2) = -b(1,2)*s;
509  b(1,2) = b(1,2)*c;
510  // Following from the definition of the Givens rotation, b(1,0) should be equal to zero.
511  b(1,0) = 0.0;
512 
513  // Calculate the third Givens rotation (on the right)
514  getGivensRotation(b(0,1), b(0,2), c, s);
515 
516  // Multiply the right v matrix
517  accumGivensRotation(1, c, s, v);
518 
519  b(0,1) = b(0,1)*c - b(0,2)*s;
520  alpha = b(1,1);
521  beta = b(1,2);
522 
523  b(1,1) = alpha*c - beta*s;
524  b(1,2) = alpha*s + beta*c;
525  b(2,1) = -b(2,2)*s;
526  b(2,2) *= c;
527  b(0,2) = 0.0;
528 
529 
530  // Calculate the fourth Givens rotation (on the left)
531  getGivensRotation(b(1,1), b(2,1), c, s);
532  // Multiply left u matrix
533  accumGivensRotation(1, c, s, u);
534  // Eliminate b(2,1)
535  b(1,1) = b(1,1)*c - b(2,1)*s;
536  alpha = b(1,2);
537  beta = b(2,2);
538  b(1,2) = alpha*c - beta*s;
539  b(2,2) = alpha*s + beta*c;
540  b(2,1) = 0.0;
541 
542  } // end of normal QR iteration
543 
544 
545  }else if( b(0,1) == 0.0 ){
546  Matrix<Real,2> m_small, u_small, v_small;
547 
548  m_small(0,0) = b(1,1);
549  m_small(0,1) = b(1,2);
550  m_small(1,0) = b(2,1);
551  m_small(1,1) = b(2,2);
552 
553  smallSVD(u_small, v_small, m_small);
554 
555  b(1,1) = m_small(0,0);
556  b(1,2) = m_small(0,1);
557  b(2,1) = m_small(1,0);
558  b(2,2) = m_small(1,1);
559 
560 
561  for(int i=0; i<3; ++i){
562  alpha = u(i,1);
563  beta = u(i,2);
564  u(i,1) = alpha*u_small(0,0) + beta*u_small(1,0);
565  u(i,2) = alpha*u_small(0,1) + beta*u_small(1,1);
566 
567  alpha = v(i,1);
568  beta = v(i,2);
569  v(i,1) = alpha*v_small(0,0) + beta*v_small(1,0);
570  v(i,2) = alpha*v_small(0,1) + beta*v_small(1,1);
571  }
572 
573 
574  }else if( b(1,2) == 0.0 ){
575  Matrix<Real,2> m_small, u_small, v_small;
576 
577  m_small(0,0) = b(0,0);
578  m_small(0,1) = b(0,1);
579  m_small(1,0) = b(1,0);
580  m_small(1,1) = b(1,1);
581 
582  smallSVD(u_small, v_small, m_small);
583 
584  b(0,0) = m_small(0,0);
585  b(0,1) = m_small(0,1);
586  b(1,0) = m_small(1,0);
587  b(1,1) = m_small(1,1);
588 
589  for(int i=0; i<3; ++i){
590  alpha = u(i,0);
591  beta = u(i,1);
592  u(i,0) = alpha*u_small(0,0) + beta*u_small(1,0);
593  u(i,1) = alpha*u_small(0,1) + beta*u_small(1,1);
594 
595  alpha = v(i,0);
596  beta = v(i,1);
597  v(i,0) = alpha*v_small(0,0) + beta*v_small(1,0);
598  v(i,1) = alpha*v_small(0,1) + beta*v_small(1,1);
599  }
600 
601  } // end if b(1,2) == 0
602  it++;
603  } while( (b(0,1) != 0.0 || b(1,2) != 0.0) && it < max_it);
604 
605  for(int i=0; i<3; ++i){
606  if( b(i,i) < 0.0) {
607  b(i,i) *= -1;
608  for(int j=0; j<3; ++j){
609  v(j,i) *= -1;
610  }
611  }
612  }
613  return;
614 }
615 
616 
617 
618 template<class Float>
619  DEVICEHOST
620 void computeSVD(const Matrix<complex<Float>,3> & m, Matrix<complex<Float>,3>& u,
621  Matrix<complex<Float>,3>& v, Float singular_values[3])
622 {
623  getRealBidiagMatrix<Float>(m, u, v);
624 
625  Matrix<Float,3> bd, u_real, v_real;
626  // Make real
627  for(int i=0; i<3; ++i){
628  for(int j=0; j<3; ++j){
629  bd(i,j) = (conj(u)*m*(v))(i,j).real();
630  }
631  }
632 
633  bdSVD(u_real, v_real, bd, 500);
634  for(int i=0; i<3; ++i){
635  singular_values[i] = bd(i,i);
636  }
637 
638  u = u*u_real;
639  v = v*v_real;
640 
641  return;
642 }
643 
644 //} // end namespace quda
645 
646 #undef INVALID_DOUBLE
647 
648 #endif // _SVD_QUDA_H
Matrix< N, std::complex< T > > conj(const Matrix< N, std::complex< T > > &mat)
#define INVALID_DOUBLE
Definition: svd_quda.h:9
__host__ __device__ ValueType exp(ValueType x)
Definition: complex_quda.h:96
__host__ __device__ ValueType sqrt(ValueType x)
Definition: complex_quda.h:120
DEVICEHOST void getGivensRotation(const Real &alpha, const Real &beta, Real &c, Real &s)
Definition: svd_quda.h:96
DEVICEHOST void getLambdaMax(const Matrix< Real, 3 > &b, Real &lambda_max)
Definition: svd_quda.h:77
DEVICEHOST float getNorm(const Array< complex< float >, 3 > &a)
Definition: svd_quda.h:38
__device__ __host__ void outerProd(const Array< T, N > &a, const Array< T, N > &b, Matrix< T, N > *m)
Definition: quda_matrix.h:805
DEVICEHOST void computeSVD(const Matrix< complex< Float >, 3 > &m, Matrix< complex< Float >, 3 > &u, Matrix< complex< Float >, 3 > &v, Float singular_values[3])
Definition: svd_quda.h:620
DEVICEHOST void smallSVD(Matrix< Real, 2 > &u, Matrix< Real, 2 > &v, Matrix< Real, 2 > &m)
Definition: svd_quda.h:148
DEVICEHOST std::remove_reference< decltype(Cmplx::x)>::type cabs(const Cmplx &z)
Definition: svd_quda.h:17
#define DEVICEHOST
Definition: svd_quda.h:6
__device__ __host__ void copyColumn(const Matrix< T, N > &m, int c, Array< T, N > *a)
Definition: quda_matrix.h:793
double norm2(Float *v, int len)
std::complex< double > Complex
Definition: quda_internal.h:46
__device__ __host__ void setIdentity(Matrix< T, N > *m)
Definition: quda_matrix.h:653
DEVICEHOST void getRealBidiagMatrix(const Matrix< complex< Float >, 3 > &mat, Matrix< complex< Float >, 3 > &u, Matrix< complex< Float >, 3 > &v)
Definition: svd_quda.h:260
DEVICEHOST PromoteTypeId< T, U >::Type quadSum(const T &a, const U &b)
Definition: svd_quda.h:28
#define SVDPREC
Definition: svd_quda.h:7
__host__ __device__ ValueType log(ValueType x)
Definition: complex_quda.h:101
static int index(int ndim, const int *dims, const int *x)
Definition: comm_common.cpp:32
double norm1(const ColorSpinorField &b)
Definition: reduce_quda.cu:714
DEVICEHOST void swap(Real &a, Real &b)
Definition: svd_quda.h:139
DEVICEHOST void assignGivensRotation(const Real &c, const Real &s, Matrix< Real, 2 > &m)
Definition: svd_quda.h:129
__shared__ float s[]
DEVICEHOST void constructHHMat(const T &tau, const Array< T, 3 > &v, Matrix< T, 3 > &hh)
Definition: svd_quda.h:61
#define LOG2
Definition: svd_quda.h:8
DEVICEHOST void bdSVD(Matrix< Real, 3 > &u, Matrix< Real, 3 > &v, Matrix< Real, 3 > &b, int max_it)
Definition: svd_quda.h:392
void mat(void *out, void **link, void *in, int dagger_bit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision)
DEVICEHOST void accumGivensRotation(int index, const Real &c, const Real &s, Matrix< Real, 3 > &m)
Definition: svd_quda.h:115