MFEM  v4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
densemat.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 
13 // Implementation of data types dense matrix, inverse dense matrix
14 
15 
16 #include "vector.hpp"
17 #include "matrix.hpp"
18 #include "densemat.hpp"
19 #include "../general/table.hpp"
20 #include "../general/globals.hpp"
21 
22 #include <iostream>
23 #include <iomanip>
24 #include <limits>
25 #include <algorithm>
26 #include <cstdlib>
27 #if defined(_MSC_VER) && (_MSC_VER < 1800)
28 #include <float.h>
29 #define copysign _copysign
30 #endif
31 
32 
33 #ifdef MFEM_USE_LAPACK
34 extern "C" void
35 dgemm_(char *, char *, int *, int *, int *, double *, double *,
36  int *, double *, int *, double *, double *, int *);
37 extern "C" void
38 dgetrf_(int *, int *, double *, int *, int *, int *);
39 extern "C" void
40 dgetrs_(char *, int *, int *, double *, int *, int *, double *, int *, int *);
41 extern "C" void
42 dgetri_(int *N, double *A, int *LDA, int *IPIV, double *WORK,
43  int *LWORK, int *INFO);
44 extern "C" void
45 dsyevr_(char *JOBZ, char *RANGE, char *UPLO, int *N, double *A, int *LDA,
46  double *VL, double *VU, int *IL, int *IU, double *ABSTOL, int *M,
47  double *W, double *Z, int *LDZ, int *ISUPPZ, double *WORK, int *LWORK,
48  int *IWORK, int *LIWORK, int *INFO);
49 extern "C" void
50 dsyev_(char *JOBZ, char *UPLO, int *N, double *A, int *LDA, double *W,
51  double *WORK, int *LWORK, int *INFO);
52 extern "C" void
53 dsygv_ (int *ITYPE, char *JOBZ, char *UPLO, int * N, double *A, int *LDA,
54  double *B, int *LDB, double *W, double *WORK, int *LWORK, int *INFO);
55 extern "C" void
56 dgesvd_(char *JOBU, char *JOBVT, int *M, int *N, double *A, int *LDA,
57  double *S, double *U, int *LDU, double *VT, int *LDVT, double *WORK,
58  int *LWORK, int *INFO);
59 #endif
60 
61 
62 namespace mfem
63 {
64 
65 using namespace std;
66 
68 {
69  data = NULL;
70  capacity = 0;
71 }
72 
73 DenseMatrix::DenseMatrix(const DenseMatrix &m) : Matrix(m.height, m.width)
74 {
75  const int hw = height * width;
76  if (hw > 0)
77  {
78  MFEM_ASSERT(m.data, "invalid source matrix");
79  data = new double[hw];
80  capacity = hw;
81  std::memcpy(data, m.data, sizeof(double)*hw);
82  }
83  else
84  {
85  data = NULL;
86  capacity = 0;
87  }
88 }
89 
91 {
92  MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
93  capacity = s*s;
94  if (capacity > 0)
95  {
96  data = new double[capacity](); // init with zeroes
97  }
98  else
99  {
100  data = NULL;
101  }
102 }
103 
104 DenseMatrix::DenseMatrix(int m, int n) : Matrix(m, n)
105 {
106  MFEM_ASSERT(m >= 0 && n >= 0,
107  "invalid DenseMatrix size: " << m << " x " << n);
108  capacity = m*n;
109  if (capacity > 0)
110  {
111  data = new double[capacity](); // init with zeroes
112  }
113  else
114  {
115  data = NULL;
116  }
117 }
118 
120  : Matrix(mat.width, mat.height)
121 {
122  capacity = height*width;
123  if (capacity > 0)
124  {
125  data = new double[capacity];
126 
127  for (int i = 0; i < height; i++)
128  {
129  for (int j = 0; j < width; j++)
130  {
131  (*this)(i,j) = mat(j,i);
132  }
133  }
134  }
135  else
136  {
137  data = NULL;
138  }
139 }
140 
141 void DenseMatrix::SetSize(int h, int w)
142 {
143  MFEM_ASSERT(h >= 0 && w >= 0,
144  "invalid DenseMatrix size: " << h << " x " << w);
145  if (Height() == h && Width() == w)
146  {
147  return;
148  }
149  height = h;
150  width = w;
151  const int hw = h*w;
152  if (hw > std::abs(capacity))
153  {
154  if (capacity > 0)
155  {
156  delete [] data;
157  }
158  capacity = hw;
159  data = new double[hw](); // init with zeroes
160  }
161 }
162 
163 double &DenseMatrix::Elem(int i, int j)
164 {
165  return (*this)(i,j);
166 }
167 
168 const double &DenseMatrix::Elem(int i, int j) const
169 {
170  return (*this)(i,j);
171 }
172 
173 void DenseMatrix::Mult(const double *x, double *y) const
174 {
175  if (width == 0)
176  {
177  for (int row = 0; row < height; row++)
178  {
179  y[row] = 0.0;
180  }
181  return;
182  }
183  double *d_col = data;
184  double x_col = x[0];
185  for (int row = 0; row < height; row++)
186  {
187  y[row] = x_col*d_col[row];
188  }
189  d_col += height;
190  for (int col = 1; col < width; col++)
191  {
192  x_col = x[col];
193  for (int row = 0; row < height; row++)
194  {
195  y[row] += x_col*d_col[row];
196  }
197  d_col += height;
198  }
199 }
200 
201 void DenseMatrix::Mult(const Vector &x, Vector &y) const
202 {
203  MFEM_ASSERT(height == y.Size() && width == x.Size(),
204  "incompatible dimensions");
205 
206  Mult((const double *)x, (double *)y);
207 }
208 
209 double DenseMatrix::operator *(const DenseMatrix &m) const
210 {
211  MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
212  "incompatible dimensions");
213 
214  const int hw = height * width;
215  double a = 0.0;
216  for (int i = 0; i < hw; i++)
217  {
218  a += data[i] * m.data[i];
219  }
220 
221  return a;
222 }
223 
224 void DenseMatrix::MultTranspose(const double *x, double *y) const
225 {
226  double *d_col = data;
227  for (int col = 0; col < width; col++)
228  {
229  double y_col = 0.0;
230  for (int row = 0; row < height; row++)
231  {
232  y_col += x[row]*d_col[row];
233  }
234  y[col] = y_col;
235  d_col += height;
236  }
237 }
238 
239 void DenseMatrix::MultTranspose(const Vector &x, Vector &y) const
240 {
241  MFEM_ASSERT(height == x.Size() && width == y.Size(),
242  "incompatible dimensions");
243 
244  MultTranspose((const double *)x, (double *)y);
245 }
246 
247 void DenseMatrix::AddMult(const Vector &x, Vector &y) const
248 {
249  MFEM_ASSERT(height == y.Size() && width == x.Size(),
250  "incompatible dimensions");
251 
252  const double *xp = x, *d_col = data;
253  double *yp = y;
254  for (int col = 0; col < width; col++)
255  {
256  double x_col = xp[col];
257  for (int row = 0; row < height; row++)
258  {
259  yp[row] += x_col*d_col[row];
260  }
261  d_col += height;
262  }
263 }
264 
266 {
267  MFEM_ASSERT(height == x.Size() && width == y.Size(),
268  "incompatible dimensions");
269 
270  const double *d_col = data;
271  for (int col = 0; col < width; col++)
272  {
273  double y_col = 0.0;
274  for (int row = 0; row < height; row++)
275  {
276  y_col += x[row]*d_col[row];
277  }
278  y[col] += y_col;
279  d_col += height;
280  }
281 }
282 
283 void DenseMatrix::AddMult_a(double a, const Vector &x, Vector &y) const
284 {
285  MFEM_ASSERT(height == y.Size() && width == x.Size(),
286  "incompatible dimensions");
287 
288  const double *xp = x, *d_col = data;
289  double *yp = y;
290  for (int col = 0; col < width; col++)
291  {
292  const double x_col = a*xp[col];
293  for (int row = 0; row < height; row++)
294  {
295  yp[row] += x_col*d_col[row];
296  }
297  d_col += height;
298  }
299 }
300 
301 void DenseMatrix::AddMultTranspose_a(double a, const Vector &x,
302  Vector &y) const
303 {
304  MFEM_ASSERT(height == x.Size() && width == y.Size(),
305  "incompatible dimensions");
306 
307  const double *d_col = data;
308  for (int col = 0; col < width; col++)
309  {
310  double y_col = 0.0;
311  for (int row = 0; row < height; row++)
312  {
313  y_col += x[row]*d_col[row];
314  }
315  y[col] += a * y_col;
316  d_col += height;
317  }
318 }
319 
320 double DenseMatrix::InnerProduct(const double *x, const double *y) const
321 {
322  double prod = 0.0;
323 
324  for (int i = 0; i < height; i++)
325  {
326  double Axi = 0.0;
327  for (int j = 0; j < width; j++)
328  {
329  Axi += (*this)(i,j) * x[j];
330  }
331  prod += y[i] * Axi;
332  }
333 
334  return prod;
335 }
336 
337 // LeftScaling this = diag(s) * this
339 {
340  double * it_data = data;
341  for (int j = 0; j < width; ++j)
342  {
343  for (int i = 0; i < height; ++i)
344  {
345  *(it_data++) *= s(i);
346  }
347  }
348 }
349 
350 // InvLeftScaling this = diag(1./s) * this
352 {
353  double * it_data = data;
354  for (int j = 0; j < width; ++j)
355  {
356  for (int i = 0; i < height; ++i)
357  {
358  *(it_data++) /= s(i);
359  }
360  }
361 }
362 
363 // RightScaling: this = this * diag(s);
365 {
366  double sj;
367  double * it_data = data;
368  for (int j = 0; j < width; ++j)
369  {
370  sj = s(j);
371  for (int i = 0; i < height; ++i)
372  {
373  *(it_data++) *= sj;
374  }
375  }
376 }
377 
378 // InvRightScaling: this = this * diag(1./s);
380 {
381  double * it_data = data;
382  for (int j = 0; j < width; ++j)
383  {
384  const double sj = 1./s(j);
385  for (int i = 0; i < height; ++i)
386  {
387  *(it_data++) *= sj;
388  }
389  }
390 }
391 
392 // SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
394 {
395  if (height != width || s.Size() != height)
396  {
397  mfem_error("DenseMatrix::SymmetricScaling");
398  }
399 
400  double * ss = new double[width];
401  double * it_s = s.GetData();
402  double * it_ss = ss;
403  for ( double * end_s = it_s + width; it_s != end_s; ++it_s)
404  {
405  *(it_ss++) = sqrt(*it_s);
406  }
407 
408  double * it_data = data;
409  for (int j = 0; j < width; ++j)
410  {
411  for (int i = 0; i < height; ++i)
412  {
413  *(it_data++) *= ss[i]*ss[j];
414  }
415  }
416 
417  delete[] ss;
418 }
419 
420 // InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
422 {
423  if (height != width || s.Size() != width)
424  {
425  mfem_error("DenseMatrix::SymmetricScaling");
426  }
427 
428  double * ss = new double[width];
429  double * it_s = s.GetData();
430  double * it_ss = ss;
431  for (double * end_s = it_s + width; it_s != end_s; ++it_s)
432  {
433  *(it_ss++) = 1./sqrt(*it_s);
434  }
435 
436  double * it_data = data;
437  for (int j = 0; j < width; ++j)
438  {
439  for (int i = 0; i < height; ++i)
440  {
441  *(it_data++) *= ss[i]*ss[j];
442  }
443  }
444 
445  delete[] ss;
446 }
447 
448 double DenseMatrix::Trace() const
449 {
450 #ifdef MFEM_DEBUG
451  if (Width() != Height())
452  {
453  mfem_error("DenseMatrix::Trace() : not a square matrix!");
454  }
455 #endif
456 
457  double t = 0.0;
458 
459  for (int i = 0; i < width; i++)
460  {
461  t += (*this)(i, i);
462  }
463 
464  return t;
465 }
466 
468 {
469  return new DenseMatrixInverse(*this);
470 }
471 
472 double DenseMatrix::Det() const
473 {
474  MFEM_ASSERT(Height() == Width() && Height() > 0,
475  "The matrix must be square and "
476  << "sized larger than zero to compute the determinant."
477  << " Height() = " << Height()
478  << ", Width() = " << Width());
479 
480  switch (Height())
481  {
482  case 1:
483  return data[0];
484 
485  case 2:
486  return data[0] * data[3] - data[1] * data[2];
487 
488  case 3:
489  {
490  const double *d = data;
491  return
492  d[0] * (d[4] * d[8] - d[5] * d[7]) +
493  d[3] * (d[2] * d[7] - d[1] * d[8]) +
494  d[6] * (d[1] * d[5] - d[2] * d[4]);
495  }
496  case 4:
497  {
498  const double *d = data;
499  return
500  d[ 0] * (d[ 5] * (d[10] * d[15] - d[11] * d[14]) -
501  d[ 9] * (d[ 6] * d[15] - d[ 7] * d[14]) +
502  d[13] * (d[ 6] * d[11] - d[ 7] * d[10])
503  ) -
504  d[ 4] * (d[ 1] * (d[10] * d[15] - d[11] * d[14]) -
505  d[ 9] * (d[ 2] * d[15] - d[ 3] * d[14]) +
506  d[13] * (d[ 2] * d[11] - d[ 3] * d[10])
507  ) +
508  d[ 8] * (d[ 1] * (d[ 6] * d[15] - d[ 7] * d[14]) -
509  d[ 5] * (d[ 2] * d[15] - d[ 3] * d[14]) +
510  d[13] * (d[ 2] * d[ 7] - d[ 3] * d[ 6])
511  ) -
512  d[12] * (d[ 1] * (d[ 6] * d[11] - d[ 7] * d[10]) -
513  d[ 5] * (d[ 2] * d[11] - d[ 3] * d[10]) +
514  d[ 9] * (d[ 2] * d[ 7] - d[ 3] * d[ 6])
515  );
516  }
517  default:
518  {
519  // In the general case we compute the determinant from the LU
520  // decomposition.
521  DenseMatrixInverse lu_factors(*this);
522 
523  return lu_factors.Det();
524  }
525  }
526  // not reachable
527 }
528 
529 double DenseMatrix::Weight() const
530 {
531  if (Height() == Width())
532  {
533  // return fabs(Det());
534  return Det();
535  }
536  else if ((Height() == 2) && (Width() == 1))
537  {
538  return sqrt(data[0] * data[0] + data[1] * data[1]);
539  }
540  else if ((Height() == 3) && (Width() == 1))
541  {
542  return sqrt(data[0] * data[0] + data[1] * data[1] + data[2] * data[2]);
543  }
544  else if ((Height() == 3) && (Width() == 2))
545  {
546  const double *d = data;
547  double E = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
548  double G = d[3] * d[3] + d[4] * d[4] + d[5] * d[5];
549  double F = d[0] * d[3] + d[1] * d[4] + d[2] * d[5];
550  return sqrt(E * G - F * F);
551  }
552  mfem_error("DenseMatrix::Weight()");
553  return 0.0;
554 }
555 
556 void DenseMatrix::Set(double alpha, const double *A)
557 {
558  const int s = Width()*Height();
559  for (int i = 0; i < s; i++)
560  {
561  data[i] = alpha*A[i];
562  }
563 }
564 
565 void DenseMatrix::Add(const double c, const DenseMatrix &A)
566 {
567  for (int j = 0; j < Width(); j++)
568  {
569  for (int i = 0; i < Height(); i++)
570  {
571  (*this)(i,j) += c * A(i,j);
572  }
573  }
574 }
575 
577 {
578  const int s = Height()*Width();
579  for (int i = 0; i < s; i++)
580  {
581  data[i] = c;
582  }
583  return *this;
584 }
585 
587 {
588  const int s = Height()*Width();
589  for (int i = 0; i < s; i++)
590  {
591  data[i] = d[i];
592  }
593  return *this;
594 }
595 
597 {
598  SetSize(m.height, m.width);
599 
600  const int hw = height * width;
601  for (int i = 0; i < hw; i++)
602  {
603  data[i] = m.data[i];
604  }
605 
606  return *this;
607 }
608 
610 {
611  const int hw = Height()*Width();
612  for (int i = 0; i < hw; i++)
613  {
614  data[i] += m[i];
615  }
616  return *this;
617 }
618 
620 {
621  MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
622  "incompatible matrix sizes.");
623  return *this += m.GetData();
624 }
625 
627 {
628  for (int j = 0; j < width; j++)
629  {
630  for (int i = 0; i < height; i++)
631  {
632  (*this)(i, j) -= m(i, j);
633  }
634  }
635 
636  return *this;
637 }
638 
640 {
641  int s = Height()*Width();
642  for (int i = 0; i < s; i++)
643  {
644  data[i] *= c;
645  }
646  return *this;
647 }
648 
650 {
651  const int hw = Height() * Width();
652  for (int i = 0; i < hw; i++)
653  {
654  data[i] = -data[i];
655  }
656 }
657 
659 {
660 #ifdef MFEM_DEBUG
661  if (Height() <= 0 || Height() != Width())
662  {
663  mfem_error("DenseMatrix::Invert()");
664  }
665 #endif
666 
667 #ifdef MFEM_USE_LAPACK
668  int *ipiv = new int[width];
669  int lwork = -1;
670  double qwork, *work;
671  int info;
672 
673  dgetrf_(&width, &width, data, &width, ipiv, &info);
674 
675  if (info)
676  {
677  mfem_error("DenseMatrix::Invert() : Error in DGETRF");
678  }
679 
680  dgetri_(&width, data, &width, ipiv, &qwork, &lwork, &info);
681 
682  lwork = (int) qwork;
683  work = new double[lwork];
684 
685  dgetri_(&width, data, &width, ipiv, work, &lwork, &info);
686 
687  if (info)
688  {
689  mfem_error("DenseMatrix::Invert() : Error in DGETRI");
690  }
691 
692  delete [] work;
693  delete [] ipiv;
694 #else
695  int c, i, j, n = Width();
696  double a, b;
697  Array<int> piv(n);
698 
699  for (c = 0; c < n; c++)
700  {
701  a = fabs((*this)(c, c));
702  i = c;
703  for (j = c + 1; j < n; j++)
704  {
705  b = fabs((*this)(j, c));
706  if (a < b)
707  {
708  a = b;
709  i = j;
710  }
711  }
712  if (a == 0.0)
713  {
714  mfem_error("DenseMatrix::Invert() : singular matrix");
715  }
716  piv[c] = i;
717  for (j = 0; j < n; j++)
718  {
719  Swap<double>((*this)(c, j), (*this)(i, j));
720  }
721 
722  a = (*this)(c, c) = 1.0 / (*this)(c, c);
723  for (j = 0; j < c; j++)
724  {
725  (*this)(c, j) *= a;
726  }
727  for (j++; j < n; j++)
728  {
729  (*this)(c, j) *= a;
730  }
731  for (i = 0; i < c; i++)
732  {
733  (*this)(i, c) = a * (b = -(*this)(i, c));
734  for (j = 0; j < c; j++)
735  {
736  (*this)(i, j) += b * (*this)(c, j);
737  }
738  for (j++; j < n; j++)
739  {
740  (*this)(i, j) += b * (*this)(c, j);
741  }
742  }
743  for (i++; i < n; i++)
744  {
745  (*this)(i, c) = a * (b = -(*this)(i, c));
746  for (j = 0; j < c; j++)
747  {
748  (*this)(i, j) += b * (*this)(c, j);
749  }
750  for (j++; j < n; j++)
751  {
752  (*this)(i, j) += b * (*this)(c, j);
753  }
754  }
755  }
756 
757  for (c = n - 1; c >= 0; c--)
758  {
759  j = piv[c];
760  for (i = 0; i < n; i++)
761  {
762  Swap<double>((*this)(i, c), (*this)(i, j));
763  }
764  }
765 #endif
766 }
767 
769 {
770  // Square root inverse using Denman--Beavers
771 #ifdef MFEM_DEBUG
772  if (Height() <= 0 || Height() != Width())
773  {
774  mfem_error("DenseMatrix::SquareRootInverse() matrix not square.");
775  }
776 #endif
777 
778  DenseMatrix tmp1(Height());
779  DenseMatrix tmp2(Height());
780  DenseMatrix tmp3(Height());
781 
782  tmp1 = (*this);
783  (*this) = 0.0;
784  for (int v = 0; v < Height() ; v++) { (*this)(v,v) = 1.0; }
785 
786  for (int j = 0; j < 10; j++)
787  {
788  for (int i = 0; i < 10; i++)
789  {
790  tmp2 = tmp1;
791  tmp3 = (*this);
792 
793  tmp2.Invert();
794  tmp3.Invert();
795 
796  tmp1 += tmp3;
797  (*this) += tmp2;
798 
799  tmp1 *= 0.5;
800  (*this) *= 0.5;
801  }
802  mfem::Mult((*this), tmp1, tmp2);
803  for (int v = 0; v < Height() ; v++) { tmp2(v,v) -= 1.0; }
804  if (tmp2.FNorm() < 1e-10) { break; }
805  }
806 
807  if (tmp2.FNorm() > 1e-10)
808  {
809  mfem_error("DenseMatrix::SquareRootInverse not converged");
810  }
811 }
812 
813 void DenseMatrix::Norm2(double *v) const
814 {
815  for (int j = 0; j < Width(); j++)
816  {
817  v[j] = 0.0;
818  for (int i = 0; i < Height(); i++)
819  {
820  v[j] += (*this)(i,j)*(*this)(i,j);
821  }
822  v[j] = sqrt(v[j]);
823  }
824 }
825 
827 {
828  int hw = Height()*Width();
829  const double *d = data;
830  double norm = 0.0, abs_entry;
831 
832  for (int i = 0; i < hw; i++)
833  {
834  abs_entry = fabs(d[i]);
835  if (norm < abs_entry)
836  {
837  norm = abs_entry;
838  }
839  }
840 
841  return norm;
842 }
843 
844 void DenseMatrix::FNorm(double &scale_factor, double &scaled_fnorm2) const
845 {
846  int i, hw = Height() * Width();
847  double max_norm = 0.0, entry, fnorm2;
848 
849  for (i = 0; i < hw; i++)
850  {
851  entry = fabs(data[i]);
852  if (entry > max_norm)
853  {
854  max_norm = entry;
855  }
856  }
857 
858  if (max_norm == 0.0)
859  {
860  scale_factor = scaled_fnorm2 = 0.0;
861  return;
862  }
863 
864  fnorm2 = 0.0;
865  for (i = 0; i < hw; i++)
866  {
867  entry = data[i] / max_norm;
868  fnorm2 += entry * entry;
869  }
870 
871  scale_factor = max_norm;
872  scaled_fnorm2 = fnorm2;
873 }
874 
876 {
877 #ifdef MFEM_USE_LAPACK
878  ev.SetSize(a.Width());
879 
880  char JOBZ = 'N';
881  char RANGE = 'A';
882  char UPLO = 'U';
883  int N = a.Width();
884  double *A = new double[N*N];
885  int LDA = N;
886  double VL = 0.0;
887  double VU = 1.0;
888  int IL = 0;
889  int IU = 1;
890  double ABSTOL = 0.0;
891  int M;
892  double *W = ev.GetData();
893  double *Z = NULL;
894  int LDZ = 1;
895  int *ISUPPZ = new int[2*N];
896  int LWORK = -1; // query optimal (double) workspace size
897  double QWORK;
898  double *WORK = NULL;
899  int LIWORK = -1; // query optimal (int) workspace size
900  int QIWORK;
901  int *IWORK = NULL;
902  int INFO;
903 
904  if (evect) // Compute eigenvectors too
905  {
906  evect->SetSize(N);
907 
908  JOBZ = 'V';
909  Z = evect->Data();
910  LDZ = N;
911  }
912 
913  int hw = a.Height() * a.Width();
914  double *data = a.Data();
915 
916  for (int i = 0; i < hw; i++)
917  {
918  A[i] = data[i];
919  }
920 
921  dsyevr_( &JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU,
922  &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, &QWORK, &LWORK,
923  &QIWORK, &LIWORK, &INFO );
924 
925  LWORK = (int) QWORK;
926  LIWORK = QIWORK;
927 
928  WORK = new double[LWORK];
929  IWORK = new int[LIWORK];
930 
931  dsyevr_( &JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU,
932  &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, WORK, &LWORK,
933  IWORK, &LIWORK, &INFO );
934 
935  if (INFO != 0)
936  {
937  mfem::err << "dsyevr_Eigensystem(...): DSYEVR error code: "
938  << INFO << endl;
939  mfem_error();
940  }
941 
942 #ifdef MFEM_DEBUG
943  if (M < N)
944  {
945  mfem::err << "dsyevr_Eigensystem(...):\n"
946  << " DSYEVR did not find all eigenvalues "
947  << M << "/" << N << endl;
948  mfem_error();
949  }
950  if (CheckFinite(W, N) > 0)
951  {
952  mfem_error("dsyevr_Eigensystem(...): inf/nan values in W");
953  }
954  if (CheckFinite(Z, N*N) > 0)
955  {
956  mfem_error("dsyevr_Eigensystem(...): inf/nan values in Z");
957  }
958  VU = 0.0;
959  for (IL = 0; IL < N; IL++)
960  for (IU = 0; IU <= IL; IU++)
961  {
962  VL = 0.0;
963  for (M = 0; M < N; M++)
964  {
965  VL += Z[M+IL*N] * Z[M+IU*N];
966  }
967  if (IU < IL)
968  {
969  VL = fabs(VL);
970  }
971  else
972  {
973  VL = fabs(VL-1.0);
974  }
975  if (VL > VU)
976  {
977  VU = VL;
978  }
979  if (VU > 0.5)
980  {
981  mfem::err << "dsyevr_Eigensystem(...):"
982  << " Z^t Z - I deviation = " << VU
983  << "\n W[max] = " << W[N-1] << ", W[min] = "
984  << W[0] << ", N = " << N << endl;
985  mfem_error();
986  }
987  }
988  if (VU > 1e-9)
989  {
990  mfem::err << "dsyevr_Eigensystem(...):"
991  << " Z^t Z - I deviation = " << VU
992  << "\n W[max] = " << W[N-1] << ", W[min] = "
993  << W[0] << ", N = " << N << endl;
994  }
995  if (VU > 1e-5)
996  {
997  mfem_error("dsyevr_Eigensystem(...): ERROR: ...");
998  }
999  VU = 0.0;
1000  for (IL = 0; IL < N; IL++)
1001  for (IU = 0; IU < N; IU++)
1002  {
1003  VL = 0.0;
1004  for (M = 0; M < N; M++)
1005  {
1006  VL += Z[IL+M*N] * W[M] * Z[IU+M*N];
1007  }
1008  VL = fabs(VL-data[IL+N*IU]);
1009  if (VL > VU)
1010  {
1011  VU = VL;
1012  }
1013  }
1014  if (VU > 1e-9)
1015  {
1016  mfem::err << "dsyevr_Eigensystem(...):"
1017  << " max matrix deviation = " << VU
1018  << "\n W[max] = " << W[N-1] << ", W[min] = "
1019  << W[0] << ", N = " << N << endl;
1020  }
1021  if (VU > 1e-5)
1022  {
1023  mfem_error("dsyevr_Eigensystem(...): ERROR: ...");
1024  }
1025 #endif
1026 
1027  delete [] IWORK;
1028  delete [] WORK;
1029  delete [] ISUPPZ;
1030  delete [] A;
1031 
1032 #endif
1033 }
1034 
1036 {
1037 #ifdef MFEM_USE_LAPACK
1038  int N = a.Width();
1039  char JOBZ = 'N';
1040  char UPLO = 'U';
1041  int LDA = N;
1042  int LWORK = -1; /* query optimal workspace size */
1043  int INFO;
1044 
1045  ev.SetSize(N);
1046 
1047  double *A = NULL;
1048  double *W = ev.GetData();
1049  double *WORK = NULL;
1050  double QWORK;
1051 
1052  if (evect)
1053  {
1054  JOBZ = 'V';
1055  evect->SetSize(N);
1056  A = evect->Data();
1057  }
1058  else
1059  {
1060  A = new double[N*N];
1061  }
1062 
1063  int hw = a.Height() * a.Width();
1064  double *data = a.Data();
1065  for (int i = 0; i < hw; i++)
1066  {
1067  A[i] = data[i];
1068  }
1069 
1070  dsyev_(&JOBZ, &UPLO, &N, A, &LDA, W, &QWORK, &LWORK, &INFO);
1071 
1072  LWORK = (int) QWORK;
1073  WORK = new double[LWORK];
1074 
1075  dsyev_(&JOBZ, &UPLO, &N, A, &LDA, W, WORK, &LWORK, &INFO);
1076 
1077  if (INFO != 0)
1078  {
1079  mfem::err << "dsyev_Eigensystem: DSYEV error code: " << INFO << endl;
1080  mfem_error();
1081  }
1082 
1083  delete [] WORK;
1084  if (evect == NULL) { delete [] A; }
1085 #endif
1086 }
1087 
1088 void DenseMatrix::Eigensystem(Vector &ev, DenseMatrix *evect)
1089 {
1090 #ifdef MFEM_USE_LAPACK
1091 
1092  // dsyevr_Eigensystem(*this, ev, evect);
1093 
1094  dsyev_Eigensystem(*this, ev, evect);
1095 
1096 #else
1097 
1098  mfem_error("DenseMatrix::Eigensystem");
1099 
1100 #endif
1101 }
1102 
1104  DenseMatrix *evect)
1105 {
1106 #ifdef MFEM_USE_LAPACK
1107  int N = a.Width();
1108  int ITYPE = 1;
1109  char JOBZ = 'N';
1110  char UPLO = 'U';
1111  int LDA = N;
1112  int LDB = N;
1113  int LWORK = -1; /* query optimal workspace size */
1114  int INFO;
1115 
1116  ev.SetSize(N);
1117 
1118  double *A = NULL;
1119  double *B = new double[N*N];
1120  double *W = ev.GetData();
1121  double *WORK = NULL;
1122  double QWORK;
1123 
1124  if (evect)
1125  {
1126  JOBZ = 'V';
1127  evect->SetSize(N);
1128  A = evect->Data();
1129  }
1130  else
1131  {
1132  A = new double[N*N];
1133  }
1134 
1135  int hw = a.Height() * a.Width();
1136  double *a_data = a.Data();
1137  double *b_data = b.Data();
1138  for (int i = 0; i < hw; i++)
1139  {
1140  A[i] = a_data[i];
1141  B[i] = b_data[i];
1142  }
1143 
1144  dsygv_(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, &QWORK, &LWORK, &INFO);
1145 
1146  LWORK = (int) QWORK;
1147  WORK = new double[LWORK];
1148 
1149  dsygv_(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, WORK, &LWORK, &INFO);
1150 
1151  if (INFO != 0)
1152  {
1153  mfem::err << "dsygv_Eigensystem: DSYGV error code: " << INFO << endl;
1154  mfem_error();
1155  }
1156 
1157  delete [] WORK;
1158  delete [] B;
1159  if (evect == NULL) { delete [] A; }
1160 #endif
1161 }
1162 
1163 void DenseMatrix::Eigensystem(DenseMatrix &b, Vector &ev,
1164  DenseMatrix *evect)
1165 {
1166 #ifdef MFEM_USE_LAPACK
1167 
1168  dsygv_Eigensystem(*this, b, ev, evect);
1169 
1170 #else
1171 
1172  mfem_error("DenseMatrix::Eigensystem for generalized eigenvalues");
1173 #endif
1174 }
1175 
1177 {
1178 #ifdef MFEM_USE_LAPACK
1179  DenseMatrix copy_of_this = *this;
1180  char jobu = 'N';
1181  char jobvt = 'N';
1182  int m = Height();
1183  int n = Width();
1184  double *a = copy_of_this.data;
1185  sv.SetSize(min(m, n));
1186  double *s = sv;
1187  double *u = NULL;
1188  double *vt = NULL;
1189  double *work = NULL;
1190  int lwork = -1;
1191  int info;
1192  double qwork;
1193 
1194  dgesvd_(&jobu, &jobvt, &m, &n, a, &m,
1195  s, u, &m, vt, &n, &qwork, &lwork, &info);
1196 
1197  lwork = (int) qwork;
1198  work = new double[lwork];
1199 
1200  dgesvd_(&jobu, &jobvt, &m, &n, a, &m,
1201  s, u, &m, vt, &n, work, &lwork, &info);
1202 
1203  delete [] work;
1204  if (info)
1205  {
1206  mfem::err << "DenseMatrix::SingularValues : info = " << info << endl;
1207  mfem_error();
1208  }
1209 #else
1210  // compiling without lapack
1211  mfem_error("DenseMatrix::SingularValues");
1212 #endif
1213 }
1214 
1215 int DenseMatrix::Rank(double tol) const
1216 {
1217  int rank=0;
1218  Vector sv(min(Height(), Width()));
1219  SingularValues(sv);
1220 
1221  for (int i=0; i < sv.Size(); ++i)
1222  if (sv(i) >= tol)
1223  {
1224  ++rank;
1225  }
1226 
1227  return rank;
1228 }
1229 
1230 static const double sqrt_1_eps = sqrt(1./numeric_limits<double>::epsilon());
1231 
1232 inline void Eigenvalues2S(const double &d12, double &d1, double &d2)
1233 {
1234  if (d12 != 0.)
1235  {
1236  // "The Symmetric Eigenvalue Problem", B. N. Parlett, pp.189-190
1237  double t, zeta = (d2 - d1)/(2*d12); // inf/inf from overflows?
1238  if (fabs(zeta) < sqrt_1_eps)
1239  {
1240  t = d12*copysign(1./(fabs(zeta) + sqrt(1. + zeta*zeta)), zeta);
1241  }
1242  else
1243  {
1244  t = d12*copysign(0.5/fabs(zeta), zeta);
1245  }
1246  d1 -= t;
1247  d2 += t;
1248  }
1249 }
1250 
1251 inline void Eigensystem2S(const double &d12, double &d1, double &d2,
1252  double &c, double &s)
1253 {
1254  if (d12 == 0.)
1255  {
1256  c = 1.;
1257  s = 0.;
1258  }
1259  else
1260  {
1261  // "The Symmetric Eigenvalue Problem", B. N. Parlett, pp.189-190
1262  double t, zeta = (d2 - d1)/(2*d12);
1263  if (fabs(zeta) < sqrt_1_eps)
1264  {
1265  t = copysign(1./(fabs(zeta) + sqrt(1. + zeta*zeta)), zeta);
1266  }
1267  else
1268  {
1269  t = copysign(0.5/fabs(zeta), zeta);
1270  }
1271  // c = 1./sqrt(1. + t*t);
1272  c = sqrt(1./(1. + t*t));
1273  s = c*t;
1274  t *= d12;
1275  d1 -= t;
1276  d2 += t;
1277  }
1278 }
1279 
1281  const double &x1, const double &x2, const double &x3,
1282  double &n1, double &n2, double &n3)
1283 {
1284  double m, t, r;
1285 
1286  m = fabs(x1);
1287  r = x2/m;
1288  t = 1. + r*r;
1289  r = x3/m;
1290  t = sqrt(1./(t + r*r));
1291  n1 = copysign(t, x1);
1292  t /= m;
1293  n2 = x2*t;
1294  n3 = x3*t;
1295 }
1296 
1297 inline void vec_normalize3(const double &x1, const double &x2, const double &x3,
1298  double &n1, double &n2, double &n3)
1299 {
1300  // should work ok when xk is the same as nk for some or all k
1301 
1302  if (fabs(x1) >= fabs(x2))
1303  {
1304  if (fabs(x1) >= fabs(x3))
1305  {
1306  if (x1 != 0.)
1307  {
1308  vec_normalize3_aux(x1, x2, x3, n1, n2, n3);
1309  }
1310  else
1311  {
1312  n1 = n2 = n3 = 0.;
1313  }
1314  return;
1315  }
1316  }
1317  else if (fabs(x2) >= fabs(x3))
1318  {
1319  vec_normalize3_aux(x2, x1, x3, n2, n1, n3);
1320  return;
1321  }
1322  vec_normalize3_aux(x3, x1, x2, n3, n1, n2);
1323 }
1324 
1325 inline bool KernelVector2G(
1326  const int &mode,
1327  double &d1, double &d12, double &d21, double &d2)
1328 {
1329  // Find a vector (z1,z2) in the "near"-kernel of the matrix
1330  // | d1 d12 |
1331  // | d21 d2 |
1332  // using QR factorization.
1333  // The vector (z1,z2) is returned in (d1,d2). Return 'true' if the matrix
1334  // is zero without setting (d1,d2).
1335  // Note: in the current implementation |z1| + |z2| = 1.
1336 
1337  // l1-norms of the columns
1338  double n1 = fabs(d1) + fabs(d21);
1339  double n2 = fabs(d2) + fabs(d12);
1340 
1341  bool swap_columns = (n2 > n1);
1342  double mu;
1343 
1344  if (!swap_columns)
1345  {
1346  if (n1 == 0.)
1347  {
1348  return true;
1349  }
1350 
1351  if (mode == 0) // eliminate the larger entry in the column
1352  {
1353  if (fabs(d1) > fabs(d21))
1354  {
1355  Swap(d1, d21);
1356  Swap(d12, d2);
1357  }
1358  }
1359  else // eliminate the smaller entry in the column
1360  {
1361  if (fabs(d1) < fabs(d21))
1362  {
1363  Swap(d1, d21);
1364  Swap(d12, d2);
1365  }
1366  }
1367  }
1368  else
1369  {
1370  // n2 > n1, swap columns 1 and 2
1371  if (mode == 0) // eliminate the larger entry in the column
1372  {
1373  if (fabs(d12) > fabs(d2))
1374  {
1375  Swap(d1, d2);
1376  Swap(d12, d21);
1377  }
1378  else
1379  {
1380  Swap(d1, d12);
1381  Swap(d21, d2);
1382  }
1383  }
1384  else // eliminate the smaller entry in the column
1385  {
1386  if (fabs(d12) < fabs(d2))
1387  {
1388  Swap(d1, d2);
1389  Swap(d12, d21);
1390  }
1391  else
1392  {
1393  Swap(d1, d12);
1394  Swap(d21, d2);
1395  }
1396  }
1397  }
1398 
1399  n1 = hypot(d1, d21);
1400 
1401  if (d21 != 0.)
1402  {
1403  // v = (n1, n2)^t, |v| = 1
1404  // Q = I - 2 v v^t, Q (d1, d21)^t = (mu, 0)^t
1405  mu = copysign(n1, d1);
1406  n1 = -d21*(d21/(d1 + mu)); // = d1 - mu
1407  d1 = mu;
1408  // normalize (n1,d21) to avoid overflow/underflow
1409  // normalize (n1,d21) by the max-norm to avoid the sqrt call
1410  if (fabs(n1) <= fabs(d21))
1411  {
1412  // (n1,n2) <-- (n1/d21,1)
1413  n1 = n1/d21;
1414  mu = (2./(1. + n1*n1))*(n1*d12 + d2);
1415  d2 = d2 - mu;
1416  d12 = d12 - mu*n1;
1417  }
1418  else
1419  {
1420  // (n1,n2) <-- (1,d21/n1)
1421  n2 = d21/n1;
1422  mu = (2./(1. + n2*n2))*(d12 + n2*d2);
1423  d2 = d2 - mu*n2;
1424  d12 = d12 - mu;
1425  }
1426  }
1427 
1428  // Solve:
1429  // | d1 d12 | | z1 | = | 0 |
1430  // | 0 d2 | | z2 | | 0 |
1431 
1432  // choose (z1,z2) to minimize |d1*z1 + d12*z2| + |d2*z2|
1433  // under the condition |z1| + |z2| = 1, z2 >= 0 (for uniqueness)
1434  // set t = z1, z2 = 1 - |t|, -1 <= t <= 1
1435  // objective function is:
1436  // |d1*t + d12*(1 - |t|)| + |d2|*(1 - |t|) -- piecewise linear with
1437  // possible minima are -1,0,1,t1 where t1: d1*t1 + d12*(1 - |t1|) = 0
1438  // values: @t=+/-1 -> |d1|, @t=0 -> |n1| + |d2|, @t=t1 -> |d2|*(1 - |t1|)
1439 
1440  // evaluate z2 @t=t1
1441  mu = -d12/d1;
1442  // note: |mu| <= 1, if using l2-norm for column pivoting
1443  // |mu| <= sqrt(2), if using l1-norm
1444  n2 = 1./(1. + fabs(mu));
1445  // check if |d1|<=|d2|*z2
1446  if (fabs(d1) <= n2*fabs(d2))
1447  {
1448  d2 = 0.;
1449  d1 = 1.;
1450  }
1451  else
1452  {
1453  d2 = n2;
1454  // d1 = (n2 < 0.5) ? copysign(1. - n2, mu) : mu*n2;
1455  d1 = mu*n2;
1456  }
1457 
1458  if (swap_columns)
1459  {
1460  Swap(d1, d2);
1461  }
1462 
1463  return false;
1464 }
1465 
1467  const int &mode,
1468  double &d1, double &d2, double &d3, double &c12, double &c13, double &c23,
1469  double &c21, double &c31, double &c32)
1470 {
1471  int kdim;
1472  double mu, n1, n2, n3, s1, s2, s3;
1473 
1474  s1 = hypot(c21, c31);
1475  n1 = hypot(d1, s1);
1476 
1477  if (s1 != 0.)
1478  {
1479  // v = (s1, s2, s3)^t, |v| = 1
1480  // Q = I - 2 v v^t, Q (d1, c12, c13)^t = (mu, 0, 0)^t
1481  mu = copysign(n1, d1);
1482  n1 = -s1*(s1/(d1 + mu)); // = d1 - mu
1483  d1 = mu;
1484 
1485  // normalize (n1,c21,c31) to avoid overflow/underflow
1486  // normalize (n1,c21,c31) by the max-norm to avoid the sqrt call
1487  if (fabs(n1) >= fabs(c21))
1488  {
1489  if (fabs(n1) >= fabs(c31))
1490  {
1491  // n1 is max, (s1,s2,s3) <-- (1,c21/n1,c31/n1)
1492  s2 = c21/n1;
1493  s3 = c31/n1;
1494  mu = 2./(1. + s2*s2 + s3*s3);
1495  n2 = mu*(c12 + s2*d2 + s3*c32);
1496  n3 = mu*(c13 + s2*c23 + s3*d3);
1497  c12 = c12 - n2;
1498  d2 = d2 - s2*n2;
1499  c32 = c32 - s3*n2;
1500  c13 = c13 - n3;
1501  c23 = c23 - s2*n3;
1502  d3 = d3 - s3*n3;
1503  goto done_column_1;
1504  }
1505  }
1506  else if (fabs(c21) >= fabs(c31))
1507  {
1508  // c21 is max, (s1,s2,s3) <-- (n1/c21,1,c31/c21)
1509  s1 = n1/c21;
1510  s3 = c31/c21;
1511  mu = 2./(1. + s1*s1 + s3*s3);
1512  n2 = mu*(s1*c12 + d2 + s3*c32);
1513  n3 = mu*(s1*c13 + c23 + s3*d3);
1514  c12 = c12 - s1*n2;
1515  d2 = d2 - n2;
1516  c32 = c32 - s3*n2;
1517  c13 = c13 - s1*n3;
1518  c23 = c23 - n3;
1519  d3 = d3 - s3*n3;
1520  goto done_column_1;
1521  }
1522  // c31 is max, (s1,s2,s3) <-- (n1/c31,c21/c31,1)
1523  s1 = n1/c31;
1524  s2 = c21/c31;
1525  mu = 2./(1. + s1*s1 + s2*s2);
1526  n2 = mu*(s1*c12 + s2*d2 + c32);
1527  n3 = mu*(s1*c13 + s2*c23 + d3);
1528  c12 = c12 - s1*n2;
1529  d2 = d2 - s2*n2;
1530  c32 = c32 - n2;
1531  c13 = c13 - s1*n3;
1532  c23 = c23 - s2*n3;
1533  d3 = d3 - n3;
1534  }
1535 
1536 done_column_1:
1537 
1538  // Solve:
1539  // | d2 c23 | | z2 | = | 0 |
1540  // | c32 d3 | | z3 | | 0 |
1541  if (KernelVector2G(mode, d2, c23, c32, d3))
1542  {
1543  // Have two solutions:
1544  // two vectors in the kernel are P (-c12/d1, 1, 0)^t and
1545  // P (-c13/d1, 0, 1)^t where P is the permutation matrix swapping
1546  // entries 1 and col.
1547 
1548  // A vector orthogonal to both these vectors is P (1, c12/d1, c13/d1)^t
1549  d2 = c12/d1;
1550  d3 = c13/d1;
1551  d1 = 1.;
1552  kdim = 2;
1553  }
1554  else
1555  {
1556  // solve for z1:
1557  // note: |z1| <= a since |z2| + |z3| = 1, and
1558  // max{|c12|,|c13|} <= max{norm(col. 2),norm(col. 3)}
1559  // <= norm(col. 1) <= a |d1|
1560  // a = 1, if using l2-norm for column pivoting
1561  // a = sqrt(3), if using l1-norm
1562  d1 = -(c12*d2 + c13*d3)/d1;
1563  kdim = 1;
1564  }
1565 
1566  vec_normalize3(d1, d2, d3, d1, d2, d3);
1567 
1568  return kdim;
1569 }
1570 
1571 inline int KernelVector3S(
1572  const int &mode,
1573  const double &d12, const double &d13, const double &d23,
1574  double &d1, double &d2, double &d3)
1575 {
1576  // Find a unit vector (z1,z2,z3) in the "near"-kernel of the matrix
1577  // | d1 d12 d13 |
1578  // | d12 d2 d23 |
1579  // | d13 d23 d3 |
1580  // using QR factorization.
1581  // The vector (z1,z2,z3) is returned in (d1,d2,d3).
1582  // Returns the dimension of the kernel, kdim, but never zero.
1583  // - if kdim == 3, then (d1,d2,d3) is not defined,
1584  // - if kdim == 2, then (d1,d2,d3) is a vector orthogonal to the kernel,
1585  // - otherwise kdim == 1 and (d1,d2,d3) is a vector in the "near"-kernel.
1586 
1587  double c12 = d12, c13 = d13, c23 = d23;
1588  double c21, c31, c32;
1589  int col, row;
1590 
1591  // l1-norms of the columns:
1592  c32 = fabs(d1) + fabs(c12) + fabs(c13);
1593  c31 = fabs(d2) + fabs(c12) + fabs(c23);
1594  c21 = fabs(d3) + fabs(c13) + fabs(c23);
1595 
1596  // column pivoting: choose the column with the largest norm
1597  if (c32 >= c21)
1598  {
1599  col = (c32 >= c31) ? 1 : 2;
1600  }
1601  else
1602  {
1603  col = (c31 >= c21) ? 2 : 3;
1604  }
1605  switch (col)
1606  {
1607  case 1:
1608  if (c32 == 0.) // zero matrix
1609  {
1610  return 3;
1611  }
1612  break;
1613 
1614  case 2:
1615  if (c31 == 0.) // zero matrix
1616  {
1617  return 3;
1618  }
1619  Swap(c13, c23);
1620  Swap(d1, d2);
1621  break;
1622 
1623  case 3:
1624  if (c21 == 0.) // zero matrix
1625  {
1626  return 3;
1627  }
1628  Swap(c12, c23);
1629  Swap(d1, d3);
1630  }
1631 
1632  // row pivoting depending on 'mode'
1633  if (mode == 0)
1634  {
1635  if (fabs(d1) <= fabs(c13))
1636  {
1637  row = (fabs(d1) <= fabs(c12)) ? 1 : 2;
1638  }
1639  else
1640  {
1641  row = (fabs(c12) <= fabs(c13)) ? 2 : 3;
1642  }
1643  }
1644  else
1645  {
1646  if (fabs(d1) >= fabs(c13))
1647  {
1648  row = (fabs(d1) >= fabs(c12)) ? 1 : 2;
1649  }
1650  else
1651  {
1652  row = (fabs(c12) >= fabs(c13)) ? 2 : 3;
1653  }
1654  }
1655  switch (row)
1656  {
1657  case 1:
1658  c21 = c12;
1659  c31 = c13;
1660  c32 = c23;
1661  break;
1662 
1663  case 2:
1664  c21 = d1;
1665  c31 = c13;
1666  c32 = c23;
1667  d1 = c12;
1668  c12 = d2;
1669  d2 = d1;
1670  c13 = c23;
1671  c23 = c31;
1672  break;
1673 
1674  case 3:
1675  c21 = c12;
1676  c31 = d1;
1677  c32 = c12;
1678  d1 = c13;
1679  c12 = c23;
1680  c13 = d3;
1681  d3 = d1;
1682  }
1683 
1684  row = KernelVector3G_aux(mode, d1, d2, d3, c12, c13, c23, c21, c31, c32);
1685  // row is kdim
1686 
1687  switch (col)
1688  {
1689  case 2:
1690  Swap(d1, d2);
1691  break;
1692 
1693  case 3:
1694  Swap(d1, d3);
1695  }
1696 
1697  return row;
1698 }
1699 
1700 inline int Reduce3S(
1701  const int &mode,
1702  double &d1, double &d2, double &d3, double &d12, double &d13, double &d23,
1703  double &z1, double &z2, double &z3, double &v1, double &v2, double &v3,
1704  double &g)
1705 {
1706  // Given the matrix
1707  // | d1 d12 d13 |
1708  // A = | d12 d2 d23 |
1709  // | d13 d23 d3 |
1710  // and a unit eigenvector z=(z1,z2,z3), transform the matrix A into the
1711  // matrix B = Q P A P Q that has the form
1712  // | b1 0 0 |
1713  // B = Q P A P Q = | 0 b2 b23 |
1714  // | 0 b23 b3 |
1715  // where P is the permutation matrix switching entries 1 and k, and
1716  // Q is the reflection matrix Q = I - g v v^t, defined by: set y = P z and
1717  // v = c(y - e_1); if y = e_1, then v = 0 and Q = I.
1718  // Note: Q y = e_1, Q e_1 = y ==> Q P A P Q e_1 = ... = lambda e_1.
1719  // The entries (b1,b2,b3,b23) are returned in (d1,d2,d3,d23), and the
1720  // return value of the function is k. The variable g = 2/(v1^2+v2^2+v3^3).
1721 
1722  int k;
1723  double s, w1, w2, w3;
1724 
1725  if (mode == 0)
1726  {
1727  // choose k such that z^t e_k = zk has the smallest absolute value, i.e.
1728  // the angle between z and e_k is closest to pi/2
1729  if (fabs(z1) <= fabs(z3))
1730  {
1731  k = (fabs(z1) <= fabs(z2)) ? 1 : 2;
1732  }
1733  else
1734  {
1735  k = (fabs(z2) <= fabs(z3)) ? 2 : 3;
1736  }
1737  }
1738  else
1739  {
1740  // choose k such that zk is the largest by absolute value
1741  if (fabs(z1) >= fabs(z3))
1742  {
1743  k = (fabs(z1) >= fabs(z2)) ? 1 : 2;
1744  }
1745  else
1746  {
1747  k = (fabs(z2) >= fabs(z3)) ? 2 : 3;
1748  }
1749  }
1750  switch (k)
1751  {
1752  case 2:
1753  Swap(d13, d23);
1754  Swap(d1, d2);
1755  Swap(z1, z2);
1756  break;
1757 
1758  case 3:
1759  Swap(d12, d23);
1760  Swap(d1, d3);
1761  Swap(z1, z3);
1762  }
1763 
1764  s = hypot(z2, z3);
1765 
1766  if (s == 0.)
1767  {
1768  // s can not be zero, if zk is the smallest (mode == 0)
1769  v1 = v2 = v3 = 0.;
1770  g = 1.;
1771  }
1772  else
1773  {
1774  g = copysign(1., z1);
1775  v1 = -s*(s/(z1 + g)); // = z1 - g
1776  // normalize (v1,z2,z3) by its max-norm, avoiding the sqrt call
1777  g = fabs(v1);
1778  if (fabs(z2) > g) { g = fabs(z2); }
1779  if (fabs(z3) > g) { g = fabs(z3); }
1780  v1 = v1/g;
1781  v2 = z2/g;
1782  v3 = z3/g;
1783  g = 2./(v1*v1 + v2*v2 + v3*v3);
1784 
1785  // Compute Q A Q = A - v w^t - w v^t, where
1786  // w = u - (g/2)(v^t u) v, and u = g A v
1787  // set w = g A v
1788  w1 = g*( d1*v1 + d12*v2 + d13*v3);
1789  w2 = g*(d12*v1 + d2*v2 + d23*v3);
1790  w3 = g*(d13*v1 + d23*v2 + d3*v3);
1791  // w := w - (g/2)(v^t w) v
1792  s = (g/2)*(v1*w1 + v2*w2 + v3*w3);
1793  w1 -= s*v1;
1794  w2 -= s*v2;
1795  w3 -= s*v3;
1796  // dij -= vi*wj + wi*vj
1797  d1 -= 2*v1*w1;
1798  d2 -= 2*v2*w2;
1799  d23 -= v2*w3 + v3*w2;
1800  d3 -= 2*v3*w3;
1801  // compute the offdiagonal entries on the first row/column of B which
1802  // should be zero (for debugging):
1803 #if 0
1804  s = d12 - v1*w2 - v2*w1; // b12 = 0
1805  s = d13 - v1*w3 - v3*w1; // b13 = 0
1806 #endif
1807  }
1808 
1809  switch (k)
1810  {
1811  case 2:
1812  Swap(z1, z2);
1813  break;
1814 
1815  case 3:
1816  Swap(z1, z3);
1817  }
1818 
1819  return k;
1820 }
1821 
1822 inline void GetScalingFactor(const double &d_max, double &mult)
1823 {
1824  int d_exp;
1825  if (d_max > 0.)
1826  {
1827  mult = frexp(d_max, &d_exp);
1828  if (d_exp == numeric_limits<double>::max_exponent)
1829  {
1830  mult *= numeric_limits<double>::radix;
1831  }
1832  mult = d_max/mult;
1833  }
1834  else
1835  {
1836  mult = 1.;
1837  }
1838  // mult = 2^d_exp is such that d_max/mult is in [0.5,1)
1839  // or in other words d_max is in the interval [0.5,1)*mult
1840 }
1841 
1842 double DenseMatrix::CalcSingularvalue(const int i) const
1843 {
1844  MFEM_ASSERT(Height() == Width() && Height() > 0 && Height() < 4,
1845  "The matrix must be square and sized 1, 2, or 3 to compute the"
1846  " singular values."
1847  << " Height() = " << Height()
1848  << ", Width() = " << Width());
1849 
1850  const int n = Height();
1851  const double *d = data;
1852 
1853  if (n == 1)
1854  {
1855  return d[0];
1856  }
1857  else if (n == 2)
1858  {
1859  double d0, d1, d2, d3;
1860  d0 = d[0];
1861  d1 = d[1];
1862  d2 = d[2];
1863  d3 = d[3];
1864  double mult;
1865  {
1866  double d_max = fabs(d0);
1867  if (d_max < fabs(d1)) { d_max = fabs(d1); }
1868  if (d_max < fabs(d2)) { d_max = fabs(d2); }
1869  if (d_max < fabs(d3)) { d_max = fabs(d3); }
1870 
1871  GetScalingFactor(d_max, mult);
1872  }
1873  d0 /= mult;
1874  d1 /= mult;
1875  d2 /= mult;
1876  d3 /= mult;
1877  // double b11 = d[0]*d[0] + d[1]*d[1];
1878  // double b12 = d[0]*d[2] + d[1]*d[3];
1879  // double b22 = d[2]*d[2] + d[3]*d[3];
1880  // t = 0.5*(a+b).(a-b) = 0.5*(|a|^2-|b|^2)
1881  // with a,b - the columns of (*this)
1882  // double t = 0.5*(b11 - b22);
1883  double t = 0.5*((d0+d2)*(d0-d2)+(d1-d3)*(d1+d3));
1884  // double s = sqrt(0.5*(b11 + b22) + sqrt(t*t + b12*b12));
1885  double s = d0*d2 + d1*d3;
1886  s = sqrt(0.5*(d0*d0 + d1*d1 + d2*d2 + d3*d3) + sqrt(t*t + s*s));
1887  if (s == 0.0)
1888  {
1889  return 0.0;
1890  }
1891  t = fabs(d0*d3 - d1*d2) / s;
1892  if (t > s)
1893  {
1894  if (i == 0)
1895  {
1896  return t*mult;
1897  }
1898  return s*mult;
1899  }
1900  if (i == 0)
1901  {
1902  return s*mult;
1903  }
1904  return t*mult;
1905  }
1906  else
1907  {
1908  double d0, d1, d2, d3, d4, d5, d6, d7, d8;
1909  d0 = d[0]; d3 = d[3]; d6 = d[6];
1910  d1 = d[1]; d4 = d[4]; d7 = d[7];
1911  d2 = d[2]; d5 = d[5]; d8 = d[8];
1912  double mult;
1913  {
1914  double d_max = fabs(d0);
1915  if (d_max < fabs(d1)) { d_max = fabs(d1); }
1916  if (d_max < fabs(d2)) { d_max = fabs(d2); }
1917  if (d_max < fabs(d3)) { d_max = fabs(d3); }
1918  if (d_max < fabs(d4)) { d_max = fabs(d4); }
1919  if (d_max < fabs(d5)) { d_max = fabs(d5); }
1920  if (d_max < fabs(d6)) { d_max = fabs(d6); }
1921  if (d_max < fabs(d7)) { d_max = fabs(d7); }
1922  if (d_max < fabs(d8)) { d_max = fabs(d8); }
1923 
1924  GetScalingFactor(d_max, mult);
1925  }
1926 
1927  d0 /= mult; d1 /= mult; d2 /= mult;
1928  d3 /= mult; d4 /= mult; d5 /= mult;
1929  d6 /= mult; d7 /= mult; d8 /= mult;
1930 
1931  double b11 = d0*d0 + d1*d1 + d2*d2;
1932  double b12 = d0*d3 + d1*d4 + d2*d5;
1933  double b13 = d0*d6 + d1*d7 + d2*d8;
1934  double b22 = d3*d3 + d4*d4 + d5*d5;
1935  double b23 = d3*d6 + d4*d7 + d5*d8;
1936  double b33 = d6*d6 + d7*d7 + d8*d8;
1937 
1938  // double a, b, c;
1939  // a = -(b11 + b22 + b33);
1940  // b = b11*(b22 + b33) + b22*b33 - b12*b12 - b13*b13 - b23*b23;
1941  // c = b11*(b23*b23 - b22*b33) + b12*(b12*b33 - 2*b13*b23) + b13*b13*b22;
1942 
1943  // double Q = (a * a - 3 * b) / 9;
1944  // double Q = (b12*b12 + b13*b13 + b23*b23 +
1945  // ((b11 - b22)*(b11 - b22) +
1946  // (b11 - b33)*(b11 - b33) +
1947  // (b22 - b33)*(b22 - b33))/6)/3;
1948  // Q = (3*(b12^2 + b13^2 + b23^2) +
1949  // ((b11 - b22)^2 + (b11 - b33)^2 + (b22 - b33)^2)/2)/9
1950  // or
1951  // Q = (1/6)*|B-tr(B)/3|_F^2
1952  // Q >= 0 and
1953  // Q = 0 <==> B = scalar * I
1954  // double R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
1955  double aa = (b11 + b22 + b33)/3; // aa = tr(B)/3
1956  double c1, c2, c3;
1957  // c1 = b11 - aa; // ((b11 - b22) + (b11 - b33))/3
1958  // c2 = b22 - aa; // ((b22 - b11) + (b22 - b33))/3
1959  // c3 = b33 - aa; // ((b33 - b11) + (b33 - b22))/3
1960  {
1961  double b11_b22 = ((d0-d3)*(d0+d3)+(d1-d4)*(d1+d4)+(d2-d5)*(d2+d5));
1962  double b22_b33 = ((d3-d6)*(d3+d6)+(d4-d7)*(d4+d7)+(d5-d8)*(d5+d8));
1963  double b33_b11 = ((d6-d0)*(d6+d0)+(d7-d1)*(d7+d1)+(d8-d2)*(d8+d2));
1964  c1 = (b11_b22 - b33_b11)/3;
1965  c2 = (b22_b33 - b11_b22)/3;
1966  c3 = (b33_b11 - b22_b33)/3;
1967  }
1968  double Q, R;
1969  Q = (2*(b12*b12 + b13*b13 + b23*b23) + c1*c1 + c2*c2 + c3*c3)/6;
1970  R = (c1*(b23*b23 - c2*c3)+ b12*(b12*c3 - 2*b13*b23) +b13*b13*c2)/2;
1971  // R = (-1/2)*det(B-(tr(B)/3)*I)
1972  // Note: 54*(det(S))^2 <= |S|_F^6, when S^t=S and tr(S)=0, S is 3x3
1973  // Therefore: R^2 <= Q^3
1974 
1975  if (Q <= 0.) { ; }
1976 
1977  // else if (fabs(R) >= sqrtQ3)
1978  // {
1979  // double det = (d[0] * (d[4] * d[8] - d[5] * d[7]) +
1980  // d[3] * (d[2] * d[7] - d[1] * d[8]) +
1981  // d[6] * (d[1] * d[5] - d[2] * d[4]));
1982  //
1983  // if (R > 0.)
1984  // {
1985  // if (i == 2)
1986  // // aa -= 2*sqrtQ;
1987  // return fabs(det)/(aa + sqrtQ);
1988  // else
1989  // aa += sqrtQ;
1990  // }
1991  // else
1992  // {
1993  // if (i != 0)
1994  // aa -= sqrtQ;
1995  // // aa = fabs(det)/sqrt(aa + 2*sqrtQ);
1996  // else
1997  // aa += 2*sqrtQ;
1998  // }
1999  // }
2000 
2001  else
2002  {
2003  double sqrtQ = sqrt(Q);
2004  double sqrtQ3 = Q*sqrtQ;
2005  // double sqrtQ3 = sqrtQ*sqrtQ*sqrtQ;
2006  // double sqrtQ3 = pow(Q, 1.5);
2007  double r;
2008 
2009  if (fabs(R) >= sqrtQ3)
2010  {
2011  if (R < 0.)
2012  {
2013  // R = -1.;
2014  r = 2*sqrtQ;
2015  }
2016  else
2017  {
2018  // R = 1.;
2019  r = -2*sqrtQ;
2020  }
2021  }
2022  else
2023  {
2024  R = R/sqrtQ3;
2025 
2026  // if (fabs(R) <= 0.95)
2027  if (fabs(R) <= 0.9)
2028  {
2029  if (i == 2)
2030  {
2031  aa -= 2*sqrtQ*cos(acos(R)/3); // min
2032  }
2033  else if (i == 0)
2034  {
2035  aa -= 2*sqrtQ*cos((acos(R) + 2.0*M_PI)/3); // max
2036  }
2037  else
2038  {
2039  aa -= 2*sqrtQ*cos((acos(R) - 2.0*M_PI)/3); // mid
2040  }
2041  goto have_aa;
2042  }
2043 
2044  if (R < 0.)
2045  {
2046  r = -2*sqrtQ*cos((acos(R) + 2.0*M_PI)/3); // max
2047  if (i == 0)
2048  {
2049  aa += r;
2050  goto have_aa;
2051  }
2052  }
2053  else
2054  {
2055  r = -2*sqrtQ*cos(acos(R)/3); // min
2056  if (i == 2)
2057  {
2058  aa += r;
2059  goto have_aa;
2060  }
2061  }
2062  }
2063 
2064  // (tr(B)/3 + r) is the root which is separated from the other
2065  // two roots which are close to each other when |R| is close to 1
2066 
2067  c1 -= r;
2068  c2 -= r;
2069  c3 -= r;
2070  // aa += r;
2071 
2072  // Type of Householder reflections: z --> mu ek, where k is the index
2073  // of the entry in z with:
2074  // mode == 0: smallest absolute value --> angle closest to pi/2
2075  // (eliminate large entries)
2076  // mode == 1: largest absolute value --> angle farthest from pi/2
2077  // (eliminate small entries)
2078  const int mode = 1;
2079 
2080  // Find a unit vector z = (z1,z2,z3) in the "near"-kernel of
2081  // | c1 b12 b13 |
2082  // | b12 c2 b23 | = B - aa*I
2083  // | b13 b23 c3 |
2084  // This vector is also an eigenvector for B corresponding to aa
2085  // The vector z overwrites (c1,c2,c3).
2086  switch (KernelVector3S(mode, b12, b13, b23, c1, c2, c3))
2087  {
2088  case 3:
2089  aa += r;
2090  goto have_aa;
2091  case 2:
2092  // ok, continue with the returned vector orthogonal to the kernel
2093  case 1:
2094  // ok, continue with the returned vector in the "near"-kernel
2095  ;
2096  }
2097 
2098  // Using the eigenvector c = (c1,c2,c3) to transform B into
2099  // | b11 0 0 |
2100  // B <-- Q P B P Q = | 0 b22 b23 |
2101  // | 0 b23 b33 |
2102  double v1, v2, v3, g;
2103  Reduce3S(mode, b11, b22, b33, b12, b13, b23,
2104  c1, c2, c3, v1, v2, v3, g);
2105  // Q = I - g v v^t
2106  // P - permutation matrix switching rows and columns 1 and k
2107 
2108  // find the eigenvalues of
2109  // | b22 b23 |
2110  // | b23 b33 |
2111  Eigenvalues2S(b23, b22, b33);
2112 
2113  if (i == 2)
2114  {
2115  aa = std::min(std::min(b11, b22), b33);
2116  }
2117  else if (i == 1)
2118  {
2119  if (b11 <= b22)
2120  {
2121  aa = (b22 <= b33) ? b22 : std::max(b11, b33);
2122  }
2123  else
2124  {
2125  aa = (b11 <= b33) ? b11 : std::max(b33, b22);
2126  }
2127  }
2128  else
2129  {
2130  aa = std::max(std::max(b11, b22), b33);
2131  }
2132  }
2133 
2134  have_aa:
2135 
2136  return sqrt(fabs(aa))*mult; // take abs before we sort?
2137  }
2138 }
2139 
2140 void DenseMatrix::CalcEigenvalues(double *lambda, double *vec) const
2141 {
2142 #ifdef MFEM_DEBUG
2143  if (Height() != Width() || Height() < 2 || Height() > 3)
2144  {
2145  mfem_error("DenseMatrix::CalcEigenvalues");
2146  }
2147 #endif
2148 
2149  const int n = Height();
2150  const double *d = data;
2151 
2152  if (n == 2)
2153  {
2154  double d0 = d[0];
2155  double d2 = d[2]; // use the upper triangular entry
2156  double d3 = d[3];
2157 
2158  double c, s;
2159  Eigensystem2S(d2, d0, d3, c, s);
2160  if (d0 <= d3)
2161  {
2162  lambda[0] = d0;
2163  lambda[1] = d3;
2164  vec[0] = c;
2165  vec[1] = -s;
2166  vec[2] = s;
2167  vec[3] = c;
2168  }
2169  else
2170  {
2171  lambda[0] = d3;
2172  lambda[1] = d0;
2173  vec[0] = s;
2174  vec[1] = c;
2175  vec[2] = c;
2176  vec[3] = -s;
2177  }
2178  }
2179  else
2180  {
2181  double d11 = d[0];
2182  double d12 = d[3]; // use the upper triangular entries
2183  double d22 = d[4];
2184  double d13 = d[6];
2185  double d23 = d[7];
2186  double d33 = d[8];
2187 
2188  double mult;
2189  {
2190  double d_max = fabs(d11);
2191  if (d_max < fabs(d22)) { d_max = fabs(d22); }
2192  if (d_max < fabs(d33)) { d_max = fabs(d33); }
2193  if (d_max < fabs(d12)) { d_max = fabs(d12); }
2194  if (d_max < fabs(d13)) { d_max = fabs(d13); }
2195  if (d_max < fabs(d23)) { d_max = fabs(d23); }
2196 
2197  GetScalingFactor(d_max, mult);
2198  }
2199 
2200  d11 /= mult; d22 /= mult; d33 /= mult;
2201  d12 /= mult; d13 /= mult; d23 /= mult;
2202 
2203  double aa = (d11 + d22 + d33)/3; // aa = tr(A)/3
2204  double c1 = d11 - aa;
2205  double c2 = d22 - aa;
2206  double c3 = d33 - aa;
2207 
2208  double Q, R;
2209 
2210  Q = (2*(d12*d12 + d13*d13 + d23*d23) + c1*c1 + c2*c2 + c3*c3)/6;
2211  R = (c1*(d23*d23 - c2*c3)+ d12*(d12*c3 - 2*d13*d23) + d13*d13*c2)/2;
2212 
2213  if (Q <= 0.)
2214  {
2215  lambda[0] = lambda[1] = lambda[2] = aa;
2216  vec[0] = 1.; vec[3] = 0.; vec[6] = 0.;
2217  vec[1] = 0.; vec[4] = 1.; vec[7] = 0.;
2218  vec[2] = 0.; vec[5] = 0.; vec[8] = 1.;
2219  }
2220  else
2221  {
2222  double sqrtQ = sqrt(Q);
2223  double sqrtQ3 = Q*sqrtQ;
2224  // double sqrtQ3 = sqrtQ*sqrtQ*sqrtQ;
2225  // double sqrtQ3 = pow(Q, 1.5);
2226  double r;
2227  if (fabs(R) >= sqrtQ3)
2228  {
2229  if (R < 0.)
2230  {
2231  // R = -1.;
2232  r = 2*sqrtQ;
2233  }
2234  else
2235  {
2236  // R = 1.;
2237  r = -2*sqrtQ;
2238  }
2239  }
2240  else
2241  {
2242  R = R/sqrtQ3;
2243 
2244  if (R < 0.)
2245  {
2246  r = -2*sqrtQ*cos((acos(R) + 2.0*M_PI)/3); // max
2247  }
2248  else
2249  {
2250  r = -2*sqrtQ*cos(acos(R)/3); // min
2251  }
2252  }
2253 
2254  aa += r;
2255  c1 = d11 - aa;
2256  c2 = d22 - aa;
2257  c3 = d33 - aa;
2258 
2259  // Type of Householder reflections: z --> mu ek, where k is the index
2260  // of the entry in z with:
2261  // mode == 0: smallest absolute value --> angle closest to pi/2
2262  // mode == 1: largest absolute value --> angle farthest from pi/2
2263  // Observations:
2264  // mode == 0 produces better eigenvectors, less accurate eigenvalues?
2265  // mode == 1 produces better eigenvalues, less accurate eigenvectors?
2266  const int mode = 0;
2267 
2268  // Find a unit vector z = (z1,z2,z3) in the "near"-kernel of
2269  // | c1 d12 d13 |
2270  // | d12 c2 d23 | = A - aa*I
2271  // | d13 d23 c3 |
2272  // This vector is also an eigenvector for A corresponding to aa.
2273  // The vector z overwrites (c1,c2,c3).
2274  switch (KernelVector3S(mode, d12, d13, d23, c1, c2, c3))
2275  {
2276  case 3:
2277  // 'aa' is a triple eigenvalue
2278  lambda[0] = lambda[1] = lambda[2] = aa;
2279  vec[0] = 1.; vec[3] = 0.; vec[6] = 0.;
2280  vec[1] = 0.; vec[4] = 1.; vec[7] = 0.;
2281  vec[2] = 0.; vec[5] = 0.; vec[8] = 1.;
2282  goto done_3d;
2283 
2284  case 2:
2285  // ok, continue with the returned vector orthogonal to the kernel
2286  case 1:
2287  // ok, continue with the returned vector in the "near"-kernel
2288  ;
2289  }
2290 
2291  // Using the eigenvector c=(c1,c2,c3) transform A into
2292  // | d11 0 0 |
2293  // A <-- Q P A P Q = | 0 d22 d23 |
2294  // | 0 d23 d33 |
2295  double v1, v2, v3, g;
2296  int k = Reduce3S(mode, d11, d22, d33, d12, d13, d23,
2297  c1, c2, c3, v1, v2, v3, g);
2298  // Q = I - 2 v v^t
2299  // P - permutation matrix switching entries 1 and k
2300 
2301  // find the eigenvalues and eigenvectors for
2302  // | d22 d23 |
2303  // | d23 d33 |
2304  double c, s;
2305  Eigensystem2S(d23, d22, d33, c, s);
2306  // d22 <-> P Q (0, c, -s), d33 <-> P Q (0, s, c)
2307 
2308  double *vec_1, *vec_2, *vec_3;
2309  if (d11 <= d22)
2310  {
2311  if (d22 <= d33)
2312  {
2313  lambda[0] = d11; vec_1 = vec;
2314  lambda[1] = d22; vec_2 = vec + 3;
2315  lambda[2] = d33; vec_3 = vec + 6;
2316  }
2317  else if (d11 <= d33)
2318  {
2319  lambda[0] = d11; vec_1 = vec;
2320  lambda[1] = d33; vec_3 = vec + 3;
2321  lambda[2] = d22; vec_2 = vec + 6;
2322  }
2323  else
2324  {
2325  lambda[0] = d33; vec_3 = vec;
2326  lambda[1] = d11; vec_1 = vec + 3;
2327  lambda[2] = d22; vec_2 = vec + 6;
2328  }
2329  }
2330  else
2331  {
2332  if (d11 <= d33)
2333  {
2334  lambda[0] = d22; vec_2 = vec;
2335  lambda[1] = d11; vec_1 = vec + 3;
2336  lambda[2] = d33; vec_3 = vec + 6;
2337  }
2338  else if (d22 <= d33)
2339  {
2340  lambda[0] = d22; vec_2 = vec;
2341  lambda[1] = d33; vec_3 = vec + 3;
2342  lambda[2] = d11; vec_1 = vec + 6;
2343  }
2344  else
2345  {
2346  lambda[0] = d33; vec_3 = vec;
2347  lambda[1] = d22; vec_2 = vec + 3;
2348  lambda[2] = d11; vec_1 = vec + 6;
2349  }
2350  }
2351 
2352  vec_1[0] = c1;
2353  vec_1[1] = c2;
2354  vec_1[2] = c3;
2355  d22 = g*(v2*c - v3*s);
2356  d33 = g*(v2*s + v3*c);
2357  vec_2[0] = - v1*d22; vec_3[0] = - v1*d33;
2358  vec_2[1] = c - v2*d22; vec_3[1] = s - v2*d33;
2359  vec_2[2] = -s - v3*d22; vec_3[2] = c - v3*d33;
2360  switch (k)
2361  {
2362  case 2:
2363  Swap(vec_2[0], vec_2[1]);
2364  Swap(vec_3[0], vec_3[1]);
2365  break;
2366 
2367  case 3:
2368  Swap(vec_2[0], vec_2[2]);
2369  Swap(vec_3[0], vec_3[2]);
2370  }
2371  }
2372 
2373  done_3d:
2374  lambda[0] *= mult;
2375  lambda[1] *= mult;
2376  lambda[2] *= mult;
2377  }
2378 }
2379 
2380 void DenseMatrix::GetRow(int r, Vector &row) const
2381 {
2382  int m = Height();
2383  int n = Width();
2384  row.SetSize(n);
2385 
2386  const double* rp = data + r;
2387  double* vp = row.GetData();
2388 
2389  for (int i = 0; i < n; i++)
2390  {
2391  vp[i] = *rp;
2392  rp += m;
2393  }
2394 }
2395 
2396 void DenseMatrix::GetColumn(int c, Vector &col) const
2397 {
2398  int m = Height();
2399  col.SetSize(m);
2400 
2401  double *cp = data + c * m;
2402  double *vp = col.GetData();
2403 
2404  for (int i = 0; i < m; i++)
2405  {
2406  vp[i] = cp[i];
2407  }
2408 }
2409 
2411 {
2412  if (height != width)
2413  {
2414  mfem_error("DenseMatrix::GetDiag\n");
2415  }
2416  d.SetSize(height);
2417 
2418  for (int i = 0; i < height; ++i)
2419  {
2420  d(i) = (*this)(i,i);
2421  }
2422 }
2423 
2425 {
2426  if (height != width)
2427  {
2428  mfem_error("DenseMatrix::Getl1Diag\n");
2429  }
2430  l.SetSize(height);
2431 
2432  l = 0.0;
2433 
2434  for (int j = 0; j < width; ++j)
2435  for (int i = 0; i < height; ++i)
2436  {
2437  l(i) += fabs((*this)(i,j));
2438  }
2439 }
2440 
2442 {
2443  l.SetSize(height);
2444  for (int i = 0; i < height; i++)
2445  {
2446  double d = 0.0;
2447  for (int j = 0; j < width; j++)
2448  {
2449  d += operator()(i, j);
2450  }
2451  l(i) = d;
2452  }
2453 }
2454 
2455 void DenseMatrix::Diag(double c, int n)
2456 {
2457  SetSize(n);
2458 
2459  const int N = n*n;
2460  for (int i = 0; i < N; i++)
2461  {
2462  data[i] = 0.0;
2463  }
2464  for (int i = 0; i < n; i++)
2465  {
2466  data[i*(n+1)] = c;
2467  }
2468 }
2469 
2470 void DenseMatrix::Diag(double *diag, int n)
2471 {
2472  SetSize(n);
2473 
2474  int i, N = n*n;
2475  for (i = 0; i < N; i++)
2476  {
2477  data[i] = 0.0;
2478  }
2479  for (i = 0; i < n; i++)
2480  {
2481  data[i*(n+1)] = diag[i];
2482  }
2483 }
2484 
2486 {
2487  int i, j;
2488  double t;
2489 
2490  if (Width() == Height())
2491  {
2492  for (i = 0; i < Height(); i++)
2493  for (j = i+1; j < Width(); j++)
2494  {
2495  t = (*this)(i,j);
2496  (*this)(i,j) = (*this)(j,i);
2497  (*this)(j,i) = t;
2498  }
2499  }
2500  else
2501  {
2502  DenseMatrix T(*this,'t');
2503  (*this) = T;
2504  }
2505 }
2506 
2508 {
2509  SetSize(A.Width(),A.Height());
2510 
2511  for (int i = 0; i < Height(); i++)
2512  for (int j = 0; j < Width(); j++)
2513  {
2514  (*this)(i,j) = A(j,i);
2515  }
2516 }
2517 
2519 {
2520 #ifdef MFEM_DEBUG
2521  if (Width() != Height())
2522  {
2523  mfem_error("DenseMatrix::Symmetrize() : not a square matrix!");
2524  }
2525 #endif
2526 
2527  for (int i = 0; i < Height(); i++)
2528  for (int j = 0; j < i; j++)
2529  {
2530  double a = 0.5 * ((*this)(i,j) + (*this)(j,i));
2531  (*this)(j,i) = (*this)(i,j) = a;
2532  }
2533 }
2534 
2536 {
2537  for (int i = 0; i < Height(); i++)
2538  {
2539  double L = 0.0;
2540  for (int j = 0; j < Width(); j++)
2541  {
2542  L += (*this)(i, j);
2543  (*this)(i, j) = 0.0;
2544  }
2545  (*this)(i, i) = L;
2546  }
2547 }
2548 
2550 {
2551  int n = Height();
2552 
2553 #ifdef MFEM_DEBUG
2554  if ((Width() != 2 || curl.Width() != 1 || 2*n != curl.Height()) &&
2555  (Width() != 3 || curl.Width() != 3 || 3*n != curl.Height()))
2556  {
2557  mfem_error("DenseMatrix::GradToCurl(...)");
2558  }
2559 #endif
2560 
2561  if (Width() == 2)
2562  {
2563  for (int i = 0; i < n; i++)
2564  {
2565  // (x,y) is grad of Ui
2566  double x = (*this)(i,0);
2567  double y = (*this)(i,1);
2568 
2569  int j = i+n;
2570 
2571  // curl of (Ui,0)
2572  curl(i,0) = -y;
2573 
2574  // curl of (0,Ui)
2575  curl(j,0) = x;
2576  }
2577  }
2578  else
2579  {
2580  for (int i = 0; i < n; i++)
2581  {
2582  // (x,y,z) is grad of Ui
2583  double x = (*this)(i,0);
2584  double y = (*this)(i,1);
2585  double z = (*this)(i,2);
2586 
2587  int j = i+n;
2588  int k = j+n;
2589 
2590  // curl of (Ui,0,0)
2591  curl(i,0) = 0.;
2592  curl(i,1) = z;
2593  curl(i,2) = -y;
2594 
2595  // curl of (0,Ui,0)
2596  curl(j,0) = -z;
2597  curl(j,1) = 0.;
2598  curl(j,2) = x;
2599 
2600  // curl of (0,0,Ui)
2601  curl(k,0) = y;
2602  curl(k,1) = -x;
2603  curl(k,2) = 0.;
2604  }
2605  }
2606 }
2607 
2609 {
2610  MFEM_ASSERT(Width()*Height() == div.Size(), "incompatible Vector 'div'!");
2611 
2612  // div(dof*j+i) <-- (*this)(i,j)
2613 
2614  const int n = height * width;
2615  double *ddata = div.GetData();
2616 
2617  for (int i = 0; i < n; i++)
2618  {
2619  ddata[i] = data[i];
2620  }
2621 }
2622 
2623 void DenseMatrix::CopyRows(const DenseMatrix &A, int row1, int row2)
2624 {
2625  SetSize(row2 - row1 + 1, A.Width());
2626 
2627  for (int j = 0; j < Width(); j++)
2628  {
2629  for (int i = row1; i <= row2; i++)
2630  {
2631  (*this)(i-row1,j) = A(i,j);
2632  }
2633  }
2634 }
2635 
2636 void DenseMatrix::CopyCols(const DenseMatrix &A, int col1, int col2)
2637 {
2638  SetSize(A.Height(), col2 - col1 + 1);
2639 
2640  for (int j = col1; j <= col2; j++)
2641  {
2642  for (int i = 0; i < Height(); i++)
2643  {
2644  (*this)(i,j-col1) = A(i,j);
2645  }
2646  }
2647 }
2648 
2649 void DenseMatrix::CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
2650 {
2651  SetSize(m,n);
2652 
2653  for (int j = 0; j < n; j++)
2654  {
2655  for (int i = 0; i < m; i++)
2656  {
2657  (*this)(i,j) = A(Aro+i,Aco+j);
2658  }
2659  }
2660 }
2661 
2662 void DenseMatrix::CopyMN(const DenseMatrix &A, int row_offset, int col_offset)
2663 {
2664  double *v = A.data;
2665 
2666  for (int j = 0; j < A.Width(); j++)
2667  {
2668  for (int i = 0; i < A.Height(); i++)
2669  {
2670  (*this)(row_offset+i,col_offset+j) = *(v++);
2671  }
2672  }
2673 }
2674 
2675 void DenseMatrix::CopyMNt(const DenseMatrix &A, int row_offset, int col_offset)
2676 {
2677  double *v = A.data;
2678 
2679  for (int i = 0; i < A.Width(); i++)
2680  {
2681  for (int j = 0; j < A.Height(); j++)
2682  {
2683  (*this)(row_offset+i,col_offset+j) = *(v++);
2684  }
2685  }
2686 }
2687 
2688 void DenseMatrix::CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco,
2689  int row_offset, int col_offset)
2690 {
2691  MFEM_VERIFY(row_offset+m <= this->Height() && col_offset+n <= this->Width(),
2692  "this DenseMatrix is too small to accomodate the submatrix. "
2693  << "row_offset = " << row_offset
2694  << ", m = " << m
2695  << ", this->Height() = " << this->Height()
2696  << ", col_offset = " << col_offset
2697  << ", n = " << n
2698  << ", this->Width() = " << this->Width()
2699  );
2700  MFEM_VERIFY(Aro+m <= A.Height() && Aco+n <= A.Width(),
2701  "The A DenseMatrix is too small to accomodate the submatrix. "
2702  << "Aro = " << Aro
2703  << ", m = " << m
2704  << ", A.Height() = " << A.Height()
2705  << ", Aco = " << Aco
2706  << ", n = " << n
2707  << ", A.Width() = " << A.Width()
2708  );
2709 
2710  for (int j = 0; j < n; j++)
2711  {
2712  for (int i = 0; i < m; i++)
2713  {
2714  (*this)(row_offset+i,col_offset+j) = A(Aro+i,Aco+j);
2715  }
2716  }
2717 }
2718 
2719 void DenseMatrix::CopyMNDiag(double c, int n, int row_offset, int col_offset)
2720 {
2721  for (int i = 0; i < n; i++)
2722  {
2723  for (int j = i+1; j < n; j++)
2724  {
2725  (*this)(row_offset+i,col_offset+j) =
2726  (*this)(row_offset+j,col_offset+i) = 0.0;
2727  }
2728  }
2729 
2730  for (int i = 0; i < n; i++)
2731  {
2732  (*this)(row_offset+i,col_offset+i) = c;
2733  }
2734 }
2735 
2736 void DenseMatrix::CopyMNDiag(double *diag, int n, int row_offset,
2737  int col_offset)
2738 {
2739  for (int i = 0; i < n; i++)
2740  {
2741  for (int j = i+1; j < n; j++)
2742  {
2743  (*this)(row_offset+i,col_offset+j) =
2744  (*this)(row_offset+j,col_offset+i) = 0.0;
2745  }
2746  }
2747 
2748  for (int i = 0; i < n; i++)
2749  {
2750  (*this)(row_offset+i,col_offset+i) = diag[i];
2751  }
2752 }
2753 
2754 void DenseMatrix::CopyExceptMN(const DenseMatrix &A, int m, int n)
2755 {
2756  SetSize(A.Width()-1,A.Height()-1);
2757 
2758  int i, j, i_off = 0, j_off = 0;
2759 
2760  for (j = 0; j < A.Width(); j++)
2761  {
2762  if ( j == n )
2763  {
2764  j_off = 1;
2765  continue;
2766  }
2767  for (i = 0; i < A.Height(); i++)
2768  {
2769  if ( i == m )
2770  {
2771  i_off = 1;
2772  continue;
2773  }
2774  (*this)(i-i_off,j-j_off) = A(i,j);
2775  }
2776  i_off = 0;
2777  }
2778 }
2779 
2780 void DenseMatrix::AddMatrix(DenseMatrix &A, int ro, int co)
2781 {
2782  int h, ah, aw;
2783  double *p, *ap;
2784 
2785  h = Height();
2786  ah = A.Height();
2787  aw = A.Width();
2788 
2789 #ifdef MFEM_DEBUG
2790  if (co+aw > Width() || ro+ah > h)
2791  {
2792  mfem_error("DenseMatrix::AddMatrix(...) 1");
2793  }
2794 #endif
2795 
2796  p = data + ro + co * h;
2797  ap = A.data;
2798 
2799  for (int c = 0; c < aw; c++)
2800  {
2801  for (int r = 0; r < ah; r++)
2802  {
2803  p[r] += ap[r];
2804  }
2805  p += h;
2806  ap += ah;
2807  }
2808 }
2809 
2810 void DenseMatrix::AddMatrix(double a, const DenseMatrix &A, int ro, int co)
2811 {
2812  int h, ah, aw;
2813  double *p, *ap;
2814 
2815  h = Height();
2816  ah = A.Height();
2817  aw = A.Width();
2818 
2819 #ifdef MFEM_DEBUG
2820  if (co+aw > Width() || ro+ah > h)
2821  {
2822  mfem_error("DenseMatrix::AddMatrix(...) 2");
2823  }
2824 #endif
2825 
2826  p = data + ro + co * h;
2827  ap = A.data;
2828 
2829  for (int c = 0; c < aw; c++)
2830  {
2831  for (int r = 0; r < ah; r++)
2832  {
2833  p[r] += a * ap[r];
2834  }
2835  p += h;
2836  ap += ah;
2837  }
2838 }
2839 
2840 void DenseMatrix::AddToVector(int offset, Vector &v) const
2841 {
2842  const int n = height * width;
2843  double *vdata = v.GetData() + offset;
2844 
2845  for (int i = 0; i < n; i++)
2846  {
2847  vdata[i] += data[i];
2848  }
2849 }
2850 
2851 void DenseMatrix::GetFromVector(int offset, const Vector &v)
2852 {
2853  const int n = height * width;
2854  const double *vdata = v.GetData() + offset;
2855 
2856  for (int i = 0; i < n; i++)
2857  {
2858  data[i] = vdata[i];
2859  }
2860 }
2861 
2863 {
2864  const int n = Height();
2865 
2866 #ifdef MFEM_DEBUG
2867  if (dofs.Size() != n || Width() != n)
2868  {
2869  mfem_error("DenseMatrix::AdjustDofDirection(...)");
2870  }
2871 #endif
2872 
2873  int *dof = dofs;
2874  for (int i = 0; i < n-1; i++)
2875  {
2876  const int s = (dof[i] < 0) ? (-1) : (1);
2877  for (int j = i+1; j < n; j++)
2878  {
2879  const int t = (dof[j] < 0) ? (-s) : (s);
2880  if (t < 0)
2881  {
2882  (*this)(i,j) = -(*this)(i,j);
2883  (*this)(j,i) = -(*this)(j,i);
2884  }
2885  }
2886  }
2887 }
2888 
2889 void DenseMatrix::SetRow(int row, double value)
2890 {
2891  for (int j = 0; j < Width(); j++)
2892  {
2893  (*this)(row, j) = value;
2894  }
2895 }
2896 
2897 void DenseMatrix::SetCol(int col, double value)
2898 {
2899  for (int i = 0; i < Height(); i++)
2900  {
2901  (*this)(i, col) = value;
2902  }
2903 }
2904 
2905 void DenseMatrix::SetRow(int r, const Vector &row)
2906 {
2907  for (int j = 0; j < Width(); j++)
2908  {
2909  (*this)(r, j) = row[j];
2910  }
2911 }
2912 
2913 void DenseMatrix::SetCol(int c, const Vector &col)
2914 {
2915  for (int i = 0; i < Height(); i++)
2916  {
2917  (*this)(i, c) = col[i];
2918  }
2919 }
2920 
2921 void DenseMatrix::Threshold(double eps)
2922 {
2923  for (int col = 0; col < Width(); col++)
2924  {
2925  for (int row = 0; row < Height(); row++)
2926  {
2927  if (std::abs(operator()(row,col)) <= eps)
2928  {
2929  operator()(row,col) = 0.0;
2930  }
2931  }
2932  }
2933 }
2934 
2935 void DenseMatrix::Print(std::ostream &out, int width_) const
2936 {
2937  // save current output flags
2938  ios::fmtflags old_flags = out.flags();
2939  // output flags = scientific + show sign
2940  out << setiosflags(ios::scientific | ios::showpos);
2941  for (int i = 0; i < height; i++)
2942  {
2943  out << "[row " << i << "]\n";
2944  for (int j = 0; j < width; j++)
2945  {
2946  out << (*this)(i,j);
2947  if (j+1 == width || (j+1) % width_ == 0)
2948  {
2949  out << '\n';
2950  }
2951  else
2952  {
2953  out << ' ';
2954  }
2955  }
2956  }
2957  // reset output flags to original values
2958  out.flags(old_flags);
2959 }
2960 
2961 void DenseMatrix::PrintMatlab(std::ostream &out) const
2962 {
2963  // save current output flags
2964  ios::fmtflags old_flags = out.flags();
2965  // output flags = scientific + show sign
2966  out << setiosflags(ios::scientific | ios::showpos);
2967  for (int i = 0; i < height; i++)
2968  {
2969  for (int j = 0; j < width; j++)
2970  {
2971  out << (*this)(i,j);
2972  out << ' ';
2973  }
2974  out << "\n";
2975  }
2976  // reset output flags to original values
2977  out.flags(old_flags);
2978 }
2979 
2980 void DenseMatrix::PrintT(std::ostream &out, int width_) const
2981 {
2982  // save current output flags
2983  ios::fmtflags old_flags = out.flags();
2984  // output flags = scientific + show sign
2985  out << setiosflags(ios::scientific | ios::showpos);
2986  for (int j = 0; j < width; j++)
2987  {
2988  out << "[col " << j << "]\n";
2989  for (int i = 0; i < height; i++)
2990  {
2991  out << (*this)(i,j);
2992  if (i+1 == height || (i+1) % width_ == 0)
2993  {
2994  out << '\n';
2995  }
2996  else
2997  {
2998  out << ' ';
2999  }
3000  }
3001  }
3002  // reset output flags to original values
3003  out.flags(old_flags);
3004 }
3005 
3007 {
3008  DenseMatrix copy(*this), C(width);
3009  Invert();
3010  mfem::Mult(*this, copy, C);
3011 
3012  for (int i = 0; i < width; i++)
3013  {
3014  C(i,i) -= 1.0;
3015  }
3016  mfem::out << "size = " << width << ", i_max = " << C.MaxMaxNorm()
3017  << ", cond_F = " << FNorm()*copy.FNorm() << endl;
3018 }
3019 
3021 {
3022  if (capacity > 0)
3023  {
3024  delete [] data;
3025  }
3026 }
3027 
3028 
3029 
3030 void Add(const DenseMatrix &A, const DenseMatrix &B,
3031  double alpha, DenseMatrix &C)
3032 {
3033  for (int j = 0; j < C.Width(); j++)
3034  {
3035  for (int i = 0; i < C.Height(); i++)
3036  {
3037  C(i,j) = A(i,j) + alpha * B(i,j);
3038  }
3039  }
3040 }
3041 
3042 void Add(double alpha, const double *A,
3043  double beta, const double *B, DenseMatrix &C)
3044 {
3045  const int m = C.Height()*C.Width();
3046  double *C_data = C.GetData();
3047  for (int i = 0; i < m; i++)
3048  {
3049  C_data[i] = alpha*A[i] + beta*B[i];
3050  }
3051 }
3052 
3053 void Add(double alpha, const DenseMatrix &A,
3054  double beta, const DenseMatrix &B, DenseMatrix &C)
3055 {
3056  MFEM_ASSERT(A.Height() == C.Height(), "");
3057  MFEM_ASSERT(B.Height() == C.Height(), "");
3058  MFEM_ASSERT(A.Width() == C.Width(), "");
3059  MFEM_ASSERT(B.Width() == C.Width(), "");
3060  Add(alpha, A.GetData(), beta, B.GetData(), C);
3061 }
3062 
3063 
3064 void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
3065 {
3066  MFEM_ASSERT(a.Height() == b.Height() && a.Width() == c.Width() &&
3067  b.Width() == c.Height(), "incompatible dimensions");
3068 
3069 #ifdef MFEM_USE_LAPACK
3070  static char transa = 'N', transb = 'N';
3071  static double alpha = 1.0, beta = 0.0;
3072  int m = b.Height(), n = c.Width(), k = b.Width();
3073 
3074  dgemm_(&transa, &transb, &m, &n, &k, &alpha, b.Data(), &m,
3075  c.Data(), &k, &beta, a.Data(), &m);
3076 #else
3077  const int ah = a.Height();
3078  const int aw = a.Width();
3079  const int bw = b.Width();
3080  double *ad = a.Data();
3081  const double *bd = b.Data();
3082  const double *cd = c.Data();
3083  for (int i = 0; i < ah*aw; i++)
3084  {
3085  ad[i] = 0.0;
3086  }
3087  for (int j = 0; j < aw; j++)
3088  {
3089  for (int k = 0; k < bw; k++)
3090  {
3091  for (int i = 0; i < ah; i++)
3092  {
3093  ad[i+j*ah] += bd[i+k*ah] * cd[k+j*bw];
3094  }
3095  }
3096  }
3097 #endif
3098 }
3099 
3100 void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
3101 {
3102  MFEM_ASSERT(a.Height() == b.Height() && a.Width() == c.Width() &&
3103  b.Width() == c.Height(), "incompatible dimensions");
3104 
3105 #ifdef MFEM_USE_LAPACK
3106  static char transa = 'N', transb = 'N';
3107  static double alpha = 1.0, beta = 1.0;
3108  int m = b.Height(), n = c.Width(), k = b.Width();
3109 
3110  dgemm_(&transa, &transb, &m, &n, &k, &alpha, b.Data(), &m,
3111  c.Data(), &k, &beta, a.Data(), &m);
3112 #else
3113  const int ah = a.Height();
3114  const int aw = a.Width();
3115  const int bw = b.Width();
3116  double *ad = a.Data();
3117  const double *bd = b.Data();
3118  const double *cd = c.Data();
3119  for (int j = 0; j < aw; j++)
3120  {
3121  for (int k = 0; k < bw; k++)
3122  {
3123  for (int i = 0; i < ah; i++)
3124  {
3125  ad[i+j*ah] += bd[i+k*ah] * cd[k+j*bw];
3126  }
3127  }
3128  }
3129 #endif
3130 }
3131 
3132 void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
3133 {
3134 #ifdef MFEM_DEBUG
3135  if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
3136  {
3137  mfem_error("CalcAdjugate(...)");
3138  }
3139  if (a.Width() != adja.Height() || a.Height() != adja.Width())
3140  {
3141  mfem_error("CalcAdjugate(...)");
3142  }
3143 #endif
3144 
3145  if (a.Width() < a.Height())
3146  {
3147  const double *d = a.Data();
3148  double *ad = adja.Data();
3149  if (a.Width() == 1)
3150  {
3151  // N x 1, N = 2,3
3152  ad[0] = d[0];
3153  ad[1] = d[1];
3154  if (a.Height() == 3)
3155  {
3156  ad[2] = d[2];
3157  }
3158  }
3159  else
3160  {
3161  // 3 x 2
3162  double e, g, f;
3163  e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
3164  g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
3165  f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
3166 
3167  ad[0] = d[0]*g - d[3]*f;
3168  ad[1] = d[3]*e - d[0]*f;
3169  ad[2] = d[1]*g - d[4]*f;
3170  ad[3] = d[4]*e - d[1]*f;
3171  ad[4] = d[2]*g - d[5]*f;
3172  ad[5] = d[5]*e - d[2]*f;
3173  }
3174  return;
3175  }
3176 
3177  if (a.Width() == 1)
3178  {
3179  adja(0,0) = 1.0;
3180  }
3181  else if (a.Width() == 2)
3182  {
3183  adja(0,0) = a(1,1);
3184  adja(0,1) = -a(0,1);
3185  adja(1,0) = -a(1,0);
3186  adja(1,1) = a(0,0);
3187  }
3188  else
3189  {
3190  adja(0,0) = a(1,1)*a(2,2)-a(1,2)*a(2,1);
3191  adja(0,1) = a(0,2)*a(2,1)-a(0,1)*a(2,2);
3192  adja(0,2) = a(0,1)*a(1,2)-a(0,2)*a(1,1);
3193 
3194  adja(1,0) = a(1,2)*a(2,0)-a(1,0)*a(2,2);
3195  adja(1,1) = a(0,0)*a(2,2)-a(0,2)*a(2,0);
3196  adja(1,2) = a(0,2)*a(1,0)-a(0,0)*a(1,2);
3197 
3198  adja(2,0) = a(1,0)*a(2,1)-a(1,1)*a(2,0);
3199  adja(2,1) = a(0,1)*a(2,0)-a(0,0)*a(2,1);
3200  adja(2,2) = a(0,0)*a(1,1)-a(0,1)*a(1,0);
3201  }
3202 }
3203 
3205 {
3206 #ifdef MFEM_DEBUG
3207  if (a.Height() != a.Width() || adjat.Height() != adjat.Width() ||
3208  a.Width() != adjat.Width() || a.Width() < 1 || a.Width() > 3)
3209  {
3210  mfem_error("CalcAdjugateTranspose(...)");
3211  }
3212 #endif
3213  if (a.Width() == 1)
3214  {
3215  adjat(0,0) = 1.0;
3216  }
3217  else if (a.Width() == 2)
3218  {
3219  adjat(0,0) = a(1,1);
3220  adjat(1,0) = -a(0,1);
3221  adjat(0,1) = -a(1,0);
3222  adjat(1,1) = a(0,0);
3223  }
3224  else
3225  {
3226  adjat(0,0) = a(1,1)*a(2,2)-a(1,2)*a(2,1);
3227  adjat(1,0) = a(0,2)*a(2,1)-a(0,1)*a(2,2);
3228  adjat(2,0) = a(0,1)*a(1,2)-a(0,2)*a(1,1);
3229 
3230  adjat(0,1) = a(1,2)*a(2,0)-a(1,0)*a(2,2);
3231  adjat(1,1) = a(0,0)*a(2,2)-a(0,2)*a(2,0);
3232  adjat(2,1) = a(0,2)*a(1,0)-a(0,0)*a(1,2);
3233 
3234  adjat(0,2) = a(1,0)*a(2,1)-a(1,1)*a(2,0);
3235  adjat(1,2) = a(0,1)*a(2,0)-a(0,0)*a(2,1);
3236  adjat(2,2) = a(0,0)*a(1,1)-a(0,1)*a(1,0);
3237  }
3238 }
3239 
3240 void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
3241 {
3242  MFEM_ASSERT(a.Width() <= a.Height() && a.Width() >= 1 && a.Height() <= 3, "");
3243  MFEM_ASSERT(inva.Height() == a.Width(), "incorrect dimensions");
3244  MFEM_ASSERT(inva.Width() == a.Height(), "incorrect dimensions");
3245 
3246  double t;
3247 
3248  if (a.Width() < a.Height())
3249  {
3250  const double *d = a.Data();
3251  double *id = inva.Data();
3252  if (a.Height() == 2)
3253  {
3254  t = 1.0 / (d[0]*d[0] + d[1]*d[1]);
3255  id[0] = d[0] * t;
3256  id[1] = d[1] * t;
3257  }
3258  else
3259  {
3260  if (a.Width() == 1)
3261  {
3262  t = 1.0 / (d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
3263  id[0] = d[0] * t;
3264  id[1] = d[1] * t;
3265  id[2] = d[2] * t;
3266  }
3267  else
3268  {
3269  double e, g, f;
3270  e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
3271  g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
3272  f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
3273  t = 1.0 / (e*g - f*f);
3274  e *= t; g *= t; f *= t;
3275 
3276  id[0] = d[0]*g - d[3]*f;
3277  id[1] = d[3]*e - d[0]*f;
3278  id[2] = d[1]*g - d[4]*f;
3279  id[3] = d[4]*e - d[1]*f;
3280  id[4] = d[2]*g - d[5]*f;
3281  id[5] = d[5]*e - d[2]*f;
3282  }
3283  }
3284  return;
3285  }
3286 
3287 #ifdef MFEM_DEBUG
3288  t = a.Det();
3289  MFEM_ASSERT(std::abs(t) > 1.0e-14 * pow(a.FNorm()/a.Width(), a.Width()),
3290  "singular matrix!");
3291  t = 1.0 / t;
3292 #else
3293  t = 1.0 / a.Det();
3294 #endif
3295 
3296  switch (a.Height())
3297  {
3298  case 1:
3299  inva(0,0) = t;
3300  break;
3301  case 2:
3302  inva(0,0) = a(1,1) * t ;
3303  inva(0,1) = -a(0,1) * t ;
3304  inva(1,0) = -a(1,0) * t ;
3305  inva(1,1) = a(0,0) * t ;
3306  break;
3307  case 3:
3308  inva(0,0) = (a(1,1)*a(2,2)-a(1,2)*a(2,1))*t;
3309  inva(0,1) = (a(0,2)*a(2,1)-a(0,1)*a(2,2))*t;
3310  inva(0,2) = (a(0,1)*a(1,2)-a(0,2)*a(1,1))*t;
3311 
3312  inva(1,0) = (a(1,2)*a(2,0)-a(1,0)*a(2,2))*t;
3313  inva(1,1) = (a(0,0)*a(2,2)-a(0,2)*a(2,0))*t;
3314  inva(1,2) = (a(0,2)*a(1,0)-a(0,0)*a(1,2))*t;
3315 
3316  inva(2,0) = (a(1,0)*a(2,1)-a(1,1)*a(2,0))*t;
3317  inva(2,1) = (a(0,1)*a(2,0)-a(0,0)*a(2,1))*t;
3318  inva(2,2) = (a(0,0)*a(1,1)-a(0,1)*a(1,0))*t;
3319  break;
3320  }
3321 }
3322 
3324 {
3325 #ifdef MFEM_DEBUG
3326  if ( (a.Width() != a.Height()) || ( (a.Height()!= 1) && (a.Height()!= 2)
3327  && (a.Height()!= 3) ) )
3328  {
3329  mfem_error("CalcInverseTranspose(...)");
3330  }
3331 #endif
3332 
3333  double t = 1. / a.Det() ;
3334 
3335  switch (a.Height())
3336  {
3337  case 1:
3338  inva(0,0) = 1.0 / a(0,0);
3339  break;
3340  case 2:
3341  inva(0,0) = a(1,1) * t ;
3342  inva(1,0) = -a(0,1) * t ;
3343  inva(0,1) = -a(1,0) * t ;
3344  inva(1,1) = a(0,0) * t ;
3345  break;
3346  case 3:
3347  inva(0,0) = (a(1,1)*a(2,2)-a(1,2)*a(2,1))*t;
3348  inva(1,0) = (a(0,2)*a(2,1)-a(0,1)*a(2,2))*t;
3349  inva(2,0) = (a(0,1)*a(1,2)-a(0,2)*a(1,1))*t;
3350 
3351  inva(0,1) = (a(1,2)*a(2,0)-a(1,0)*a(2,2))*t;
3352  inva(1,1) = (a(0,0)*a(2,2)-a(0,2)*a(2,0))*t;
3353  inva(2,1) = (a(0,2)*a(1,0)-a(0,0)*a(1,2))*t;
3354 
3355  inva(0,2) = (a(1,0)*a(2,1)-a(1,1)*a(2,0))*t;
3356  inva(1,2) = (a(0,1)*a(2,0)-a(0,0)*a(2,1))*t;
3357  inva(2,2) = (a(0,0)*a(1,1)-a(0,1)*a(1,0))*t;
3358  break;
3359  }
3360 }
3361 
3362 void CalcOrtho(const DenseMatrix &J, Vector &n)
3363 {
3364  MFEM_ASSERT( ((J.Height() == 2 && J.Width() == 1)
3365  || (J.Height() == 3 && J.Width() == 2))
3366  && (J.Height() == n.Size()),
3367  "Matrix must be 3x2 or 2x1, "
3368  << "and the Vector must be sized with the rows. "
3369  << " J.Height() = " << J.Height()
3370  << ", J.Width() = " << J.Width()
3371  << ", n.Size() = " << n.Size()
3372  );
3373 
3374  const double *d = J.Data();
3375  if (J.Height() == 2)
3376  {
3377  n(0) = d[1];
3378  n(1) = -d[0];
3379  }
3380  else
3381  {
3382  n(0) = d[1]*d[5] - d[2]*d[4];
3383  n(1) = d[2]*d[3] - d[0]*d[5];
3384  n(2) = d[0]*d[4] - d[1]*d[3];
3385  }
3386 }
3387 
3388 void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
3389 {
3390  const int height = a.Height();
3391  const int width = a.Width();
3392  for (int i = 0; i < height; i++)
3393  {
3394  for (int j = 0; j <= i; j++)
3395  {
3396  double temp = 0.;
3397  for (int k = 0; k < width; k++)
3398  {
3399  temp += a(i,k) * a(j,k);
3400  }
3401  aat(j,i) = aat(i,j) = temp;
3402  }
3403  }
3404 }
3405 
3406 void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
3407 {
3408  for (int i = 0; i < A.Height(); i++)
3409  {
3410  for (int j = 0; j < i; j++)
3411  {
3412  double t = 0.;
3413  for (int k = 0; k < A.Width(); k++)
3414  {
3415  t += D(k) * A(i, k) * A(j, k);
3416  }
3417  ADAt(i, j) += t;
3418  ADAt(j, i) += t;
3419  }
3420  }
3421 
3422  // process diagonal
3423  for (int i = 0; i < A.Height(); i++)
3424  {
3425  double t = 0.;
3426  for (int k = 0; k < A.Width(); k++)
3427  {
3428  t += D(k) * A(i, k) * A(i, k);
3429  }
3430  ADAt(i, i) += t;
3431  }
3432 }
3433 
3434 void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
3435 {
3436  for (int i = 0; i < A.Height(); i++)
3437  {
3438  for (int j = 0; j <= i; j++)
3439  {
3440  double t = 0.;
3441  for (int k = 0; k < A.Width(); k++)
3442  {
3443  t += D(k) * A(i, k) * A(j, k);
3444  }
3445  ADAt(j, i) = ADAt(i, j) = t;
3446  }
3447  }
3448 }
3449 
3450 void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
3451 {
3452 #ifdef MFEM_DEBUG
3453  if (A.Height() != ABt.Height() || B.Height() != ABt.Width() ||
3454  A.Width() != B.Width())
3455  {
3456  mfem_error("MultABt(...)");
3457  }
3458 #endif
3459 
3460 #ifdef MFEM_USE_LAPACK
3461  static char transa = 'N', transb = 'T';
3462  static double alpha = 1.0, beta = 0.0;
3463  int m = A.Height(), n = B.Height(), k = A.Width();
3464 
3465  dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &m,
3466  B.Data(), &n, &beta, ABt.Data(), &m);
3467 #elif 1
3468  const int ah = A.Height();
3469  const int bh = B.Height();
3470  const int aw = A.Width();
3471  const double *ad = A.Data();
3472  const double *bd = B.Data();
3473  double *cd = ABt.Data();
3474 
3475  for (int i = 0, s = ah*bh; i < s; i++)
3476  {
3477  cd[i] = 0.0;
3478  }
3479  for (int k = 0; k < aw; k++)
3480  {
3481  double *cp = cd;
3482  for (int j = 0; j < bh; j++)
3483  {
3484  const double bjk = bd[j];
3485  for (int i = 0; i < ah; i++)
3486  {
3487  cp[i] += ad[i] * bjk;
3488  }
3489  cp += ah;
3490  }
3491  ad += ah;
3492  bd += bh;
3493  }
3494 #elif 1
3495  const int ah = A.Height();
3496  const int bh = B.Height();
3497  const int aw = A.Width();
3498  const double *ad = A.Data();
3499  const double *bd = B.Data();
3500  double *cd = ABt.Data();
3501 
3502  for (int j = 0; j < bh; j++)
3503  for (int i = 0; i < ah; i++)
3504  {
3505  double d = 0.0;
3506  const double *ap = ad + i;
3507  const double *bp = bd + j;
3508  for (int k = 0; k < aw; k++)
3509  {
3510  d += (*ap) * (*bp);
3511  ap += ah;
3512  bp += bh;
3513  }
3514  *(cd++) = d;
3515  }
3516 #else
3517  int i, j, k;
3518  double d;
3519 
3520  for (i = 0; i < A.Height(); i++)
3521  for (j = 0; j < B.Height(); j++)
3522  {
3523  d = 0.0;
3524  for (k = 0; k < A.Width(); k++)
3525  {
3526  d += A(i, k) * B(j, k);
3527  }
3528  ABt(i, j) = d;
3529  }
3530 #endif
3531 }
3532 
3533 void MultADBt(const DenseMatrix &A, const Vector &D,
3534  const DenseMatrix &B, DenseMatrix &ADBt)
3535 {
3536 #ifdef MFEM_DEBUG
3537  if (A.Height() != ADBt.Height() || B.Height() != ADBt.Width() ||
3538  A.Width() != B.Width() || A.Width() != D.Size())
3539  {
3540  mfem_error("MultADBt(...)");
3541  }
3542 #endif
3543 
3544  const int ah = A.Height();
3545  const int bh = B.Height();
3546  const int aw = A.Width();
3547  const double *ad = A.Data();
3548  const double *bd = B.Data();
3549  const double *dd = D.GetData();
3550  double *cd = ADBt.Data();
3551 
3552  for (int i = 0, s = ah*bh; i < s; i++)
3553  {
3554  cd[i] = 0.0;
3555  }
3556  for (int k = 0; k < aw; k++)
3557  {
3558  double *cp = cd;
3559  for (int j = 0; j < bh; j++)
3560  {
3561  const double dk_bjk = dd[k] * bd[j];
3562  for (int i = 0; i < ah; i++)
3563  {
3564  cp[i] += ad[i] * dk_bjk;
3565  }
3566  cp += ah;
3567  }
3568  ad += ah;
3569  bd += bh;
3570  }
3571 }
3572 
3573 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
3574 {
3575 #ifdef MFEM_DEBUG
3576  if (A.Height() != ABt.Height() || B.Height() != ABt.Width() ||
3577  A.Width() != B.Width())
3578  {
3579  mfem_error("AddMultABt(...)");
3580  }
3581 #endif
3582 
3583 #ifdef MFEM_USE_LAPACK
3584  static char transa = 'N', transb = 'T';
3585  static double alpha = 1.0, beta = 1.0;
3586  int m = A.Height(), n = B.Height(), k = A.Width();
3587 
3588  dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &m,
3589  B.Data(), &n, &beta, ABt.Data(), &m);
3590 #elif 1
3591  const int ah = A.Height();
3592  const int bh = B.Height();
3593  const int aw = A.Width();
3594  const double *ad = A.Data();
3595  const double *bd = B.Data();
3596  double *cd = ABt.Data();
3597 
3598  for (int k = 0; k < aw; k++)
3599  {
3600  double *cp = cd;
3601  for (int j = 0; j < bh; j++)
3602  {
3603  const double bjk = bd[j];
3604  for (int i = 0; i < ah; i++)
3605  {
3606  cp[i] += ad[i] * bjk;
3607  }
3608  cp += ah;
3609  }
3610  ad += ah;
3611  bd += bh;
3612  }
3613 #else
3614  int i, j, k;
3615  double d;
3616 
3617  for (i = 0; i < A.Height(); i++)
3618  for (j = 0; j < B.Height(); j++)
3619  {
3620  d = 0.0;
3621  for (k = 0; k < A.Width(); k++)
3622  {
3623  d += A(i, k) * B(j, k);
3624  }
3625  ABt(i, j) += d;
3626  }
3627 #endif
3628 }
3629 
3630 void AddMultADBt(const DenseMatrix &A, const Vector &D,
3631  const DenseMatrix &B, DenseMatrix &ADBt)
3632 {
3633 #ifdef MFEM_DEBUG
3634  if (A.Height() != ADBt.Height() || B.Height() != ADBt.Width() ||
3635  A.Width() != B.Width() || A.Width() != D.Size())
3636  {
3637  mfem_error("AddMultADBt(...)");
3638  }
3639 #endif
3640 
3641  const int ah = A.Height();
3642  const int bh = B.Height();
3643  const int aw = A.Width();
3644  const double *ad = A.Data();
3645  const double *bd = B.Data();
3646  const double *dd = D.GetData();
3647  double *cd = ADBt.Data();
3648 
3649  for (int k = 0; k < aw; k++)
3650  {
3651  double *cp = cd;
3652  for (int j = 0; j < bh; j++)
3653  {
3654  const double dk_bjk = dd[k] * bd[j];
3655  for (int i = 0; i < ah; i++)
3656  {
3657  cp[i] += ad[i] * dk_bjk;
3658  }
3659  cp += ah;
3660  }
3661  ad += ah;
3662  bd += bh;
3663  }
3664 }
3665 
3666 void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B,
3667  DenseMatrix &ABt)
3668 {
3669 #ifdef MFEM_DEBUG
3670  if (A.Height() != ABt.Height() || B.Height() != ABt.Width() ||
3671  A.Width() != B.Width())
3672  {
3673  mfem_error("AddMult_a_ABt(...)");
3674  }
3675 #endif
3676 
3677 #ifdef MFEM_USE_LAPACK
3678  static char transa = 'N', transb = 'T';
3679  double alpha = a;
3680  static double beta = 1.0;
3681  int m = A.Height(), n = B.Height(), k = A.Width();
3682 
3683  dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &m,
3684  B.Data(), &n, &beta, ABt.Data(), &m);
3685 #elif 1
3686  const int ah = A.Height();
3687  const int bh = B.Height();
3688  const int aw = A.Width();
3689  const double *ad = A.Data();
3690  const double *bd = B.Data();
3691  double *cd = ABt.Data();
3692 
3693  for (int k = 0; k < aw; k++)
3694  {
3695  double *cp = cd;
3696  for (int j = 0; j < bh; j++)
3697  {
3698  const double bjk = a * bd[j];
3699  for (int i = 0; i < ah; i++)
3700  {
3701  cp[i] += ad[i] * bjk;
3702  }
3703  cp += ah;
3704  }
3705  ad += ah;
3706  bd += bh;
3707  }
3708 #else
3709  int i, j, k;
3710  double d;
3711 
3712  for (i = 0; i < A.Height(); i++)
3713  for (j = 0; j < B.Height(); j++)
3714  {
3715  d = 0.0;
3716  for (k = 0; k < A.Width(); k++)
3717  {
3718  d += A(i, k) * B(j, k);
3719  }
3720  ABt(i, j) += a * d;
3721  }
3722 #endif
3723 }
3724 
3725 void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
3726 {
3727 #ifdef MFEM_DEBUG
3728  if (A.Width() != AtB.Height() || B.Width() != AtB.Width() ||
3729  A.Height() != B.Height())
3730  {
3731  mfem_error("MultAtB(...)");
3732  }
3733 #endif
3734 
3735 #ifdef MFEM_USE_LAPACK
3736  static char transa = 'T', transb = 'N';
3737  static double alpha = 1.0, beta = 0.0;
3738  int m = A.Width(), n = B.Width(), k = A.Height();
3739 
3740  dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &k,
3741  B.Data(), &k, &beta, AtB.Data(), &m);
3742 #elif 1
3743  const int ah = A.Height();
3744  const int aw = A.Width();
3745  const int bw = B.Width();
3746  const double *ad = A.Data();
3747  const double *bd = B.Data();
3748  double *cd = AtB.Data();
3749 
3750  for (int j = 0; j < bw; j++)
3751  {
3752  const double *ap = ad;
3753  for (int i = 0; i < aw; i++)
3754  {
3755  double d = 0.0;
3756  for (int k = 0; k < ah; k++)
3757  {
3758  d += ap[k] * bd[k];
3759  }
3760  *(cd++) = d;
3761  ap += ah;
3762  }
3763  bd += ah;
3764  }
3765 #else
3766  int i, j, k;
3767  double d;
3768 
3769  for (i = 0; i < A.Width(); i++)
3770  for (j = 0; j < B.Width(); j++)
3771  {
3772  d = 0.0;
3773  for (k = 0; k < A.Height(); k++)
3774  {
3775  d += A(k, i) * B(k, j);
3776  }
3777  AtB(i, j) = d;
3778  }
3779 #endif
3780 }
3781 
3782 void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
3783 {
3784  double d;
3785 
3786  for (int i = 0; i < A.Height(); i++)
3787  {
3788  for (int j = 0; j < i; j++)
3789  {
3790  d = 0.;
3791  for (int k = 0; k < A.Width(); k++)
3792  {
3793  d += A(i,k) * A(j,k);
3794  }
3795  AAt(i, j) += (d *= a);
3796  AAt(j, i) += d;
3797  }
3798  d = 0.;
3799  for (int k = 0; k < A.Width(); k++)
3800  {
3801  d += A(i,k) * A(i,k);
3802  }
3803  AAt(i, i) += a * d;
3804  }
3805 }
3806 
3807 void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
3808 {
3809  for (int i = 0; i < A.Height(); i++)
3810  {
3811  for (int j = 0; j <= i; j++)
3812  {
3813  double d = 0.;
3814  for (int k = 0; k < A.Width(); k++)
3815  {
3816  d += A(i,k) * A(j,k);
3817  }
3818  AAt(i, j) = AAt(j, i) = a * d;
3819  }
3820  }
3821 }
3822 
3823 void MultVVt(const Vector &v, DenseMatrix &vvt)
3824 {
3825  for (int i = 0; i < v.Size(); i++)
3826  {
3827  for (int j = 0; j <= i; j++)
3828  {
3829  vvt(i,j) = vvt(j,i) = v(i) * v(j);
3830  }
3831  }
3832 }
3833 
3834 void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
3835 {
3836 #ifdef MFEM_DEBUG
3837  if (v.Size() != VWt.Height() || w.Size() != VWt.Width())
3838  {
3839  mfem_error("MultVWt(...)");
3840  }
3841 #endif
3842 
3843  for (int i = 0; i < v.Size(); i++)
3844  {
3845  const double vi = v(i);
3846  for (int j = 0; j < w.Size(); j++)
3847  {
3848  VWt(i, j) = vi * w(j);
3849  }
3850  }
3851 }
3852 
3853 void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
3854 {
3855  const int m = v.Size(), n = w.Size();
3856 
3857 #ifdef MFEM_DEBUG
3858  if (VWt.Height() != m || VWt.Width() != n)
3859  {
3860  mfem_error("AddMultVWt(...)");
3861  }
3862 #endif
3863 
3864  for (int i = 0; i < m; i++)
3865  {
3866  const double vi = v(i);
3867  for (int j = 0; j < n; j++)
3868  {
3869  VWt(i, j) += vi * w(j);
3870  }
3871  }
3872 }
3873 
3874 void AddMultVVt(const Vector &v, DenseMatrix &VVt)
3875 {
3876  const int n = v.Size();
3877 
3878 #ifdef MFEM_DEBUG
3879  if (VVt.Height() != n || VVt.Width() != n)
3880  {
3881  mfem_error("AddMultVVt(...)");
3882  }
3883 #endif
3884 
3885  for (int i = 0; i < n; i++)
3886  {
3887  const double vi = v(i);
3888  for (int j = 0; j < i; j++)
3889  {
3890  const double vivj = vi * v(j);
3891  VVt(i, j) += vivj;
3892  VVt(j, i) += vivj;
3893  }
3894  VVt(i, i) += vi * vi;
3895  }
3896 }
3897 
3898 void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
3899  DenseMatrix &VWt)
3900 {
3901  const int m = v.Size(), n = w.Size();
3902 
3903 #ifdef MFEM_DEBUG
3904  if (VWt.Height() != m || VWt.Width() != n)
3905  {
3906  mfem_error("AddMult_a_VWt(...)");
3907  }
3908 #endif
3909 
3910  for (int j = 0; j < n; j++)
3911  {
3912  const double awj = a * w(j);
3913  for (int i = 0; i < m; i++)
3914  {
3915  VWt(i, j) += v(i) * awj;
3916  }
3917  }
3918 }
3919 
3920 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
3921 {
3922  MFEM_ASSERT(VVt.Height() == v.Size() && VVt.Width() == v.Size(),
3923  "incompatible dimensions!");
3924 
3925  const int n = v.Size();
3926  for (int i = 0; i < n; i++)
3927  {
3928  double avi = a * v(i);
3929  for (int j = 0; j < i; j++)
3930  {
3931  const double avivj = avi * v(j);
3932  VVt(i, j) += avivj;
3933  VVt(j, i) += avivj;
3934  }
3935  VVt(i, i) += avi * v(i);
3936  }
3937 }
3938 
3939 
3941 {
3942 #ifdef MFEM_USE_LAPACK
3943  int info = 0;
3944  if (m) { dgetrf_(&m, &m, data, &m, ipiv, &info); }
3945  MFEM_VERIFY(!info, "LAPACK: error in DGETRF");
3946 #else
3947  // compiling without LAPACK
3948  double *data = this->data;
3949  for (int i = 0; i < m; i++)
3950  {
3951  // pivoting
3952  {
3953  int piv = i;
3954  double a = std::abs(data[piv+i*m]);
3955  for (int j = i+1; j < m; j++)
3956  {
3957  const double b = std::abs(data[j+i*m]);
3958  if (b > a)
3959  {
3960  a = b;
3961  piv = j;
3962  }
3963  }
3964  ipiv[i] = piv;
3965  if (piv != i)
3966  {
3967  // swap rows i and piv in both L and U parts
3968  for (int j = 0; j < m; j++)
3969  {
3970  Swap<double>(data[i+j*m], data[piv+j*m]);
3971  }
3972  }
3973  }
3974  MFEM_ASSERT(data[i+i*m] != 0.0, "division by zero");
3975  const double a_ii_inv = 1.0/data[i+i*m];
3976  for (int j = i+1; j < m; j++)
3977  {
3978  data[j+i*m] *= a_ii_inv;
3979  }
3980  for (int k = i+1; k < m; k++)
3981  {
3982  const double a_ik = data[i+k*m];
3983  for (int j = i+1; j < m; j++)
3984  {
3985  data[j+k*m] -= a_ik * data[j+i*m];
3986  }
3987  }
3988  }
3989 #endif
3990 }
3991 
3992 double LUFactors::Det(int m) const
3993 {
3994  double det = 1.0;
3995  for (int i=0; i<m; i++)
3996  {
3997  if (ipiv[i] != i-ipiv_base)
3998  {
3999  det *= -data[m * i + i];
4000  }
4001  else
4002  {
4003  det *= data[m * i + i];
4004  }
4005  }
4006  return det;
4007 }
4008 
4009 void LUFactors::Mult(int m, int n, double *X) const
4010 {
4011  const double *data = this->data;
4012  const int *ipiv = this->ipiv;
4013  double *x = X;
4014  for (int k = 0; k < n; k++)
4015  {
4016  // X <- U X
4017  for (int i = 0; i < m; i++)
4018  {
4019  double x_i = x[i] * data[i+i*m];
4020  for (int j = i+1; j < m; j++)
4021  {
4022  x_i += x[j] * data[i+j*m];
4023  }
4024  x[i] = x_i;
4025  }
4026  // X <- L X
4027  for (int i = m-1; i >= 0; i--)
4028  {
4029  double x_i = x[i];
4030  for (int j = 0; j < i; j++)
4031  {
4032  x_i += x[j] * data[i+j*m];
4033  }
4034  x[i] = x_i;
4035  }
4036  // X <- P^{-1} X
4037  for (int i = m-1; i >= 0; i--)
4038  {
4039  Swap<double>(x[i], x[ipiv[i]-ipiv_base]);
4040  }
4041  x += m;
4042  }
4043 }
4044 
4045 void LUFactors::LSolve(int m, int n, double *X) const
4046 {
4047  const double *data = this->data;
4048  const int *ipiv = this->ipiv;
4049  double *x = X;
4050  for (int k = 0; k < n; k++)
4051  {
4052  // X <- P X
4053  for (int i = 0; i < m; i++)
4054  {
4055  Swap<double>(x[i], x[ipiv[i]-ipiv_base]);
4056  }
4057  // X <- L^{-1} X
4058  for (int j = 0; j < m; j++)
4059  {
4060  const double x_j = x[j];
4061  for (int i = j+1; i < m; i++)
4062  {
4063  x[i] -= data[i+j*m] * x_j;
4064  }
4065  }
4066  x += m;
4067  }
4068 }
4069 
4070 void LUFactors::USolve(int m, int n, double *X) const
4071 {
4072  const double *data = this->data;
4073  double *x = X;
4074  // X <- U^{-1} X
4075  for (int k = 0; k < n; k++)
4076  {
4077  for (int j = m-1; j >= 0; j--)
4078  {
4079  const double x_j = ( x[j] /= data[j+j*m] );
4080  for (int i = 0; i < j; i++)
4081  {
4082  x[i] -= data[i+j*m] * x_j;
4083  }
4084  }
4085  x += m;
4086  }
4087 }
4088 
4089 void LUFactors::Solve(int m, int n, double *X) const
4090 {
4091 #ifdef MFEM_USE_LAPACK
4092  char trans = 'N';
4093  int info = 0;
4094  if (m > 0 && n > 0) { dgetrs_(&trans, &m, &n, data, &m, ipiv, X, &m, &info); }
4095  MFEM_VERIFY(!info, "LAPACK: error in DGETRS");
4096 #else
4097  // compiling without LAPACK
4098  LSolve(m, n, X);
4099  USolve(m, n, X);
4100 #endif
4101 }
4102 
4103 void LUFactors::GetInverseMatrix(int m, double *X) const
4104 {
4105  // A^{-1} = U^{-1} L^{-1} P
4106  const double *data = this->data;
4107  const int *ipiv = this->ipiv;
4108  // X <- U^{-1} (set only the upper triangular part of X)
4109  double *x = X;
4110  for (int k = 0; k < m; k++)
4111  {
4112  const double minus_x_k = -( x[k] = 1.0/data[k+k*m] );
4113  for (int i = 0; i < k; i++)
4114  {
4115  x[i] = data[i+k*m] * minus_x_k;
4116  }
4117  for (int j = k-1; j >= 0; j--)
4118  {
4119  const double x_j = ( x[j] /= data[j+j*m] );
4120  for (int i = 0; i < j; i++)
4121  {
4122  x[i] -= data[i+j*m] * x_j;
4123  }
4124  }
4125  x += m;
4126  }
4127  // X <- X L^{-1} (use input only from the upper triangular part of X)
4128  {
4129  int k = m-1;
4130  for (int j = 0; j < k; j++)
4131  {
4132  const double minus_L_kj = -data[k+j*m];
4133  for (int i = 0; i <= j; i++)
4134  {
4135  X[i+j*m] += X[i+k*m] * minus_L_kj;
4136  }
4137  for (int i = j+1; i < m; i++)
4138  {
4139  X[i+j*m] = X[i+k*m] * minus_L_kj;
4140  }
4141  }
4142  }
4143  for (int k = m-2; k >= 0; k--)
4144  {
4145  for (int j = 0; j < k; j++)
4146  {
4147  const double L_kj = data[k+j*m];
4148  for (int i = 0; i < m; i++)
4149  {
4150  X[i+j*m] -= X[i+k*m] * L_kj;
4151  }
4152  }
4153  }
4154  // X <- X P
4155  for (int k = m-1; k >= 0; k--)
4156  {
4157  const int piv_k = ipiv[k]-ipiv_base;
4158  if (k != piv_k)
4159  {
4160  for (int i = 0; i < m; i++)
4161  {
4162  Swap<double>(X[i+k*m], X[i+piv_k*m]);
4163  }
4164  }
4165  }
4166 }
4167 
4168 void LUFactors::SubMult(int m, int n, int r, const double *A21,
4169  const double *X1, double *X2)
4170 {
4171  // X2 <- X2 - A21 X1
4172  for (int k = 0; k < r; k++)
4173  {
4174  for (int j = 0; j < m; j++)
4175  {
4176  const double x1_jk = X1[j+k*m];
4177  for (int i = 0; i < n; i++)
4178  {
4179  X2[i+k*n] -= A21[i+j*n] * x1_jk;
4180  }
4181  }
4182  }
4183 }
4184 
4186  int m, int n, double *A12, double *A21, double *A22) const
4187 {
4188  const double *data = this->data;
4189  // A12 <- L^{-1} P A12
4190  LSolve(m, n, A12);
4191  // A21 <- A21 U^{-1}
4192  for (int j = 0; j < m; j++)
4193  {
4194  const double u_jj_inv = 1.0/data[j+j*m];
4195  for (int i = 0; i < n; i++)
4196  {
4197  A21[i+j*n] *= u_jj_inv;
4198  }
4199  for (int k = j+1; k < m; k++)
4200  {
4201  const double u_jk = data[j+k*m];
4202  for (int i = 0; i < n; i++)
4203  {
4204  A21[i+k*n] -= A21[i+j*n] * u_jk;
4205  }
4206  }
4207  }
4208  // A22 <- A22 - A21 A12
4209  SubMult(m, n, n, A21, A12, A22);
4210 }
4211 
4212 void LUFactors::BlockForwSolve(int m, int n, int r, const double *L21,
4213  double *B1, double *B2) const
4214 {
4215  // B1 <- L^{-1} P B1
4216  LSolve(m, r, B1);
4217  // B2 <- B2 - L21 B1
4218  SubMult(m, n, r, L21, B1, B2);
4219 }
4220 
4221 void LUFactors::BlockBackSolve(int m, int n, int r, const double *U12,
4222  const double *X2, double *Y1) const
4223 {
4224  // Y1 <- Y1 - U12 X2
4225  SubMult(n, m, r, U12, X2, Y1);
4226  // Y1 <- U^{-1} Y1
4227  USolve(m, r, Y1);
4228 }
4229 
4230 
4232  : MatrixInverse(mat)
4233 {
4234  MFEM_ASSERT(height == width, "not a square matrix");
4235  a = &mat;
4236  lu.data = new double[width*width];
4237  lu.ipiv = new int[width];
4238  Factor();
4239 }
4240 
4242  : MatrixInverse(*mat)
4243 {
4244  MFEM_ASSERT(height == width, "not a square matrix");
4245  a = mat;
4246  lu.data = new double[width*width];
4247  lu.ipiv = new int[width];
4248 }
4249 
4251 {
4252  MFEM_ASSERT(a, "DenseMatrix is not given");
4253  const double *adata = a->data;
4254  const int s = width*width;
4255  for (int i = 0; i < s; i++)
4256  {
4257  lu.data[i] = adata[i];
4258  }
4259  lu.Factor(width);
4260 }
4261 
4263 {
4264  Ainv.SetSize(width);
4265  lu.GetInverseMatrix(width, Ainv.Data());
4266 }
4267 
4269 {
4270  MFEM_VERIFY(mat.height == mat.width, "DenseMatrix is not square!");
4271  if (width != mat.width)
4272  {
4273  height = width = mat.width;
4274  delete [] lu.data;
4275  lu.data = new double[width*width];
4276  delete [] lu.ipiv;
4277  lu.ipiv = new int[width];
4278  }
4279  a = &mat;
4280  Factor();
4281 }
4282 
4284 {
4285  const DenseMatrix *p = dynamic_cast<const DenseMatrix*>(&op);
4286  MFEM_VERIFY(p != NULL, "Operator is not a DenseMatrix!");
4287  Factor(*p);
4288 }
4289 
4290 void DenseMatrixInverse::Mult(const Vector &x, Vector &y) const
4291 {
4292  y = x;
4293  lu.Solve(width, 1, y.GetData());
4294 }
4295 
4297 {
4298  X = B;
4299  lu.Solve(width, X.Width(), X.Data());
4300 }
4301 
4303 {
4304  DenseMatrix C(width);
4305  Mult(*a, C);
4306  for (int i = 0; i < width; i++)
4307  {
4308  C(i,i) -= 1.0;
4309  }
4310  mfem::out << "size = " << width << ", i_max = " << C.MaxMaxNorm() << endl;
4311 }
4312 
4314 {
4315  delete [] lu.data;
4316  delete [] lu.ipiv;
4317 }
4318 
4319 
4321  : mat(m)
4322 {
4323  n = mat.Width();
4324  EVal.SetSize(n);
4325  EVect.SetSize(n);
4326  ev.SetDataAndSize(NULL, n);
4327 
4328 #ifdef MFEM_USE_LAPACK
4329  jobz = 'V';
4330  uplo = 'U';
4331  lwork = -1;
4332  double qwork;
4333  dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
4334  &qwork, &lwork, &info);
4335 
4336  lwork = (int) qwork;
4337  work = new double[lwork];
4338 #endif
4339 }
4340 
4342  const DenseMatrixEigensystem &other)
4343  : mat(other.mat), EVal(other.EVal), EVect(other.EVect), ev(NULL, other.n),
4344  n(other.n)
4345 {
4346 #ifdef MFEM_USE_LAPACK
4347  jobz = other.jobz;
4348  uplo = other.uplo;
4349  lwork = other.lwork;
4350 
4351  work = new double[lwork];
4352 #endif
4353 }
4354 
4356 {
4357 #ifdef MFEM_DEBUG
4358  if (mat.Width() != n)
4359  {
4360  mfem_error("DenseMatrixEigensystem::Eval()");
4361  }
4362 #endif
4363 
4364 #ifdef MFEM_USE_LAPACK
4365  EVect = mat;
4366  dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
4367  work, &lwork, &info);
4368 
4369  if (info != 0)
4370  {
4371  mfem::err << "DenseMatrixEigensystem::Eval(): DSYEV error code: "
4372  << info << endl;
4373  mfem_error();
4374  }
4375 #else
4376  mfem_error("DenseMatrixEigensystem::Eval(): Compiled without LAPACK");
4377 #endif
4378 }
4379 
4381 {
4382 #ifdef MFEM_USE_LAPACK
4383  delete [] work;
4384 #endif
4385 }
4386 
4387 
4389 {
4390  m = M.Height();
4391  n = M.Width();
4392  Init();
4393 }
4394 
4396 {
4397  m = h;
4398  n = w;
4399  Init();
4400 }
4401 
4402 void DenseMatrixSVD::Init()
4403 {
4404 #ifdef MFEM_USE_LAPACK
4405  sv.SetSize(min(m, n));
4406 
4407  jobu = 'N';
4408  jobvt = 'N';
4409 
4410  double qwork;
4411  lwork = -1;
4412  dgesvd_(&jobu, &jobvt, &m, &n, NULL, &m, sv.GetData(), NULL, &m,
4413  NULL, &n, &qwork, &lwork, &info);
4414 
4415  lwork = (int) qwork;
4416  work = new double[lwork];
4417 #else
4418  mfem_error("DenseMatrixSVD::Init(): Compiled without LAPACK");
4419 #endif
4420 }
4421 
4423 {
4424 #ifdef MFEM_DEBUG
4425  if (M.Height() != m || M.Width() != n)
4426  {
4427  mfem_error("DenseMatrixSVD::Eval()");
4428  }
4429 #endif
4430 
4431 #ifdef MFEM_USE_LAPACK
4432  dgesvd_(&jobu, &jobvt, &m, &n, M.Data(), &m, sv.GetData(), NULL, &m,
4433  NULL, &n, work, &lwork, &info);
4434 
4435  if (info)
4436  {
4437  mfem::err << "DenseMatrixSVD::Eval() : info = " << info << endl;
4438  mfem_error();
4439  }
4440 #else
4441  mfem_error("DenseMatrixSVD::Eval(): Compiled without LAPACK");
4442 #endif
4443 }
4444 
4446 {
4447 #ifdef MFEM_USE_LAPACK
4448  delete [] work;
4449 #endif
4450 }
4451 
4452 
4453 void DenseTensor::AddMult(const Table &elem_dof, const Vector &x, Vector &y)
4454 const
4455 {
4456  int n = SizeI(), ne = SizeK();
4457  const int *I = elem_dof.GetI(), *J = elem_dof.GetJ(), *dofs;
4458  const double *d_col = tdata;
4459  double *yp = y, x_col;
4460  const double *xp = x;
4461  // the '4' here can be tuned for given platform and compiler
4462  if (n <= 4)
4463  {
4464  for (int i = 0; i < ne; i++)
4465  {
4466  dofs = J + I[i];
4467  for (int col = 0; col < n; col++)
4468  {
4469  x_col = xp[dofs[col]];
4470  for (int row = 0; row < n; row++)
4471  {
4472  yp[dofs[row]] += x_col*d_col[row];
4473  }
4474  d_col += n;
4475  }
4476  }
4477  }
4478  else
4479  {
4480  Vector ye(n);
4481  for (int i = 0; i < ne; i++)
4482  {
4483  dofs = J + I[i];
4484  x_col = xp[dofs[0]];
4485  for (int row = 0; row < n; row++)
4486  {
4487  ye(row) = x_col*d_col[row];
4488  }
4489  d_col += n;
4490  for (int col = 1; col < n; col++)
4491  {
4492  x_col = xp[dofs[col]];
4493  for (int row = 0; row < n; row++)
4494  {
4495  ye(row) += x_col*d_col[row];
4496  }
4497  d_col += n;
4498  }
4499  for (int row = 0; row < n; row++)
4500  {
4501  yp[dofs[row]] += ye(row);
4502  }
4503  }
4504  }
4505 }
4506 
4508 {
4509  int s = SizeI() * SizeJ() * SizeK();
4510  for (int i=0; i<s; i++)
4511  {
4512  tdata[i] = c;
4513  }
4514  return *this;
4515 }
4516 
4517 }
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: densemat.cpp:2518
void dsyevr_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
Definition: densemat.cpp:875
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
Definition: densemat.cpp:3450
int Size() const
Logical size of the array.
Definition: array.hpp:118
DenseMatrix & operator-=(const DenseMatrix &m)
Definition: densemat.cpp:626
void SymmetricScaling(const Vector &s)
SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
Definition: densemat.cpp:393
void SquareRootInverse()
Replaces the current matrix with its square root inverse.
Definition: densemat.cpp:768
int CheckFinite(const double *v, const int n)
Definition: vector.hpp:376
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3853
int * GetJ()
Definition: table.hpp:114
DenseMatrix & operator*=(double c)
Definition: densemat.cpp:639
void GetDiag(Vector &d) const
Returns the diagonal of the matrix.
Definition: densemat.cpp:2410
bool KernelVector2G(const int &mode, double &d1, double &d12, double &d21, double &d2)
Definition: densemat.cpp:1325
DenseTensor & operator=(double c)
Sets the tensor elements equal to constant c.
Definition: densemat.cpp:4507
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:3834
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition: densemat.cpp:379
void Eigenvalues2S(const double &d12, double &d1, double &d2)
Definition: densemat.cpp:1232
void GetScalingFactor(const double &d_max, double &mult)
Definition: densemat.cpp:1822
void SingularValues(Vector &sv) const
Definition: densemat.cpp:1176
void dsyev_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
Definition: densemat.cpp:1035
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:400
double Det() const
Definition: densemat.cpp:472
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
int SizeK() const
Definition: densemat.hpp:695
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
Definition: densemat.cpp:4185
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:4221
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:320
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:3132
void dgetri_(int *N, double *A, int *LDA, int *IPIV, double *WORK, int *LWORK, int *INFO)
void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const
Definition: densemat.cpp:4453
void TestInversion()
Invert and print the numerical conditioning of the inversion.
Definition: densemat.cpp:3006
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:150
void CopyRows(const DenseMatrix &A, int row1, int row2)
Copy rows row1 through row2 from A to *this.
Definition: densemat.cpp:2623
double Det(int m) const
Definition: densemat.cpp:3992
void Eval(DenseMatrix &M)
Definition: densemat.cpp:4422
Abstract data type for matrix inverse.
Definition: matrix.hpp:66
void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
Definition: densemat.cpp:3666
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:4262
void Factor(int m)
Definition: densemat.cpp:3940
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:4250
void GetInverseMatrix(int m, double *X) const
Assuming L.U = P.A factored data of size (m x m), compute X &lt;- A^{-1}.
Definition: densemat.cpp:4103
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:96
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:159
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:3362
DenseMatrix & operator=(double c)
Sets the matrix elements equal to constant c.
Definition: densemat.cpp:576
void vec_normalize3_aux(const double &x1, const double &x2, const double &x3, double &n1, double &n2, double &n3)
Definition: densemat.cpp:1280
void Set(double alpha, const double *A)
Set the matrix to alpha * A, assuming that A has the same dimensions as the matrix and uses column-ma...
Definition: densemat.cpp:556
void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt = a * A * A^t.
Definition: densemat.cpp:3807
void dgesvd_(char *JOBU, char *JOBVT, int *M, int *N, double *A, int *LDA, double *S, double *U, int *LDU, double *VT, int *LDVT, double *WORK, int *LWORK, int *INFO)
void dsygv_Eigensystem(DenseMatrix &a, DenseMatrix &b, Vector &ev, DenseMatrix *evect)
Definition: densemat.cpp:1103
static void SubMult(int m, int n, int r, const double *A21, const double *X1, double *X2)
Definition: densemat.cpp:4168
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: densemat.cpp:2935
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:4290
int KernelVector3G_aux(const int &mode, double &d1, double &d2, double &d3, double &c12, double &c13, double &c23, double &c21, double &c31, double &c32)
Definition: densemat.cpp:1466
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:3030
double & operator()(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.hpp:752
double Weight() const
Definition: densemat.cpp:529
void USolve(int m, int n, double *X) const
Definition: densemat.cpp:4070
double FNorm() const
Compute the Frobenius norm of the matrix.
Definition: densemat.hpp:216
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:224
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:3204
friend class DenseMatrixInverse
Definition: densemat.hpp:26
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:3100
double operator*(const DenseMatrix &m) const
Matrix inner product: tr(A^t B)
Definition: densemat.cpp:209
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:565
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:146
void AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
Definition: densemat.cpp:3898
void InvSymmetricScaling(const Vector &s)
InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
Definition: densemat.cpp:421
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:2380
void BlockForwSolve(int m, int n, int r, const double *L21, double *B1, double *B2) const
Definition: densemat.cpp:4212
DenseMatrixSVD(DenseMatrix &M)
Definition: densemat.cpp:4388
Abstract data type matrix.
Definition: matrix.hpp:27
void Norm2(double *v) const
Take the 2-norm of the columns of A and store in v.
Definition: densemat.cpp:813
void MultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
Definition: densemat.cpp:3533
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:658
virtual ~DenseMatrixInverse()
Destroys dense inverse matrix.
Definition: densemat.cpp:4313
void LSolve(int m, int n, double *X) const
Definition: densemat.cpp:4045
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
Definition: densemat.cpp:338
double Det() const
Compute the determinant of the original DenseMatrix using the LU factors.
Definition: densemat.hpp:590
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
Definition: densemat.cpp:3874
void CopyMNDiag(double c, int n, int row_offset, int col_offset)
Copy c on the diagonal of size n to *this at row_offset, col_offset.
Definition: densemat.cpp:2719
void vec_normalize3(const double &x1, const double &x2, const double &x3, double &n1, double &n2, double &n3)
Definition: densemat.cpp:1297
void dgetrf_(int *, int *, double *, int *, int *, int *)
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
Definition: densemat.cpp:3920
void Neg()
(*this) = -(*this)
Definition: densemat.cpp:649
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: densemat.cpp:4283
void Solve(int m, int n, double *X) const
Definition: densemat.cpp:4089
int KernelVector3S(const int &mode, const double &d12, const double &d13, const double &d23, double &d1, double &d2, double &d3)
Definition: densemat.cpp:1571
void SetRow(int r, const Vector &row)
Definition: densemat.cpp:2905
void Getl1Diag(Vector &l) const
Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|.
Definition: densemat.cpp:2424
void AddToVector(int offset, Vector &v) const
Add the matrix &#39;data&#39; to the Vector &#39;v&#39; at the given &#39;offset&#39;.
Definition: densemat.cpp:2840
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:2396
void AddMult(const Vector &x, Vector &y) const
y += A.x
Definition: densemat.cpp:247
void trans(const Vector &x, Vector &p)
Definition: toroid.cpp:239
void dgemm_(char *, char *, int *, int *, int *, double *, double *, int *, double *, int *, double *, double *, int *)
void Threshold(double eps)
Replace small entries, abs(a_ij) &lt;= eps, with zero.
Definition: densemat.cpp:2921
int SizeI() const
Definition: densemat.hpp:693
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:3240
void TestInversion()
Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
Definition: densemat.cpp:4302
void Eigensystem2S(const double &d12, double &d1, double &d2, double &c, double &s)
Definition: densemat.cpp:1251
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:826
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:94
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:575
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:2485
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:3823
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:448
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:3573
void dgetrs_(char *, int *, int *, double *, int *, int *, double *, int *, int *)
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:69
int Reduce3S(const int &mode, double &d1, double &d2, double &d3, double &d12, double &d13, double &d23, double &z1, double &z2, double &z3, double &v1, double &v2, double &v3, double &g)
Definition: densemat.cpp:1700
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:3388
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0&lt;=i&lt;A.Height, 0&lt;=j&lt;A.Width.
Definition: densemat.cpp:2780
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:3323
A class to initialize the size of a Tensor.
Definition: dtensor.hpp:55
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:125
void dsygv_(int *ITYPE, char *JOBZ, char *UPLO, int *N, double *A, int *LDA, double *B, int *LDB, double *W, double *WORK, int *LWORK, int *INFO)
void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt = A D A^t, where D is diagonal.
Definition: densemat.cpp:3434
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
virtual void PrintT(std::ostream &out=mfem::out, int width_=4) const
Prints the transpose matrix to stream out.
Definition: densemat.cpp:2980
void CopyCols(const DenseMatrix &A, int col1, int col2)
Copy columns col1 through col2 from A to *this.
Definition: densemat.cpp:2636
int SizeJ() const
Definition: densemat.hpp:694
virtual MatrixInverse * Inverse() const
Returns a pointer to the inverse matrix.
Definition: densemat.cpp:467
void AddMultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
Definition: densemat.cpp:3630
virtual ~DenseMatrix()
Destroys dense matrix.
Definition: densemat.cpp:3020
void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset)
Copy matrix A^t to the location in *this at row_offset, col_offset.
Definition: densemat.cpp:2675
void AddMultTranspose(const Vector &x, Vector &y) const
y += A^t x
Definition: densemat.cpp:265
void CopyExceptMN(const DenseMatrix &A, int m, int n)
Copy All rows and columns except m and n from A.
Definition: densemat.cpp:2754
void Diag(double c, int n)
Creates n x n diagonal matrix with diagonal elements c.
Definition: densemat.cpp:2455
void Mult(int m, int n, double *X) const
Definition: densemat.cpp:4009
static const int ipiv_base
Definition: densemat.hpp:462
void GradToCurl(DenseMatrix &curl)
Definition: densemat.cpp:2549
DenseMatrixInverse()
Default constructor.
Definition: densemat.hpp:557
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1842
void GetRowSums(Vector &l) const
Compute the row sums of the DenseMatrix.
Definition: densemat.cpp:2441
void CalcEigenvalues(double *lambda, double *vec) const
Definition: densemat.cpp:2140
DenseMatrixEigensystem(DenseMatrix &m)
Definition: densemat.cpp:4320
double epsilon(const Vector &x)
Definition: maxwell.cpp:90
void dsyevr_(char *JOBZ, char *RANGE, char *UPLO, int *N, double *A, int *LDA, double *VL, double *VU, int *IL, int *IU, double *ABSTOL, int *M, double *W, double *Z, int *LDZ, int *ISUPPZ, double *WORK, int *LWORK, int *IWORK, int *LIWORK, int *INFO)
int Rank(double tol) const
Definition: densemat.cpp:1215
const double alpha
Definition: ex15.cpp:339
void AddMult_a(double a, const Vector &x, Vector &y) const
y += a * A.x
Definition: densemat.cpp:283
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
Definition: densemat.cpp:364
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:3725
Vector data type.
Definition: vector.hpp:48
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
void AddMultTranspose_a(double a, const Vector &x, Vector &y) const
y += a * A^t x
Definition: densemat.cpp:301
void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt += A D A^t, where D is diagonal.
Definition: densemat.cpp:3406
void GetFromVector(int offset, const Vector &v)
Get the matrix &#39;data&#39; from the Vector &#39;v&#39; at the given &#39;offset&#39;.
Definition: densemat.cpp:2851
int * GetI()
Definition: table.hpp:113
void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this.
Definition: densemat.cpp:2649
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
Definition: densemat.cpp:351
void SetCol(int c, const Vector &col)
Definition: densemat.cpp:2913
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:88
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
Abstract operator.
Definition: operator.hpp:21
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:656
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.cpp:163
void AdjustDofDirection(Array< int > &dofs)
Definition: densemat.cpp:2862
void GradToDiv(Vector &div)
Definition: densemat.cpp:2608
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:3782
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
void dsyev_(char *JOBZ, char *UPLO, int *N, double *A, int *LDA, double *W, double *WORK, int *LWORK, int *INFO)
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Definition: densemat.cpp:2961
double * data
Definition: densemat.hpp:459
DenseMatrix & operator+=(const double *m)
Definition: densemat.cpp:609