MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
kernels.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_LINALG_KERNELS_HPP
13 #define MFEM_LINALG_KERNELS_HPP
14 
15 #include "../config/config.hpp"
16 #include "../general/backends.hpp"
17 #include "../general/globals.hpp"
18 
19 #include "matrix.hpp"
20 #include "tmatrix.hpp"
21 #include "tlayout.hpp"
22 #include "ttensor.hpp"
23 
24 // This header contains stand-alone functions for "small" dense linear algebra
25 // (at quadrature point or element-level) designed to be inlined directly into
26 // device kernels.
27 
28 // Many methods of the DenseMatrix class and some of the Vector class call these
29 // kernels directly on the host, see the implementations in linalg/densemat.cpp
30 // and linalg/vector.cpp.
31 
32 namespace mfem
33 {
34 
35 namespace kernels
36 {
37 
38 /// Compute the square of the Euclidean distance to another vector
39 template<int dim>
40 MFEM_HOST_DEVICE inline double DistanceSquared(const double *x, const double *y)
41 {
42  double d = 0.0;
43  for (int i = 0; i < dim; i++) { d += (x[i]-y[i])*(x[i]-y[i]); }
44  return d;
45 }
46 
47 /// Creates n x n diagonal matrix with diagonal elements c
48 template<int dim>
49 MFEM_HOST_DEVICE inline void Diag(const double c, double *data)
50 {
51  const int N = dim*dim;
52  for (int i = 0; i < N; i++) { data[i] = 0.0; }
53  for (int i = 0; i < dim; i++) { data[i*(dim+1)] = c; }
54 }
55 
56 /// Vector subtraction operation: z = a * (x - y)
57 template<int dim>
58 MFEM_HOST_DEVICE inline void Subtract(const double a,
59  const double *x, const double *y,
60  double *z)
61 {
62  for (int i = 0; i < dim; i++) { z[i] = a * (x[i] - y[i]); }
63 }
64 
65 /// Dense matrix operation: VWt += v w^t
66 template<int dim>
67 MFEM_HOST_DEVICE inline void AddMultVWt(const double *v, const double *w,
68  double *VWt)
69 {
70  for (int i = 0; i < dim; i++)
71  {
72  const double vi = v[i];
73  for (int j = 0; j < dim; j++) { VWt[i*dim+j] += vi * w[j]; }
74  }
75 }
76 
77 template<int H, int W, typename T>
78 MFEM_HOST_DEVICE inline
79 void FNorm(double &scale_factor, double &scaled_fnorm2, const T *data)
80 {
81  int i, hw = H * W;
82  T max_norm = 0.0, entry, fnorm2;
83 
84  for (i = 0; i < hw; i++)
85  {
86  entry = fabs(data[i]);
87  if (entry > max_norm)
88  {
89  max_norm = entry;
90  }
91  }
92 
93  if (max_norm == 0.0)
94  {
95  scale_factor = scaled_fnorm2 = 0.0;
96  return;
97  }
98 
99  fnorm2 = 0.0;
100  for (i = 0; i < hw; i++)
101  {
102  entry = data[i] / max_norm;
103  fnorm2 += entry * entry;
104  }
105 
106  scale_factor = max_norm;
107  scaled_fnorm2 = fnorm2;
108 }
109 
110 /// Compute the Frobenius norm of the matrix
111 template<int H, int W, typename T>
112 MFEM_HOST_DEVICE inline
113 double FNorm(const T *data)
114 {
115  double s, n2;
116  kernels::FNorm<H,W>(s, n2, data);
117  return s*sqrt(n2);
118 }
119 
120 /// Compute the square of the Frobenius norm of the matrix
121 template<int H, int W, typename T>
122 MFEM_HOST_DEVICE inline
123 double FNorm2(const T *data)
124 {
125  double s, n2;
126  kernels::FNorm<H,W>(s, n2, data);
127  return s*s*n2;
128 }
129 
130 /// Returns the l2 norm of the Vector with given @a size and @a data.
131 template<typename T>
132 MFEM_HOST_DEVICE inline
133 double Norml2(const int size, const T *data)
134 {
135  if (0 == size) { return 0.0; }
136  if (1 == size) { return std::abs(data[0]); }
137  T scale = 0.0;
138  T sum = 0.0;
139  for (int i = 0; i < size; i++)
140  {
141  if (data[i] != 0.0)
142  {
143  const T absdata = fabs(data[i]);
144  if (scale <= absdata)
145  {
146  const T sqr_arg = scale / absdata;
147  sum = 1.0 + sum * (sqr_arg * sqr_arg);
148  scale = absdata;
149  continue;
150  } // end if scale <= absdata
151  const T sqr_arg = absdata / scale;
152  sum += (sqr_arg * sqr_arg); // else scale > absdata
153  } // end if data[i] != 0
154  }
155  return scale * sqrt(sum);
156 }
157 
158 /** @brief Matrix vector multiplication: y = A x, where the matrix A is of size
159  @a height x @a width with given @a data, while @a x and @a y specify the
160  data of the input and output vectors. */
161 template<typename TA, typename TX, typename TY>
162 MFEM_HOST_DEVICE inline
163 void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
164 {
165  if (width == 0)
166  {
167  for (int row = 0; row < height; row++)
168  {
169  y[row] = 0.0;
170  }
171  return;
172  }
173  const TA *d_col = data;
174  TX x_col = x[0];
175  for (int row = 0; row < height; row++)
176  {
177  y[row] = x_col*d_col[row];
178  }
179  d_col += height;
180  for (int col = 1; col < width; col++)
181  {
182  x_col = x[col];
183  for (int row = 0; row < height; row++)
184  {
185  y[row] += x_col*d_col[row];
186  }
187  d_col += height;
188  }
189 }
190 
191 /** @brief Matrix transpose vector multiplication: y = At x, where the matrix A
192  is of size @a height x @a width with given @a data, while @a x and @a y
193  specify the data of the input and output vectors. */
194 template<typename TA, typename TX, typename TY>
195 MFEM_HOST_DEVICE inline
196 void MultTranspose(const int height, const int width, const TA *data,
197  const TX *x, TY *y)
198 {
199  if (height == 0)
200  {
201  for (int row = 0; row < width; row++)
202  {
203  y[row] = 0.0;
204  }
205  return;
206  }
207  TY *y_off = y;
208  for (int i = 0; i < width; ++i)
209  {
210  TY val = 0.0;
211  for (int j = 0; j < height; ++j)
212  {
213  val += x[j] * data[i * height + j];
214  }
215  *y_off = val;
216  y_off++;
217  }
218 }
219 
220 /// Symmetrize a square matrix with given @a size and @a data: A -> (A+A^T)/2.
221 template<typename T>
222 MFEM_HOST_DEVICE inline
223 void Symmetrize(const int size, T *data)
224 {
225  for (int i = 0; i < size; i++)
226  {
227  for (int j = 0; j < i; j++)
228  {
229  const T a = 0.5 * (data[i*size+j] + data[j*size+i]);
230  data[j*size+i] = data[i*size+j] = a;
231  }
232  }
233 }
234 
235 /// Compute the determinant of a square matrix of size dim with given @a data.
236 template<int dim, typename T>
237 MFEM_HOST_DEVICE inline T Det(const T *data)
238 {
239  return TDetHD<T>(ColumnMajorLayout2D<dim,dim>(), data);
240 }
241 
242 /** @brief Return the inverse of a matrix with given @a size and @a data into
243  the matrix with data @a inv_data. */
244 template<int dim, typename T>
245 MFEM_HOST_DEVICE inline
246 void CalcInverse(const T *data, T *inv_data)
247 {
248  typedef ColumnMajorLayout2D<dim,dim> layout_t;
249  const T det = TAdjDetHD<T>(layout_t(), data, layout_t(), inv_data);
250  TAssignHD<AssignOp::Mult>(layout_t(), inv_data, static_cast<T>(1.0)/det);
251 }
252 
253 /** @brief Return the adjugate of a matrix */
254 template<int dim, typename T>
255 MFEM_HOST_DEVICE inline
256 void CalcAdjugate(const T *data, T *adj_data)
257 {
258  typedef ColumnMajorLayout2D<dim,dim> layout_t;
259  TAdjugateHD<T>(layout_t(), data, layout_t(), adj_data);
260 }
261 
262 /** @brief Compute C = A + alpha*B, where the matrices A, B and C are of size @a
263  height x @a width with data @a Adata, @a Bdata and @a Cdata. */
264 template<typename TALPHA, typename TA, typename TB, typename TC>
265 MFEM_HOST_DEVICE inline
266 void Add(const int height, const int width, const TALPHA alpha,
267  const TA *Adata, const TB *Bdata, TC *Cdata)
268 {
269  for (int j = 0; j < width; j++)
270  {
271  for (int i = 0; i < height; i++)
272  {
273  const int n = i*width+j;
274  Cdata[n] = Adata[n] + alpha * Bdata[n];
275  }
276  }
277 }
278 
279 /** @brief Compute C = alpha*A + beta*B, where the matrices A, B and C are of
280  size @a height x @a width with data @a Adata, @a Bdata and @a Cdata. */
281 template<typename TALPHA, typename TBETA, typename TA, typename TB, typename TC>
282 MFEM_HOST_DEVICE inline
283 void Add(const int height, const int width,
284  const TALPHA alpha, const TA *Adata,
285  const TBETA beta, const TB *Bdata,
286  TC *Cdata)
287 {
288  const int m = height * width;
289  for (int i = 0; i < m; i++)
290  {
291  Cdata[i] = alpha * Adata[i] + beta * Bdata[i];
292  }
293 }
294 
295 /** @brief Compute B += A, where the matrices A and B are of size
296  @a height x @a width with data @a Adata and @a Bdata. */
297 template<typename TA, typename TB>
298 MFEM_HOST_DEVICE inline
299 void Add(const int height, const int width, const TA *Adata, TB *Bdata)
300 {
301  const int m = height * width;
302  for (int i = 0; i < m; i++)
303  {
304  Bdata[i] += Adata[i];
305  }
306 }
307 
308 /** @brief Compute B +=alpha*A, where the matrices A and B are of size
309  @a height x @a width with data @a Adata and @a Bdata. */
310 template<typename TA, typename TB>
311 MFEM_HOST_DEVICE inline
312 void Add(const int height, const int width,
313  const double alpha, const TA *Adata, TB *Bdata)
314 {
315  const int m = height * width;
316  for (int i = 0; i < m; i++)
317  {
318  Bdata[i] += alpha * Adata[i];
319  }
320 }
321 
322 /** @brief Compute B = alpha*A, where the matrices A and B are of size
323  @a height x @a width with data @a Adata and @a Bdata. */
324 template<typename TA, typename TB>
325 MFEM_HOST_DEVICE inline
326 void Set(const int height, const int width,
327  const double alpha, const TA *Adata, TB *Bdata)
328 {
329  const int m = height * width;
330  for (int i = 0; i < m; i++)
331  {
332  Bdata[i] = alpha * Adata[i];
333  }
334 }
335 
336 /** @brief Matrix-matrix multiplication: A = B * C, where the matrices A, B and
337  C are of sizes @a Aheight x @a Awidth, @a Aheight x @a Bwidth and @a Bwidth
338  x @a Awidth, respectively. */
339 template<typename TA, typename TB, typename TC>
340 MFEM_HOST_DEVICE inline
341 void Mult(const int Aheight, const int Awidth, const int Bwidth,
342  const TB *Bdata, const TC *Cdata, TA *Adata)
343 {
344  const int ah_x_aw = Aheight * Awidth;
345  for (int i = 0; i < ah_x_aw; i++) { Adata[i] = 0.0; }
346  for (int j = 0; j < Awidth; j++)
347  {
348  for (int k = 0; k < Bwidth; k++)
349  {
350  for (int i = 0; i < Aheight; i++)
351  {
352  Adata[i+j*Aheight] += Bdata[i+k*Aheight] * Cdata[k+j*Bwidth];
353  }
354  }
355  }
356 }
357 
358 /** @brief Multiply a matrix of size @a Aheight x @a Awidth and data @a Adata
359  with the transpose of a matrix of size @a Bheight x @a Awidth and data @a
360  Bdata: A * Bt. Return the result in a matrix with data @a ABtdata. */
361 template<typename TA, typename TB, typename TC>
362 MFEM_HOST_DEVICE inline
363 void MultABt(const int Aheight, const int Awidth, const int Bheight,
364  const TA *Adata, const TB *Bdata, TC *ABtdata)
365 {
366  const int ah_x_bh = Aheight * Bheight;
367  for (int i = 0; i < ah_x_bh; i++) { ABtdata[i] = 0.0; }
368  for (int k = 0; k < Awidth; k++)
369  {
370  TC *c = ABtdata;
371  for (int j = 0; j < Bheight; j++)
372  {
373  const double bjk = Bdata[j];
374  for (int i = 0; i < Aheight; i++)
375  {
376  c[i] += Adata[i] * bjk;
377  }
378  c += Aheight;
379  }
380  Adata += Aheight;
381  Bdata += Bheight;
382  }
383 }
384 
385 /** @brief Multiply the transpose of a matrix of size @a Aheight x @a Awidth
386  and data @a Adata with a matrix of size @a Aheight x @a Bwidth and data @a
387  Bdata: At * B. Return the result in a matrix with data @a AtBdata. */
388 template<typename TA, typename TB, typename TC>
389 MFEM_HOST_DEVICE inline
390 void MultAtB(const int Aheight, const int Awidth, const int Bwidth,
391  const TA *Adata, const TB *Bdata, TC *AtBdata)
392 {
393  TC *c = AtBdata;
394  for (int i = 0; i < Bwidth; ++i)
395  {
396  for (int j = 0; j < Awidth; ++j)
397  {
398  TC val = 0.0;
399  for (int k = 0; k < Aheight; ++k)
400  {
401  val += Adata[j * Aheight + k] * Bdata[i * Aheight + k];
402  }
403  *c = val;
404  c++;
405  }
406  }
407 }
408 
409 /// Compute the spectrum of the matrix of size dim with given @a data, returning
410 /// the eigenvalues in the array @a lambda and the eigenvectors in the array @a
411 /// vec (listed consecutively).
412 template<int dim> MFEM_HOST_DEVICE
413 void CalcEigenvalues(const double *data, double *lambda, double *vec);
414 
415 /// Return the i'th singular value of the matrix of size dim with given @a data.
416 template<int dim> MFEM_HOST_DEVICE
417 double CalcSingularvalue(const double *data, const int i);
418 
419 
420 // Utility functions for CalcEigenvalues and CalcSingularvalue
421 namespace internal
422 {
423 
424 /// Utility function to swap the values of @a a and @a b.
425 template<typename T>
426 MFEM_HOST_DEVICE static inline
427 void Swap(T &a, T &b)
428 {
429  T tmp = a;
430  a = b;
431  b = tmp;
432 }
433 
434 const double Epsilon = std::numeric_limits<double>::epsilon();
435 
436 /// Utility function used in CalcSingularvalue<3>.
437 MFEM_HOST_DEVICE static inline
438 void Eigenvalues2S(const double &d12, double &d1, double &d2)
439 {
440  const double sqrt_1_eps = sqrt(1./Epsilon);
441  if (d12 != 0.)
442  {
443  // "The Symmetric Eigenvalue Problem", B. N. Parlett, pp.189-190
444  double t;
445  const double zeta = (d2 - d1)/(2*d12); // inf/inf from overflows?
446  if (fabs(zeta) < sqrt_1_eps)
447  {
448  t = d12*copysign(1./(fabs(zeta) + sqrt(1. + zeta*zeta)), zeta);
449  }
450  else
451  {
452  t = d12*copysign(0.5/fabs(zeta), zeta);
453  }
454  d1 -= t;
455  d2 += t;
456  }
457 }
458 
459 /// Utility function used in CalcEigenvalues().
460 MFEM_HOST_DEVICE static inline
461 void Eigensystem2S(const double &d12, double &d1, double &d2,
462  double &c, double &s)
463 {
464  const double sqrt_1_eps = sqrt(1./Epsilon);
465  if (d12 == 0.0)
466  {
467  c = 1.;
468  s = 0.;
469  }
470  else
471  {
472  // "The Symmetric Eigenvalue Problem", B. N. Parlett, pp.189-190
473  double t;
474  const double zeta = (d2 - d1)/(2*d12);
475  const double azeta = fabs(zeta);
476  if (azeta < sqrt_1_eps)
477  {
478  t = copysign(1./(azeta + sqrt(1. + zeta*zeta)), zeta);
479  }
480  else
481  {
482  t = copysign(0.5/azeta, zeta);
483  }
484  c = sqrt(1./(1. + t*t));
485  s = c*t;
486  t *= d12;
487  d1 -= t;
488  d2 += t;
489  }
490 }
491 
492 
493 /// Utility function used in CalcEigenvalues<3>.
494 MFEM_HOST_DEVICE static inline
495 void GetScalingFactor(const double &d_max, double &mult)
496 {
497  int d_exp;
498  if (d_max > 0.)
499  {
500  mult = frexp(d_max, &d_exp);
501  if (d_exp == std::numeric_limits<double>::max_exponent)
502  {
503  mult *= std::numeric_limits<double>::radix;
504  }
505  mult = d_max/mult;
506  }
507  else
508  {
509  mult = 1.;
510  }
511  // mult = 2^d_exp is such that d_max/mult is in [0.5,1) or in other words
512  // d_max is in the interval [0.5,1)*mult
513 }
514 
515 /// Utility function used in CalcEigenvalues<3>.
516 MFEM_HOST_DEVICE static inline
517 bool KernelVector2G(const int &mode,
518  double &d1, double &d12, double &d21, double &d2)
519 {
520  // Find a vector (z1,z2) in the "near"-kernel of the matrix
521  // | d1 d12 |
522  // | d21 d2 |
523  // using QR factorization.
524  // The vector (z1,z2) is returned in (d1,d2). Return 'true' if the matrix
525  // is zero without setting (d1,d2).
526  // Note: in the current implementation |z1| + |z2| = 1.
527 
528  // l1-norms of the columns
529  double n1 = fabs(d1) + fabs(d21);
530  double n2 = fabs(d2) + fabs(d12);
531 
532  bool swap_columns = (n2 > n1);
533  double mu;
534 
535  if (!swap_columns)
536  {
537  if (n1 == 0.)
538  {
539  return true;
540  }
541 
542  if (mode == 0) // eliminate the larger entry in the column
543  {
544  if (fabs(d1) > fabs(d21))
545  {
546  Swap(d1, d21);
547  Swap(d12, d2);
548  }
549  }
550  else // eliminate the smaller entry in the column
551  {
552  if (fabs(d1) < fabs(d21))
553  {
554  Swap(d1, d21);
555  Swap(d12, d2);
556  }
557  }
558  }
559  else
560  {
561  // n2 > n1, swap columns 1 and 2
562  if (mode == 0) // eliminate the larger entry in the column
563  {
564  if (fabs(d12) > fabs(d2))
565  {
566  Swap(d1, d2);
567  Swap(d12, d21);
568  }
569  else
570  {
571  Swap(d1, d12);
572  Swap(d21, d2);
573  }
574  }
575  else // eliminate the smaller entry in the column
576  {
577  if (fabs(d12) < fabs(d2))
578  {
579  Swap(d1, d2);
580  Swap(d12, d21);
581  }
582  else
583  {
584  Swap(d1, d12);
585  Swap(d21, d2);
586  }
587  }
588  }
589 
590  n1 = hypot(d1, d21);
591 
592  if (d21 != 0.)
593  {
594  // v = (n1, n2)^t, |v| = 1
595  // Q = I - 2 v v^t, Q (d1, d21)^t = (mu, 0)^t
596  mu = copysign(n1, d1);
597  n1 = -d21*(d21/(d1 + mu)); // = d1 - mu
598  d1 = mu;
599  // normalize (n1,d21) to avoid overflow/underflow
600  // normalize (n1,d21) by the max-norm to avoid the sqrt call
601  if (fabs(n1) <= fabs(d21))
602  {
603  // (n1,n2) <-- (n1/d21,1)
604  n1 = n1/d21;
605  mu = (2./(1. + n1*n1))*(n1*d12 + d2);
606  d2 = d2 - mu;
607  d12 = d12 - mu*n1;
608  }
609  else
610  {
611  // (n1,n2) <-- (1,d21/n1)
612  n2 = d21/n1;
613  mu = (2./(1. + n2*n2))*(d12 + n2*d2);
614  d2 = d2 - mu*n2;
615  d12 = d12 - mu;
616  }
617  }
618 
619  // Solve:
620  // | d1 d12 | | z1 | = | 0 |
621  // | 0 d2 | | z2 | | 0 |
622 
623  // choose (z1,z2) to minimize |d1*z1 + d12*z2| + |d2*z2|
624  // under the condition |z1| + |z2| = 1, z2 >= 0 (for uniqueness)
625  // set t = z1, z2 = 1 - |t|, -1 <= t <= 1
626  // objective function is:
627  // |d1*t + d12*(1 - |t|)| + |d2|*(1 - |t|) -- piecewise linear with
628  // possible minima are -1,0,1,t1 where t1: d1*t1 + d12*(1 - |t1|) = 0
629  // values: @t=+/-1 -> |d1|, @t=0 -> |n1| + |d2|, @t=t1 -> |d2|*(1 - |t1|)
630 
631  // evaluate z2 @t=t1
632  mu = -d12/d1;
633  // note: |mu| <= 1, if using l2-norm for column pivoting
634  // |mu| <= sqrt(2), if using l1-norm
635  n2 = 1./(1. + fabs(mu));
636  // check if |d1|<=|d2|*z2
637  if (fabs(d1) <= n2*fabs(d2))
638  {
639  d2 = 0.;
640  d1 = 1.;
641  }
642  else
643  {
644  d2 = n2;
645  // d1 = (n2 < 0.5) ? copysign(1. - n2, mu) : mu*n2;
646  d1 = mu*n2;
647  }
648 
649  if (swap_columns)
650  {
651  Swap(d1, d2);
652  }
653 
654  return false;
655 }
656 
657 /// Utility function used in CalcEigenvalues<3>.
658 MFEM_HOST_DEVICE static inline
659 void Vec_normalize3_aux(const double &x1, const double &x2,
660  const double &x3,
661  double &n1, double &n2, double &n3)
662 {
663  double t, r;
664 
665  const double m = fabs(x1);
666  r = x2/m;
667  t = 1. + r*r;
668  r = x3/m;
669  t = sqrt(1./(t + r*r));
670  n1 = copysign(t, x1);
671  t /= m;
672  n2 = x2*t;
673  n3 = x3*t;
674 }
675 
676 /// Utility function used in CalcEigenvalues<3>.
677 MFEM_HOST_DEVICE static inline
678 void Vec_normalize3(const double &x1, const double &x2, const double &x3,
679  double &n1, double &n2, double &n3)
680 {
681  // should work ok when xk is the same as nk for some or all k
682  if (fabs(x1) >= fabs(x2))
683  {
684  if (fabs(x1) >= fabs(x3))
685  {
686  if (x1 != 0.)
687  {
688  Vec_normalize3_aux(x1, x2, x3, n1, n2, n3);
689  }
690  else
691  {
692  n1 = n2 = n3 = 0.;
693  }
694  return;
695  }
696  }
697  else if (fabs(x2) >= fabs(x3))
698  {
699  Vec_normalize3_aux(x2, x1, x3, n2, n1, n3);
700  return;
701  }
702  Vec_normalize3_aux(x3, x1, x2, n3, n1, n2);
703 }
704 
705 /// Utility function used in CalcEigenvalues<3>.
706 MFEM_HOST_DEVICE static inline
707 int KernelVector3G_aux(const int &mode,
708  double &d1, double &d2, double &d3,
709  double &c12, double &c13, double &c23,
710  double &c21, double &c31, double &c32)
711 {
712  int kdim;
713  double mu, n1, n2, n3, s1, s2, s3;
714 
715  s1 = hypot(c21, c31);
716  n1 = hypot(d1, s1);
717 
718  if (s1 != 0.)
719  {
720  // v = (s1, s2, s3)^t, |v| = 1
721  // Q = I - 2 v v^t, Q (d1, c12, c13)^t = (mu, 0, 0)^t
722  mu = copysign(n1, d1);
723  n1 = -s1*(s1/(d1 + mu)); // = d1 - mu
724  d1 = mu;
725 
726  // normalize (n1,c21,c31) to avoid overflow/underflow
727  // normalize (n1,c21,c31) by the max-norm to avoid the sqrt call
728  if (fabs(n1) >= fabs(c21))
729  {
730  if (fabs(n1) >= fabs(c31))
731  {
732  // n1 is max, (s1,s2,s3) <-- (1,c21/n1,c31/n1)
733  s2 = c21/n1;
734  s3 = c31/n1;
735  mu = 2./(1. + s2*s2 + s3*s3);
736  n2 = mu*(c12 + s2*d2 + s3*c32);
737  n3 = mu*(c13 + s2*c23 + s3*d3);
738  c12 = c12 - n2;
739  d2 = d2 - s2*n2;
740  c32 = c32 - s3*n2;
741  c13 = c13 - n3;
742  c23 = c23 - s2*n3;
743  d3 = d3 - s3*n3;
744  goto done_column_1;
745  }
746  }
747  else if (fabs(c21) >= fabs(c31))
748  {
749  // c21 is max, (s1,s2,s3) <-- (n1/c21,1,c31/c21)
750  s1 = n1/c21;
751  s3 = c31/c21;
752  mu = 2./(1. + s1*s1 + s3*s3);
753  n2 = mu*(s1*c12 + d2 + s3*c32);
754  n3 = mu*(s1*c13 + c23 + s3*d3);
755  c12 = c12 - s1*n2;
756  d2 = d2 - n2;
757  c32 = c32 - s3*n2;
758  c13 = c13 - s1*n3;
759  c23 = c23 - n3;
760  d3 = d3 - s3*n3;
761  goto done_column_1;
762  }
763  // c31 is max, (s1,s2,s3) <-- (n1/c31,c21/c31,1)
764  s1 = n1/c31;
765  s2 = c21/c31;
766  mu = 2./(1. + s1*s1 + s2*s2);
767  n2 = mu*(s1*c12 + s2*d2 + c32);
768  n3 = mu*(s1*c13 + s2*c23 + d3);
769  c12 = c12 - s1*n2;
770  d2 = d2 - s2*n2;
771  c32 = c32 - n2;
772  c13 = c13 - s1*n3;
773  c23 = c23 - s2*n3;
774  d3 = d3 - n3;
775  }
776 
777 done_column_1:
778 
779  // Solve:
780  // | d2 c23 | | z2 | = | 0 |
781  // | c32 d3 | | z3 | | 0 |
782  if (KernelVector2G(mode, d2, c23, c32, d3))
783  {
784  // Have two solutions:
785  // two vectors in the kernel are P (-c12/d1, 1, 0)^t and
786  // P (-c13/d1, 0, 1)^t where P is the permutation matrix swapping
787  // entries 1 and col.
788 
789  // A vector orthogonal to both these vectors is P (1, c12/d1, c13/d1)^t
790  d2 = c12/d1;
791  d3 = c13/d1;
792  d1 = 1.;
793  kdim = 2;
794  }
795  else
796  {
797  // solve for z1:
798  // note: |z1| <= a since |z2| + |z3| = 1, and
799  // max{|c12|,|c13|} <= max{norm(col. 2),norm(col. 3)}
800  // <= norm(col. 1) <= a |d1|
801  // a = 1, if using l2-norm for column pivoting
802  // a = sqrt(3), if using l1-norm
803  d1 = -(c12*d2 + c13*d3)/d1;
804  kdim = 1;
805  }
806 
807  Vec_normalize3(d1, d2, d3, d1, d2, d3);
808 
809  return kdim;
810 }
811 
812 /// Utility function used in CalcEigenvalues<3>.
813 MFEM_HOST_DEVICE static inline
814 int KernelVector3S(const int &mode, const double &d12,
815  const double &d13, const double &d23,
816  double &d1, double &d2, double &d3)
817 {
818  // Find a unit vector (z1,z2,z3) in the "near"-kernel of the matrix
819  // | d1 d12 d13 |
820  // | d12 d2 d23 |
821  // | d13 d23 d3 |
822  // using QR factorization.
823  // The vector (z1,z2,z3) is returned in (d1,d2,d3).
824  // Returns the dimension of the kernel, kdim, but never zero.
825  // - if kdim == 3, then (d1,d2,d3) is not defined,
826  // - if kdim == 2, then (d1,d2,d3) is a vector orthogonal to the kernel,
827  // - otherwise kdim == 1 and (d1,d2,d3) is a vector in the "near"-kernel.
828 
829  double c12 = d12, c13 = d13, c23 = d23;
830  double c21, c31, c32;
831  int col, row;
832 
833  // l1-norms of the columns:
834  c32 = fabs(d1) + fabs(c12) + fabs(c13);
835  c31 = fabs(d2) + fabs(c12) + fabs(c23);
836  c21 = fabs(d3) + fabs(c13) + fabs(c23);
837 
838  // column pivoting: choose the column with the largest norm
839  if (c32 >= c21)
840  {
841  col = (c32 >= c31) ? 1 : 2;
842  }
843  else
844  {
845  col = (c31 >= c21) ? 2 : 3;
846  }
847  switch (col)
848  {
849  case 1:
850  if (c32 == 0.) // zero matrix
851  {
852  return 3;
853  }
854  break;
855 
856  case 2:
857  if (c31 == 0.) // zero matrix
858  {
859  return 3;
860  }
861  Swap(c13, c23);
862  Swap(d1, d2);
863  break;
864 
865  case 3:
866  if (c21 == 0.) // zero matrix
867  {
868  return 3;
869  }
870  Swap(c12, c23);
871  Swap(d1, d3);
872  }
873 
874  // row pivoting depending on 'mode'
875  if (mode == 0)
876  {
877  if (fabs(d1) <= fabs(c13))
878  {
879  row = (fabs(d1) <= fabs(c12)) ? 1 : 2;
880  }
881  else
882  {
883  row = (fabs(c12) <= fabs(c13)) ? 2 : 3;
884  }
885  }
886  else
887  {
888  if (fabs(d1) >= fabs(c13))
889  {
890  row = (fabs(d1) >= fabs(c12)) ? 1 : 2;
891  }
892  else
893  {
894  row = (fabs(c12) >= fabs(c13)) ? 2 : 3;
895  }
896  }
897  switch (row)
898  {
899  case 1:
900  c21 = c12;
901  c31 = c13;
902  c32 = c23;
903  break;
904 
905  case 2:
906  c21 = d1;
907  c31 = c13;
908  c32 = c23;
909  d1 = c12;
910  c12 = d2;
911  d2 = d1;
912  c13 = c23;
913  c23 = c31;
914  break;
915 
916  case 3:
917  c21 = c12;
918  c31 = d1;
919  c32 = c12;
920  d1 = c13;
921  c12 = c23;
922  c13 = d3;
923  d3 = d1;
924  }
925  row = KernelVector3G_aux(mode, d1, d2, d3, c12, c13, c23, c21, c31, c32);
926  // row is kdim
927 
928  switch (col)
929  {
930  case 2:
931  Swap(d1, d2);
932  break;
933 
934  case 3:
935  Swap(d1, d3);
936  }
937  return row;
938 }
939 
940 /// Utility function used in CalcEigenvalues<3>.
941 MFEM_HOST_DEVICE static inline
942 int Reduce3S(const int &mode,
943  double &d1, double &d2, double &d3,
944  double &d12, double &d13, double &d23,
945  double &z1, double &z2, double &z3,
946  double &v1, double &v2, double &v3,
947  double &g)
948 {
949  // Given the matrix
950  // | d1 d12 d13 |
951  // A = | d12 d2 d23 |
952  // | d13 d23 d3 |
953  // and a unit eigenvector z=(z1,z2,z3), transform the matrix A into the
954  // matrix B = Q P A P Q that has the form
955  // | b1 0 0 |
956  // B = Q P A P Q = | 0 b2 b23 |
957  // | 0 b23 b3 |
958  // where P is the permutation matrix switching entries 1 and k, and
959  // Q is the reflection matrix Q = I - g v v^t, defined by: set y = P z and
960  // v = c(y - e_1); if y = e_1, then v = 0 and Q = I.
961  // Note: Q y = e_1, Q e_1 = y ==> Q P A P Q e_1 = ... = lambda e_1.
962  // The entries (b1,b2,b3,b23) are returned in (d1,d2,d3,d23), and the
963  // return value of the function is k. The variable g = 2/(v1^2+v2^2+v3^3).
964 
965  int k;
966  double s, w1, w2, w3;
967 
968  if (mode == 0)
969  {
970  // choose k such that z^t e_k = zk has the smallest absolute value, i.e.
971  // the angle between z and e_k is closest to pi/2
972  if (fabs(z1) <= fabs(z3))
973  {
974  k = (fabs(z1) <= fabs(z2)) ? 1 : 2;
975  }
976  else
977  {
978  k = (fabs(z2) <= fabs(z3)) ? 2 : 3;
979  }
980  }
981  else
982  {
983  // choose k such that zk is the largest by absolute value
984  if (fabs(z1) >= fabs(z3))
985  {
986  k = (fabs(z1) >= fabs(z2)) ? 1 : 2;
987  }
988  else
989  {
990  k = (fabs(z2) >= fabs(z3)) ? 2 : 3;
991  }
992  }
993  switch (k)
994  {
995  case 2:
996  Swap(d13, d23);
997  Swap(d1, d2);
998  Swap(z1, z2);
999  break;
1000 
1001  case 3:
1002  Swap(d12, d23);
1003  Swap(d1, d3);
1004  Swap(z1, z3);
1005  }
1006 
1007  s = hypot(z2, z3);
1008 
1009  if (s == 0.)
1010  {
1011  // s can not be zero, if zk is the smallest (mode == 0)
1012  v1 = v2 = v3 = 0.;
1013  g = 1.;
1014  }
1015  else
1016  {
1017  g = copysign(1., z1);
1018  v1 = -s*(s/(z1 + g)); // = z1 - g
1019  // normalize (v1,z2,z3) by its max-norm, avoiding the sqrt call
1020  g = fabs(v1);
1021  if (fabs(z2) > g) { g = fabs(z2); }
1022  if (fabs(z3) > g) { g = fabs(z3); }
1023  v1 = v1/g;
1024  v2 = z2/g;
1025  v3 = z3/g;
1026  g = 2./(v1*v1 + v2*v2 + v3*v3);
1027 
1028  // Compute Q A Q = A - v w^t - w v^t, where
1029  // w = u - (g/2)(v^t u) v, and u = g A v
1030  // set w = g A v
1031  w1 = g*( d1*v1 + d12*v2 + d13*v3);
1032  w2 = g*(d12*v1 + d2*v2 + d23*v3);
1033  w3 = g*(d13*v1 + d23*v2 + d3*v3);
1034  // w := w - (g/2)(v^t w) v
1035  s = (g/2)*(v1*w1 + v2*w2 + v3*w3);
1036  w1 -= s*v1;
1037  w2 -= s*v2;
1038  w3 -= s*v3;
1039  // dij -= vi*wj + wi*vj
1040  d1 -= 2*v1*w1;
1041  d2 -= 2*v2*w2;
1042  d23 -= v2*w3 + v3*w2;
1043  d3 -= 2*v3*w3;
1044  // compute the off-diagonal entries on the first row/column of B which
1045  // should be zero (for debugging):
1046 #if 0
1047  s = d12 - v1*w2 - v2*w1; // b12 = 0
1048  s = d13 - v1*w3 - v3*w1; // b13 = 0
1049 #endif
1050  }
1051 
1052  switch (k)
1053  {
1054  case 2:
1055  Swap(z1, z2);
1056  break;
1057  case 3:
1058  Swap(z1, z3);
1059  }
1060  return k;
1061 }
1062 
1063 } // namespace kernels::internal
1064 
1065 
1066 // Implementations of CalcEigenvalues and CalcSingularvalue for dim = 2, 3.
1067 
1068 /// Compute the spectrum of the matrix of size 2 with given @a data, returning
1069 /// the eigenvalues in the array @a lambda and the eigenvectors in the array @a
1070 /// vec (listed consecutively).
1071 template<> MFEM_HOST_DEVICE inline
1072 void CalcEigenvalues<2>(const double *data, double *lambda, double *vec)
1073 {
1074  double d0 = data[0];
1075  double d2 = data[2]; // use the upper triangular entry
1076  double d3 = data[3];
1077  double c, s;
1078  internal::Eigensystem2S(d2, d0, d3, c, s);
1079  if (d0 <= d3)
1080  {
1081  lambda[0] = d0;
1082  lambda[1] = d3;
1083  vec[0] = c;
1084  vec[1] = -s;
1085  vec[2] = s;
1086  vec[3] = c;
1087  }
1088  else
1089  {
1090  lambda[0] = d3;
1091  lambda[1] = d0;
1092  vec[0] = s;
1093  vec[1] = c;
1094  vec[2] = c;
1095  vec[3] = -s;
1096  }
1097 }
1098 
1099 /// Compute the spectrum of the matrix of size 3 with given @a data, returning
1100 /// the eigenvalues in the array @a lambda and the eigenvectors in the array @a
1101 /// vec (listed consecutively).
1102 template<> MFEM_HOST_DEVICE inline
1103 void CalcEigenvalues<3>(const double *data, double *lambda, double *vec)
1104 {
1105  double d11 = data[0];
1106  double d12 = data[3]; // use the upper triangular entries
1107  double d22 = data[4];
1108  double d13 = data[6];
1109  double d23 = data[7];
1110  double d33 = data[8];
1111 
1112  double mult;
1113  {
1114  double d_max = fabs(d11);
1115  if (d_max < fabs(d22)) { d_max = fabs(d22); }
1116  if (d_max < fabs(d33)) { d_max = fabs(d33); }
1117  if (d_max < fabs(d12)) { d_max = fabs(d12); }
1118  if (d_max < fabs(d13)) { d_max = fabs(d13); }
1119  if (d_max < fabs(d23)) { d_max = fabs(d23); }
1120 
1121  internal::GetScalingFactor(d_max, mult);
1122  }
1123 
1124  d11 /= mult; d22 /= mult; d33 /= mult;
1125  d12 /= mult; d13 /= mult; d23 /= mult;
1126 
1127  double aa = (d11 + d22 + d33)/3; // aa = tr(A)/3
1128  double c1 = d11 - aa;
1129  double c2 = d22 - aa;
1130  double c3 = d33 - aa;
1131 
1132  double Q, R;
1133 
1134  Q = (2*(d12*d12 + d13*d13 + d23*d23) + c1*c1 + c2*c2 + c3*c3)/6;
1135  R = (c1*(d23*d23 - c2*c3)+ d12*(d12*c3 - 2*d13*d23) + d13*d13*c2)/2;
1136 
1137  if (Q <= 0.)
1138  {
1139  lambda[0] = lambda[1] = lambda[2] = aa;
1140  vec[0] = 1.; vec[3] = 0.; vec[6] = 0.;
1141  vec[1] = 0.; vec[4] = 1.; vec[7] = 0.;
1142  vec[2] = 0.; vec[5] = 0.; vec[8] = 1.;
1143  }
1144  else
1145  {
1146  double sqrtQ = sqrt(Q);
1147  double sqrtQ3 = Q*sqrtQ;
1148  // double sqrtQ3 = sqrtQ*sqrtQ*sqrtQ;
1149  // double sqrtQ3 = pow(Q, 1.5);
1150  double r;
1151  if (fabs(R) >= sqrtQ3)
1152  {
1153  if (R < 0.)
1154  {
1155  // R = -1.;
1156  r = 2*sqrtQ;
1157  }
1158  else
1159  {
1160  // R = 1.;
1161  r = -2*sqrtQ;
1162  }
1163  }
1164  else
1165  {
1166  R = R/sqrtQ3;
1167 
1168  if (R < 0.)
1169  {
1170  r = -2*sqrtQ*cos((acos(R) + 2.0*M_PI)/3); // max
1171  }
1172  else
1173  {
1174  r = -2*sqrtQ*cos(acos(R)/3); // min
1175  }
1176  }
1177 
1178  aa += r;
1179  c1 = d11 - aa;
1180  c2 = d22 - aa;
1181  c3 = d33 - aa;
1182 
1183  // Type of Householder reflections: z --> mu ek, where k is the index
1184  // of the entry in z with:
1185  // mode == 0: smallest absolute value --> angle closest to pi/2
1186  // mode == 1: largest absolute value --> angle farthest from pi/2
1187  // Observations:
1188  // mode == 0 produces better eigenvectors, less accurate eigenvalues?
1189  // mode == 1 produces better eigenvalues, less accurate eigenvectors?
1190  const int mode = 0;
1191 
1192  // Find a unit vector z = (z1,z2,z3) in the "near"-kernel of
1193  // | c1 d12 d13 |
1194  // | d12 c2 d23 | = A - aa*I
1195  // | d13 d23 c3 |
1196  // This vector is also an eigenvector for A corresponding to aa.
1197  // The vector z overwrites (c1,c2,c3).
1198  switch (internal::KernelVector3S(mode, d12, d13, d23, c1, c2, c3))
1199  {
1200  case 3:
1201  // 'aa' is a triple eigenvalue
1202  lambda[0] = lambda[1] = lambda[2] = aa;
1203  vec[0] = 1.; vec[3] = 0.; vec[6] = 0.;
1204  vec[1] = 0.; vec[4] = 1.; vec[7] = 0.;
1205  vec[2] = 0.; vec[5] = 0.; vec[8] = 1.;
1206  goto done_3d;
1207 
1208  case 2:
1209  // ok, continue with the returned vector orthogonal to the kernel
1210  case 1:
1211  // ok, continue with the returned vector in the "near"-kernel
1212  ;
1213  }
1214 
1215  // Using the eigenvector c=(c1,c2,c3) transform A into
1216  // | d11 0 0 |
1217  // A <-- Q P A P Q = | 0 d22 d23 |
1218  // | 0 d23 d33 |
1219  double v1, v2, v3, g;
1220  int k = internal::Reduce3S(mode, d11, d22, d33, d12, d13, d23,
1221  c1, c2, c3, v1, v2, v3, g);
1222  // Q = I - 2 v v^t
1223  // P - permutation matrix switching entries 1 and k
1224 
1225  // find the eigenvalues and eigenvectors for
1226  // | d22 d23 |
1227  // | d23 d33 |
1228  double c, s;
1229  internal::Eigensystem2S(d23, d22, d33, c, s);
1230  // d22 <-> P Q (0, c, -s), d33 <-> P Q (0, s, c)
1231 
1232  double *vec_1, *vec_2, *vec_3;
1233  if (d11 <= d22)
1234  {
1235  if (d22 <= d33)
1236  {
1237  lambda[0] = d11; vec_1 = vec;
1238  lambda[1] = d22; vec_2 = vec + 3;
1239  lambda[2] = d33; vec_3 = vec + 6;
1240  }
1241  else if (d11 <= d33)
1242  {
1243  lambda[0] = d11; vec_1 = vec;
1244  lambda[1] = d33; vec_3 = vec + 3;
1245  lambda[2] = d22; vec_2 = vec + 6;
1246  }
1247  else
1248  {
1249  lambda[0] = d33; vec_3 = vec;
1250  lambda[1] = d11; vec_1 = vec + 3;
1251  lambda[2] = d22; vec_2 = vec + 6;
1252  }
1253  }
1254  else
1255  {
1256  if (d11 <= d33)
1257  {
1258  lambda[0] = d22; vec_2 = vec;
1259  lambda[1] = d11; vec_1 = vec + 3;
1260  lambda[2] = d33; vec_3 = vec + 6;
1261  }
1262  else if (d22 <= d33)
1263  {
1264  lambda[0] = d22; vec_2 = vec;
1265  lambda[1] = d33; vec_3 = vec + 3;
1266  lambda[2] = d11; vec_1 = vec + 6;
1267  }
1268  else
1269  {
1270  lambda[0] = d33; vec_3 = vec;
1271  lambda[1] = d22; vec_2 = vec + 3;
1272  lambda[2] = d11; vec_1 = vec + 6;
1273  }
1274  }
1275 
1276  vec_1[0] = c1;
1277  vec_1[1] = c2;
1278  vec_1[2] = c3;
1279  d22 = g*(v2*c - v3*s);
1280  d33 = g*(v2*s + v3*c);
1281  vec_2[0] = - v1*d22; vec_3[0] = - v1*d33;
1282  vec_2[1] = c - v2*d22; vec_3[1] = s - v2*d33;
1283  vec_2[2] = -s - v3*d22; vec_3[2] = c - v3*d33;
1284  switch (k)
1285  {
1286  case 2:
1287  internal::Swap(vec_2[0], vec_2[1]);
1288  internal::Swap(vec_3[0], vec_3[1]);
1289  break;
1290 
1291  case 3:
1292  internal::Swap(vec_2[0], vec_2[2]);
1293  internal::Swap(vec_3[0], vec_3[2]);
1294  }
1295  }
1296 
1297 done_3d:
1298  lambda[0] *= mult;
1299  lambda[1] *= mult;
1300  lambda[2] *= mult;
1301 }
1302 
1303 /// Return the i'th singular value of the matrix of size 2 with given @a data.
1304 template<> MFEM_HOST_DEVICE inline
1305 double CalcSingularvalue<2>(const double *data, const int i)
1306 {
1307  double d0, d1, d2, d3;
1308  d0 = data[0];
1309  d1 = data[1];
1310  d2 = data[2];
1311  d3 = data[3];
1312  double mult;
1313 
1314  {
1315  double d_max = fabs(d0);
1316  if (d_max < fabs(d1)) { d_max = fabs(d1); }
1317  if (d_max < fabs(d2)) { d_max = fabs(d2); }
1318  if (d_max < fabs(d3)) { d_max = fabs(d3); }
1319  internal::GetScalingFactor(d_max, mult);
1320  }
1321 
1322  d0 /= mult;
1323  d1 /= mult;
1324  d2 /= mult;
1325  d3 /= mult;
1326 
1327  double t = 0.5*((d0+d2)*(d0-d2)+(d1-d3)*(d1+d3));
1328  double s = d0*d2 + d1*d3;
1329  s = sqrt(0.5*(d0*d0 + d1*d1 + d2*d2 + d3*d3) + sqrt(t*t + s*s));
1330 
1331  if (s == 0.0)
1332  {
1333  return 0.0;
1334  }
1335  t = fabs(d0*d3 - d1*d2) / s;
1336  if (t > s)
1337  {
1338  if (i == 0)
1339  {
1340  return t*mult;
1341  }
1342  return s*mult;
1343  }
1344  if (i == 0)
1345  {
1346  return s*mult;
1347  }
1348  return t*mult;
1349 }
1350 
1351 /// Return the i'th singular value of the matrix of size 3 with given @a data.
1352 template<> MFEM_HOST_DEVICE inline
1353 double CalcSingularvalue<3>(const double *data, const int i)
1354 {
1355  double d0, d1, d2, d3, d4, d5, d6, d7, d8;
1356  d0 = data[0]; d3 = data[3]; d6 = data[6];
1357  d1 = data[1]; d4 = data[4]; d7 = data[7];
1358  d2 = data[2]; d5 = data[5]; d8 = data[8];
1359  double mult;
1360  {
1361  double d_max = fabs(d0);
1362  if (d_max < fabs(d1)) { d_max = fabs(d1); }
1363  if (d_max < fabs(d2)) { d_max = fabs(d2); }
1364  if (d_max < fabs(d3)) { d_max = fabs(d3); }
1365  if (d_max < fabs(d4)) { d_max = fabs(d4); }
1366  if (d_max < fabs(d5)) { d_max = fabs(d5); }
1367  if (d_max < fabs(d6)) { d_max = fabs(d6); }
1368  if (d_max < fabs(d7)) { d_max = fabs(d7); }
1369  if (d_max < fabs(d8)) { d_max = fabs(d8); }
1370  internal::GetScalingFactor(d_max, mult);
1371  }
1372 
1373  d0 /= mult; d1 /= mult; d2 /= mult;
1374  d3 /= mult; d4 /= mult; d5 /= mult;
1375  d6 /= mult; d7 /= mult; d8 /= mult;
1376 
1377  double b11 = d0*d0 + d1*d1 + d2*d2;
1378  double b12 = d0*d3 + d1*d4 + d2*d5;
1379  double b13 = d0*d6 + d1*d7 + d2*d8;
1380  double b22 = d3*d3 + d4*d4 + d5*d5;
1381  double b23 = d3*d6 + d4*d7 + d5*d8;
1382  double b33 = d6*d6 + d7*d7 + d8*d8;
1383 
1384  // double a, b, c;
1385  // a = -(b11 + b22 + b33);
1386  // b = b11*(b22 + b33) + b22*b33 - b12*b12 - b13*b13 - b23*b23;
1387  // c = b11*(b23*b23 - b22*b33) + b12*(b12*b33 - 2*b13*b23) + b13*b13*b22;
1388 
1389  // double Q = (a * a - 3 * b) / 9;
1390  // double Q = (b12*b12 + b13*b13 + b23*b23 +
1391  // ((b11 - b22)*(b11 - b22) +
1392  // (b11 - b33)*(b11 - b33) +
1393  // (b22 - b33)*(b22 - b33))/6)/3;
1394  // Q = (3*(b12^2 + b13^2 + b23^2) +
1395  // ((b11 - b22)^2 + (b11 - b33)^2 + (b22 - b33)^2)/2)/9
1396  // or
1397  // Q = (1/6)*|B-tr(B)/3|_F^2
1398  // Q >= 0 and
1399  // Q = 0 <==> B = scalar * I
1400  // double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
1401  double aa = (b11 + b22 + b33)/3; // aa = tr(B)/3
1402  double c1, c2, c3;
1403  // c1 = b11 - aa; // ((b11 - b22) + (b11 - b33))/3
1404  // c2 = b22 - aa; // ((b22 - b11) + (b22 - b33))/3
1405  // c3 = b33 - aa; // ((b33 - b11) + (b33 - b22))/3
1406  {
1407  double b11_b22 = ((d0-d3)*(d0+d3)+(d1-d4)*(d1+d4)+(d2-d5)*(d2+d5));
1408  double b22_b33 = ((d3-d6)*(d3+d6)+(d4-d7)*(d4+d7)+(d5-d8)*(d5+d8));
1409  double b33_b11 = ((d6-d0)*(d6+d0)+(d7-d1)*(d7+d1)+(d8-d2)*(d8+d2));
1410  c1 = (b11_b22 - b33_b11)/3;
1411  c2 = (b22_b33 - b11_b22)/3;
1412  c3 = (b33_b11 - b22_b33)/3;
1413  }
1414  double Q, R;
1415  Q = (2*(b12*b12 + b13*b13 + b23*b23) + c1*c1 + c2*c2 + c3*c3)/6;
1416  R = (c1*(b23*b23 - c2*c3)+ b12*(b12*c3 - 2*b13*b23) +b13*b13*c2)/2;
1417  // R = (-1/2)*det(B-(tr(B)/3)*I)
1418  // Note: 54*(det(S))^2 <= |S|_F^6, when S^t=S and tr(S)=0, S is 3x3
1419  // Therefore: R^2 <= Q^3
1420 
1421  if (Q <= 0.) { ; }
1422 
1423  // else if (fabs(R) >= sqrtQ3)
1424  // {
1425  // double det = (d[0] * (d[4] * d[8] - d[5] * d[7]) +
1426  // d[3] * (d[2] * d[7] - d[1] * d[8]) +
1427  // d[6] * (d[1] * d[5] - d[2] * d[4]));
1428  //
1429  // if (R > 0.)
1430  // {
1431  // if (i == 2)
1432  // // aa -= 2*sqrtQ;
1433  // return fabs(det)/(aa + sqrtQ);
1434  // else
1435  // aa += sqrtQ;
1436  // }
1437  // else
1438  // {
1439  // if (i != 0)
1440  // aa -= sqrtQ;
1441  // // aa = fabs(det)/sqrt(aa + 2*sqrtQ);
1442  // else
1443  // aa += 2*sqrtQ;
1444  // }
1445  // }
1446 
1447  else
1448  {
1449  double sqrtQ = sqrt(Q);
1450  double sqrtQ3 = Q*sqrtQ;
1451  // double sqrtQ3 = sqrtQ*sqrtQ*sqrtQ;
1452  // double sqrtQ3 = pow(Q, 1.5);
1453  double r;
1454 
1455  if (fabs(R) >= sqrtQ3)
1456  {
1457  if (R < 0.)
1458  {
1459  // R = -1.;
1460  r = 2*sqrtQ;
1461  }
1462  else
1463  {
1464  // R = 1.;
1465  r = -2*sqrtQ;
1466  }
1467  }
1468  else
1469  {
1470  R = R/sqrtQ3;
1471 
1472  // if (fabs(R) <= 0.95)
1473  if (fabs(R) <= 0.9)
1474  {
1475  if (i == 2)
1476  {
1477  aa -= 2*sqrtQ*cos(acos(R)/3); // min
1478  }
1479  else if (i == 0)
1480  {
1481  aa -= 2*sqrtQ*cos((acos(R) + 2.0*M_PI)/3); // max
1482  }
1483  else
1484  {
1485  aa -= 2*sqrtQ*cos((acos(R) - 2.0*M_PI)/3); // mid
1486  }
1487  goto have_aa;
1488  }
1489 
1490  if (R < 0.)
1491  {
1492  r = -2*sqrtQ*cos((acos(R) + 2.0*M_PI)/3); // max
1493  if (i == 0)
1494  {
1495  aa += r;
1496  goto have_aa;
1497  }
1498  }
1499  else
1500  {
1501  r = -2*sqrtQ*cos(acos(R)/3); // min
1502  if (i == 2)
1503  {
1504  aa += r;
1505  goto have_aa;
1506  }
1507  }
1508  }
1509 
1510  // (tr(B)/3 + r) is the root which is separated from the other
1511  // two roots which are close to each other when |R| is close to 1
1512 
1513  c1 -= r;
1514  c2 -= r;
1515  c3 -= r;
1516  // aa += r;
1517 
1518  // Type of Householder reflections: z --> mu ek, where k is the index
1519  // of the entry in z with:
1520  // mode == 0: smallest absolute value --> angle closest to pi/2
1521  // (eliminate large entries)
1522  // mode == 1: largest absolute value --> angle farthest from pi/2
1523  // (eliminate small entries)
1524  const int mode = 1;
1525 
1526  // Find a unit vector z = (z1,z2,z3) in the "near"-kernel of
1527  // | c1 b12 b13 |
1528  // | b12 c2 b23 | = B - aa*I
1529  // | b13 b23 c3 |
1530  // This vector is also an eigenvector for B corresponding to aa
1531  // The vector z overwrites (c1,c2,c3).
1532  switch (internal::KernelVector3S(mode, b12, b13, b23, c1, c2, c3))
1533  {
1534  case 3:
1535  aa += r;
1536  goto have_aa;
1537  case 2:
1538  // ok, continue with the returned vector orthogonal to the kernel
1539  case 1:
1540  // ok, continue with the returned vector in the "near"-kernel
1541  ;
1542  }
1543 
1544  // Using the eigenvector c = (c1,c2,c3) to transform B into
1545  // | b11 0 0 |
1546  // B <-- Q P B P Q = | 0 b22 b23 |
1547  // | 0 b23 b33 |
1548  double v1, v2, v3, g;
1549  internal::Reduce3S(mode, b11, b22, b33, b12, b13, b23,
1550  c1, c2, c3, v1, v2, v3, g);
1551  // Q = I - g v v^t
1552  // P - permutation matrix switching rows and columns 1 and k
1553 
1554  // find the eigenvalues of
1555  // | b22 b23 |
1556  // | b23 b33 |
1557  internal::Eigenvalues2S(b23, b22, b33);
1558 
1559  if (i == 2)
1560  {
1561  aa = fmin(fmin(b11, b22), b33);
1562  }
1563  else if (i == 1)
1564  {
1565  if (b11 <= b22)
1566  {
1567  aa = (b22 <= b33) ? b22 : fmax(b11, b33);
1568  }
1569  else
1570  {
1571  aa = (b11 <= b33) ? b11 : fmax(b33, b22);
1572  }
1573  }
1574  else
1575  {
1576  aa = fmax(fmax(b11, b22), b33);
1577  }
1578  }
1579 
1580 have_aa:
1581 
1582  return sqrt(fabs(aa))*mult; // take abs before we sort?
1583 }
1584 
1585 
1586 /// Assuming L.U = P.A for a factored matrix (m x m),
1587 // compute x <- A x
1588 //
1589 // @param [in] data LU factorization of A
1590 // @param [in] m square matrix height
1591 // @param [in] ipiv array storing pivot information
1592 // @param [in, out] x vector storing right-hand side and then solution
1593 MFEM_HOST_DEVICE
1594 inline void LUSolve(const double *data, const int m, const int *ipiv,
1595  double *x)
1596 {
1597  // X <- P X
1598  for (int i = 0; i < m; i++)
1599  {
1600  internal::Swap<double>(x[i], x[ipiv[i]]);
1601  }
1602 
1603  // X <- L^{-1} X
1604  for (int j = 0; j < m; j++)
1605  {
1606  const double x_j = x[j];
1607  for (int i = j + 1; i < m; i++)
1608  {
1609  x[i] -= data[i + j * m] * x_j;
1610  }
1611  }
1612 
1613  // X <- U^{-1} X
1614  for (int j = m - 1; j >= 0; j--)
1615  {
1616  const double x_j = (x[j] /= data[j + j * m]);
1617  for (int i = 0; i < j; i++)
1618  {
1619  x[i] -= data[i + j * m] * x_j;
1620  }
1621  }
1622 }
1623 
1624 } // namespace kernels
1625 
1626 } // namespace mfem
1627 
1628 #endif // MFEM_LINALG_KERNELS_HPP
MFEM_HOST_DEVICE T Det(const T *data)
Compute the determinant of a square matrix of size dim with given data.
Definition: kernels.hpp:237
MFEM_HOST_DEVICE double DistanceSquared(const double *x, const double *y)
Compute the square of the Euclidean distance to another vector.
Definition: kernels.hpp:40
MFEM_HOST_DEVICE void FNorm(double &scale_factor, double &scaled_fnorm2, const T *data)
Definition: kernels.hpp:79
MFEM_HOST_DEVICE void MultABt(const int Aheight, const int Awidth, const int Bheight, const TA *Adata, const TB *Bdata, TC *ABtdata)
Multiply a matrix of size Aheight x Awidth and data Adata with the transpose of a matrix of size Bhei...
Definition: kernels.hpp:363
double epsilon
Definition: ex25.cpp:140
MFEM_HOST_DEVICE void Diag(const double c, double *data)
Creates n x n diagonal matrix with diagonal elements c.
Definition: kernels.hpp:49
MFEM_HOST_DEVICE void Add(const int height, const int width, const TALPHA alpha, const TA *Adata, const TB *Bdata, TC *Cdata)
Compute C = A + alpha*B, where the matrices A, B and C are of size height x width with data Adata...
Definition: kernels.hpp:266
FDualNumber< tbase > acos(const FDualNumber< tbase > &f)
acos([dual number])
Definition: fdual.hpp:423
MFEM_HOST_DEVICE void MultAtB(const int Aheight, const int Awidth, const int Bwidth, const TA *Adata, const TB *Bdata, TC *AtBdata)
Multiply the transpose of a matrix of size Aheight x Awidth and data Adata with a matrix of size Ahei...
Definition: kernels.hpp:390
MFEM_HOST_DEVICE double CalcSingularvalue< 3 >(const double *data, const int i)
Return the i&#39;th singular value of the matrix of size 3 with given data.
Definition: kernels.hpp:1353
MFEM_HOST_DEVICE void CalcEigenvalues(const double *data, double *lambda, double *vec)
MFEM_HOST_DEVICE void CalcEigenvalues< 2 >(const double *data, double *lambda, double *vec)
Definition: kernels.hpp:1072
double b
Definition: lissajous.cpp:42
MFEM_HOST_DEVICE void LUSolve(const double *data, const int m, const int *ipiv, double *x)
Assuming L.U = P.A for a factored matrix (m x m),.
Definition: kernels.hpp:1594
MFEM_HOST_DEVICE double Norml2(const int size, const T *data)
Returns the l2 norm of the Vector with given size and data.
Definition: kernels.hpp:133
MFEM_HOST_DEVICE void CalcAdjugate(const T *data, T *adj_data)
Return the adjugate of a matrix.
Definition: kernels.hpp:256
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
MFEM_HOST_DEVICE void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
Matrix vector multiplication: y = A x, where the matrix A is of size height x width with given data...
Definition: kernels.hpp:163
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:630
MFEM_HOST_DEVICE void MultTranspose(const int height, const int width, const TA *data, const TX *x, TY *y)
Matrix transpose vector multiplication: y = At x, where the matrix A is of size height x width with g...
Definition: kernels.hpp:196
MFEM_HOST_DEVICE double CalcSingularvalue(const double *data, const int i)
Return the i&#39;th singular value of the matrix of size dim with given data.
MFEM_HOST_DEVICE double CalcSingularvalue< 2 >(const double *data, const int i)
Return the i&#39;th singular value of the matrix of size 2 with given data.
Definition: kernels.hpp:1305
double a
Definition: lissajous.cpp:41
MFEM_HOST_DEVICE void CalcEigenvalues< 3 >(const double *data, double *lambda, double *vec)
Definition: kernels.hpp:1103
double mu
Definition: ex25.cpp:139
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
MFEM_HOST_DEVICE void CalcInverse(const T *data, T *inv_data)
Return the inverse of a matrix with given size and data into the matrix with data inv_data...
Definition: kernels.hpp:246
int dim
Definition: ex24.cpp:53
MFEM_HOST_DEVICE void Subtract(const double a, const double *x, const double *y, double *z)
Vector subtraction operation: z = a * (x - y)
Definition: kernels.hpp:58
RefCoord t[3]
MFEM_HOST_DEVICE void Set(const int height, const int width, const double alpha, const TA *Adata, TB *Bdata)
Compute B = alpha*A, where the matrices A and B are of size height x width with data Adata and Bdata...
Definition: kernels.hpp:326
const double alpha
Definition: ex15.cpp:369
RefCoord s[3]
MFEM_HOST_DEVICE void AddMultVWt(const double *v, const double *w, double *VWt)
Dense matrix operation: VWt += v w^t.
Definition: kernels.hpp:67
MFEM_HOST_DEVICE void Symmetrize(const int size, T *data)
Symmetrize a square matrix with given size and data: A -&gt; (A+A^T)/2.
Definition: kernels.hpp:223
MFEM_HOST_DEVICE double FNorm2(const T *data)
Compute the square of the Frobenius norm of the matrix.
Definition: kernels.hpp:123