MFEM  v4.4.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-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 
13 // Implementation of data types dense matrix, inverse dense matrix
14 
15 
16 #include "kernels.hpp"
17 #include "vector.hpp"
18 #include "matrix.hpp"
19 #include "densemat.hpp"
20 #include "kernels.hpp"
21 #include "../general/forall.hpp"
22 #include "../general/table.hpp"
23 #include "../general/globals.hpp"
24 
25 #include <iostream>
26 #include <iomanip>
27 #include <limits>
28 #include <algorithm>
29 #include <cstdlib>
30 #if defined(_MSC_VER) && (_MSC_VER < 1800)
31 #include <float.h>
32 #define copysign _copysign
33 #endif
34 
35 
36 #ifdef MFEM_USE_LAPACK
37 extern "C" void
38 dgemm_(char *, char *, int *, int *, int *, double *, double *,
39  int *, double *, int *, double *, double *, int *);
40 extern "C" void
41 dgetrf_(int *, int *, double *, int *, int *, int *);
42 extern "C" void
43 dgetrs_(char *, int *, int *, double *, int *, int *, double *, int *, int *);
44 extern "C" void
45 dgetri_(int *N, double *A, int *LDA, int *IPIV, double *WORK,
46  int *LWORK, int *INFO);
47 extern "C" void
48 dsyevr_(char *JOBZ, char *RANGE, char *UPLO, int *N, double *A, int *LDA,
49  double *VL, double *VU, int *IL, int *IU, double *ABSTOL, int *M,
50  double *W, double *Z, int *LDZ, int *ISUPPZ, double *WORK, int *LWORK,
51  int *IWORK, int *LIWORK, int *INFO);
52 extern "C" void
53 dsyev_(char *JOBZ, char *UPLO, int *N, double *A, int *LDA, double *W,
54  double *WORK, int *LWORK, int *INFO);
55 extern "C" void
56 dsygv_ (int *ITYPE, char *JOBZ, char *UPLO, int * N, double *A, int *LDA,
57  double *B, int *LDB, double *W, double *WORK, int *LWORK, int *INFO);
58 extern "C" void
59 dgesvd_(char *JOBU, char *JOBVT, int *M, int *N, double *A, int *LDA,
60  double *S, double *U, int *LDU, double *VT, int *LDVT, double *WORK,
61  int *LWORK, int *INFO);
62 extern "C" void
63 dtrsm_(char *side, char *uplo, char *transa, char *diag, int *m, int *n,
64  double *alpha, double *a, int *lda, double *b, int *ldb);
65 extern "C" void
66 dggev_(char *jobvl, char *jobvr, int *n, double *a, int *lda, double *B,
67  int *ldb, double *alphar, double *alphai, double *beta, double *vl,
68  int * ldvl, double * vr, int * ldvr, double * work, int * lwork, int* info);
69 #endif
70 
71 
72 namespace mfem
73 {
74 
75 using namespace std;
76 
78 
79 DenseMatrix::DenseMatrix(const DenseMatrix &m) : Matrix(m.height, m.width)
80 {
81  const int hw = height * width;
82  if (hw > 0)
83  {
84  MFEM_ASSERT(m.data, "invalid source matrix");
85  data.New(hw);
86  std::memcpy(data, m.data, sizeof(double)*hw);
87  }
88 }
89 
91 {
92  MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
93  if (s > 0)
94  {
95  data.New(s*s);
96  *this = 0.0; // init with zeroes
97  }
98 }
99 
100 DenseMatrix::DenseMatrix(int m, int n) : Matrix(m, n)
101 {
102  MFEM_ASSERT(m >= 0 && n >= 0,
103  "invalid DenseMatrix size: " << m << " x " << n);
104  const int capacity = m*n;
105  if (capacity > 0)
106  {
107  data.New(capacity);
108  *this = 0.0; // init with zeroes
109  }
110 }
111 
113  : Matrix(mat.width, mat.height)
114 {
115  MFEM_CONTRACT_VAR(ch);
116  const int capacity = height*width;
117  if (capacity > 0)
118  {
119  data.New(capacity);
120 
121  for (int i = 0; i < height; i++)
122  {
123  for (int j = 0; j < width; j++)
124  {
125  (*this)(i,j) = mat(j,i);
126  }
127  }
128  }
129 }
130 
131 void DenseMatrix::SetSize(int h, int w)
132 {
133  MFEM_ASSERT(h >= 0 && w >= 0,
134  "invalid DenseMatrix size: " << h << " x " << w);
135  if (Height() == h && Width() == w)
136  {
137  return;
138  }
139  height = h;
140  width = w;
141  const int hw = h*w;
142  if (hw > data.Capacity())
143  {
144  data.Delete();
145  data.New(hw);
146  *this = 0.0; // init with zeroes
147  }
148 }
149 
150 double &DenseMatrix::Elem(int i, int j)
151 {
152  return (*this)(i,j);
153 }
154 
155 const double &DenseMatrix::Elem(int i, int j) const
156 {
157  return (*this)(i,j);
158 }
159 
160 void DenseMatrix::Mult(const double *x, double *y) const
161 {
162  kernels::Mult(height, width, Data(), x, y);
163 }
164 
165 void DenseMatrix::Mult(const Vector &x, Vector &y) const
166 {
167  MFEM_ASSERT(height == y.Size() && width == x.Size(),
168  "incompatible dimensions");
169 
170  Mult((const double *)x, (double *)y);
171 }
172 
173 double DenseMatrix::operator *(const DenseMatrix &m) const
174 {
175  MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
176  "incompatible dimensions");
177 
178  const int hw = height * width;
179  double a = 0.0;
180  for (int i = 0; i < hw; i++)
181  {
182  a += data[i] * m.data[i];
183  }
184 
185  return a;
186 }
187 
188 void DenseMatrix::MultTranspose(const double *x, double *y) const
189 {
190  double *d_col = Data();
191  for (int col = 0; col < width; col++)
192  {
193  double y_col = 0.0;
194  for (int row = 0; row < height; row++)
195  {
196  y_col += x[row]*d_col[row];
197  }
198  y[col] = y_col;
199  d_col += height;
200  }
201 }
202 
203 void DenseMatrix::MultTranspose(const Vector &x, Vector &y) const
204 {
205  MFEM_ASSERT(height == x.Size() && width == y.Size(),
206  "incompatible dimensions");
207 
208  MultTranspose((const double *)x, (double *)y);
209 }
210 
211 void DenseMatrix::AddMult(const Vector &x, Vector &y) const
212 {
213  MFEM_ASSERT(height == y.Size() && width == x.Size(),
214  "incompatible dimensions");
215 
216  const double *xp = x, *d_col = data;
217  double *yp = y;
218  for (int col = 0; col < width; col++)
219  {
220  double x_col = xp[col];
221  for (int row = 0; row < height; row++)
222  {
223  yp[row] += x_col*d_col[row];
224  }
225  d_col += height;
226  }
227 }
228 
230 {
231  MFEM_ASSERT(height == x.Size() && width == y.Size(),
232  "incompatible dimensions");
233 
234  const double *d_col = data;
235  for (int col = 0; col < width; col++)
236  {
237  double y_col = 0.0;
238  for (int row = 0; row < height; row++)
239  {
240  y_col += x[row]*d_col[row];
241  }
242  y[col] += y_col;
243  d_col += height;
244  }
245 }
246 
247 void DenseMatrix::AddMult_a(double a, 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  const double x_col = a*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  Vector &y) const
267 {
268  MFEM_ASSERT(height == x.Size() && width == y.Size(),
269  "incompatible dimensions");
270 
271  const double *d_col = data;
272  for (int col = 0; col < width; col++)
273  {
274  double y_col = 0.0;
275  for (int row = 0; row < height; row++)
276  {
277  y_col += x[row]*d_col[row];
278  }
279  y[col] += a * y_col;
280  d_col += height;
281  }
282 }
283 
284 double DenseMatrix::InnerProduct(const double *x, const double *y) const
285 {
286  double prod = 0.0;
287 
288  for (int i = 0; i < height; i++)
289  {
290  double Axi = 0.0;
291  for (int j = 0; j < width; j++)
292  {
293  Axi += (*this)(i,j) * x[j];
294  }
295  prod += y[i] * Axi;
296  }
297 
298  return prod;
299 }
300 
301 // LeftScaling this = diag(s) * this
303 {
304  double * it_data = data;
305  for (int j = 0; j < width; ++j)
306  {
307  for (int i = 0; i < height; ++i)
308  {
309  *(it_data++) *= s(i);
310  }
311  }
312 }
313 
314 // InvLeftScaling this = diag(1./s) * this
316 {
317  double * it_data = data;
318  for (int j = 0; j < width; ++j)
319  {
320  for (int i = 0; i < height; ++i)
321  {
322  *(it_data++) /= s(i);
323  }
324  }
325 }
326 
327 // RightScaling: this = this * diag(s);
329 {
330  double sj;
331  double * it_data = data;
332  for (int j = 0; j < width; ++j)
333  {
334  sj = s(j);
335  for (int i = 0; i < height; ++i)
336  {
337  *(it_data++) *= sj;
338  }
339  }
340 }
341 
342 // InvRightScaling: this = this * diag(1./s);
344 {
345  double * it_data = data;
346  for (int j = 0; j < width; ++j)
347  {
348  const double sj = 1./s(j);
349  for (int i = 0; i < height; ++i)
350  {
351  *(it_data++) *= sj;
352  }
353  }
354 }
355 
356 // SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
358 {
359  if (height != width || s.Size() != height)
360  {
361  mfem_error("DenseMatrix::SymmetricScaling: dimension mismatch");
362  }
363 
364  double * ss = new double[width];
365  double * it_s = s.GetData();
366  double * it_ss = ss;
367  for ( double * end_s = it_s + width; it_s != end_s; ++it_s)
368  {
369  *(it_ss++) = sqrt(*it_s);
370  }
371 
372  double * it_data = data;
373  for (int j = 0; j < width; ++j)
374  {
375  for (int i = 0; i < height; ++i)
376  {
377  *(it_data++) *= ss[i]*ss[j];
378  }
379  }
380 
381  delete[] ss;
382 }
383 
384 // InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
386 {
387  if (height != width || s.Size() != width)
388  {
389  mfem_error("DenseMatrix::InvSymmetricScaling: dimension mismatch");
390  }
391 
392  double * ss = new double[width];
393  double * it_s = s.GetData();
394  double * it_ss = ss;
395  for (double * end_s = it_s + width; it_s != end_s; ++it_s)
396  {
397  *(it_ss++) = 1./sqrt(*it_s);
398  }
399 
400  double * it_data = data;
401  for (int j = 0; j < width; ++j)
402  {
403  for (int i = 0; i < height; ++i)
404  {
405  *(it_data++) *= ss[i]*ss[j];
406  }
407  }
408 
409  delete[] ss;
410 }
411 
412 double DenseMatrix::Trace() const
413 {
414 #ifdef MFEM_DEBUG
415  if (Width() != Height())
416  {
417  mfem_error("DenseMatrix::Trace() : not a square matrix!");
418  }
419 #endif
420 
421  double t = 0.0;
422 
423  for (int i = 0; i < width; i++)
424  {
425  t += (*this)(i, i);
426  }
427 
428  return t;
429 }
430 
432 {
433  return new DenseMatrixInverse(*this);
434 }
435 
436 double DenseMatrix::Det() const
437 {
438  MFEM_ASSERT(Height() == Width() && Height() > 0,
439  "The matrix must be square and "
440  << "sized larger than zero to compute the determinant."
441  << " Height() = " << Height()
442  << ", Width() = " << Width());
443 
444  switch (Height())
445  {
446  case 1:
447  return data[0];
448 
449  case 2:
450  return data[0] * data[3] - data[1] * data[2];
451 
452  case 3:
453  {
454  const double *d = data;
455  return
456  d[0] * (d[4] * d[8] - d[5] * d[7]) +
457  d[3] * (d[2] * d[7] - d[1] * d[8]) +
458  d[6] * (d[1] * d[5] - d[2] * d[4]);
459  }
460  case 4:
461  {
462  const double *d = data;
463  return
464  d[ 0] * (d[ 5] * (d[10] * d[15] - d[11] * d[14]) -
465  d[ 9] * (d[ 6] * d[15] - d[ 7] * d[14]) +
466  d[13] * (d[ 6] * d[11] - d[ 7] * d[10])
467  ) -
468  d[ 4] * (d[ 1] * (d[10] * d[15] - d[11] * d[14]) -
469  d[ 9] * (d[ 2] * d[15] - d[ 3] * d[14]) +
470  d[13] * (d[ 2] * d[11] - d[ 3] * d[10])
471  ) +
472  d[ 8] * (d[ 1] * (d[ 6] * d[15] - d[ 7] * d[14]) -
473  d[ 5] * (d[ 2] * d[15] - d[ 3] * d[14]) +
474  d[13] * (d[ 2] * d[ 7] - d[ 3] * d[ 6])
475  ) -
476  d[12] * (d[ 1] * (d[ 6] * d[11] - d[ 7] * d[10]) -
477  d[ 5] * (d[ 2] * d[11] - d[ 3] * d[10]) +
478  d[ 9] * (d[ 2] * d[ 7] - d[ 3] * d[ 6])
479  );
480  }
481  default:
482  {
483  // In the general case we compute the determinant from the LU
484  // decomposition.
485  DenseMatrixInverse lu_factors(*this);
486 
487  return lu_factors.Det();
488  }
489  }
490  // not reachable
491 }
492 
493 double DenseMatrix::Weight() const
494 {
495  if (Height() == Width())
496  {
497  // return fabs(Det());
498  return Det();
499  }
500  else if ((Height() == 2) && (Width() == 1))
501  {
502  return sqrt(data[0] * data[0] + data[1] * data[1]);
503  }
504  else if ((Height() == 3) && (Width() == 1))
505  {
506  return sqrt(data[0] * data[0] + data[1] * data[1] + data[2] * data[2]);
507  }
508  else if ((Height() == 3) && (Width() == 2))
509  {
510  const double *d = data;
511  double E = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
512  double G = d[3] * d[3] + d[4] * d[4] + d[5] * d[5];
513  double F = d[0] * d[3] + d[1] * d[4] + d[2] * d[5];
514  return sqrt(E * G - F * F);
515  }
516  mfem_error("DenseMatrix::Weight(): mismatched or unsupported dimensions");
517  return 0.0;
518 }
519 
520 void DenseMatrix::Set(double alpha, const double *A)
521 {
522  const int s = Width()*Height();
523  for (int i = 0; i < s; i++)
524  {
525  data[i] = alpha*A[i];
526  }
527 }
528 
529 void DenseMatrix::Add(const double c, const DenseMatrix &A)
530 {
531  for (int j = 0; j < Width(); j++)
532  {
533  for (int i = 0; i < Height(); i++)
534  {
535  (*this)(i,j) += c * A(i,j);
536  }
537  }
538 }
539 
541 {
542  const int s = Height()*Width();
543  for (int i = 0; i < s; i++)
544  {
545  data[i] = c;
546  }
547  return *this;
548 }
549 
551 {
552  const int s = Height()*Width();
553  for (int i = 0; i < s; i++)
554  {
555  data[i] = d[i];
556  }
557  return *this;
558 }
559 
561 {
562  SetSize(m.height, m.width);
563 
564  const int hw = height * width;
565  for (int i = 0; i < hw; i++)
566  {
567  data[i] = m.data[i];
568  }
569 
570  return *this;
571 }
572 
574 {
575  kernels::Add(Height(), Width(), m, (double*)data);
576  return *this;
577 }
578 
580 {
581  MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
582  "incompatible matrix sizes.");
583  return *this += m.GetData();
584 }
585 
587 {
588  for (int j = 0; j < width; j++)
589  {
590  for (int i = 0; i < height; i++)
591  {
592  (*this)(i, j) -= m(i, j);
593  }
594  }
595 
596  return *this;
597 }
598 
600 {
601  int s = Height()*Width();
602  for (int i = 0; i < s; i++)
603  {
604  data[i] *= c;
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] = -data[i];
615  }
616 }
617 
619 {
620 #ifdef MFEM_DEBUG
621  if (Height() <= 0 || Height() != Width())
622  {
623  mfem_error("DenseMatrix::Invert(): dimension mismatch");
624  }
625 #endif
626 
627 #ifdef MFEM_USE_LAPACK
628  int *ipiv = new int[width];
629  int lwork = -1;
630  double qwork, *work;
631  int info;
632 
633  dgetrf_(&width, &width, data, &width, ipiv, &info);
634 
635  if (info)
636  {
637  mfem_error("DenseMatrix::Invert() : Error in DGETRF");
638  }
639 
640  dgetri_(&width, data, &width, ipiv, &qwork, &lwork, &info);
641 
642  lwork = (int) qwork;
643  work = new double[lwork];
644 
645  dgetri_(&width, data, &width, ipiv, work, &lwork, &info);
646 
647  if (info)
648  {
649  mfem_error("DenseMatrix::Invert() : Error in DGETRI");
650  }
651 
652  delete [] work;
653  delete [] ipiv;
654 #else
655  int c, i, j, n = Width();
656  double a, b;
657  Array<int> piv(n);
658 
659  for (c = 0; c < n; c++)
660  {
661  a = fabs((*this)(c, c));
662  i = c;
663  for (j = c + 1; j < n; j++)
664  {
665  b = fabs((*this)(j, c));
666  if (a < b)
667  {
668  a = b;
669  i = j;
670  }
671  }
672  if (a == 0.0)
673  {
674  mfem_error("DenseMatrix::Invert() : singular matrix");
675  }
676  piv[c] = i;
677  for (j = 0; j < n; j++)
678  {
679  mfem::Swap<double>((*this)(c, j), (*this)(i, j));
680  }
681 
682  a = (*this)(c, c) = 1.0 / (*this)(c, c);
683  for (j = 0; j < c; j++)
684  {
685  (*this)(c, j) *= a;
686  }
687  for (j++; j < n; j++)
688  {
689  (*this)(c, j) *= a;
690  }
691  for (i = 0; i < c; i++)
692  {
693  (*this)(i, c) = a * (b = -(*this)(i, c));
694  for (j = 0; j < c; j++)
695  {
696  (*this)(i, j) += b * (*this)(c, j);
697  }
698  for (j++; j < n; j++)
699  {
700  (*this)(i, j) += b * (*this)(c, j);
701  }
702  }
703  for (i++; i < n; i++)
704  {
705  (*this)(i, c) = a * (b = -(*this)(i, c));
706  for (j = 0; j < c; j++)
707  {
708  (*this)(i, j) += b * (*this)(c, j);
709  }
710  for (j++; j < n; j++)
711  {
712  (*this)(i, j) += b * (*this)(c, j);
713  }
714  }
715  }
716 
717  for (c = n - 1; c >= 0; c--)
718  {
719  j = piv[c];
720  for (i = 0; i < n; i++)
721  {
722  mfem::Swap<double>((*this)(i, c), (*this)(i, j));
723  }
724  }
725 #endif
726 }
727 
729 {
730  // Square root inverse using Denman--Beavers
731 #ifdef MFEM_DEBUG
732  if (Height() <= 0 || Height() != Width())
733  {
734  mfem_error("DenseMatrix::SquareRootInverse() matrix not square.");
735  }
736 #endif
737 
738  DenseMatrix tmp1(Height());
739  DenseMatrix tmp2(Height());
740  DenseMatrix tmp3(Height());
741 
742  tmp1 = (*this);
743  (*this) = 0.0;
744  for (int v = 0; v < Height() ; v++) { (*this)(v,v) = 1.0; }
745 
746  for (int j = 0; j < 10; j++)
747  {
748  for (int i = 0; i < 10; i++)
749  {
750  tmp2 = tmp1;
751  tmp3 = (*this);
752 
753  tmp2.Invert();
754  tmp3.Invert();
755 
756  tmp1 += tmp3;
757  (*this) += tmp2;
758 
759  tmp1 *= 0.5;
760  (*this) *= 0.5;
761  }
762  mfem::Mult((*this), tmp1, tmp2);
763  for (int v = 0; v < Height() ; v++) { tmp2(v,v) -= 1.0; }
764  if (tmp2.FNorm() < 1e-10) { break; }
765  }
766 
767  if (tmp2.FNorm() > 1e-10)
768  {
769  mfem_error("DenseMatrix::SquareRootInverse not converged");
770  }
771 }
772 
773 void DenseMatrix::Norm2(double *v) const
774 {
775  for (int j = 0; j < Width(); j++)
776  {
777  v[j] = 0.0;
778  for (int i = 0; i < Height(); i++)
779  {
780  v[j] += (*this)(i,j)*(*this)(i,j);
781  }
782  v[j] = sqrt(v[j]);
783  }
784 }
785 
787 {
788  int hw = Height()*Width();
789  const double *d = data;
790  double norm = 0.0, abs_entry;
791 
792  for (int i = 0; i < hw; i++)
793  {
794  abs_entry = fabs(d[i]);
795  if (norm < abs_entry)
796  {
797  norm = abs_entry;
798  }
799  }
800 
801  return norm;
802 }
803 
804 void DenseMatrix::FNorm(double &scale_factor, double &scaled_fnorm2) const
805 {
806  int i, hw = Height() * Width();
807  double max_norm = 0.0, entry, fnorm2;
808 
809  for (i = 0; i < hw; i++)
810  {
811  entry = fabs(data[i]);
812  if (entry > max_norm)
813  {
814  max_norm = entry;
815  }
816  }
817 
818  if (max_norm == 0.0)
819  {
820  scale_factor = scaled_fnorm2 = 0.0;
821  return;
822  }
823 
824  fnorm2 = 0.0;
825  for (i = 0; i < hw; i++)
826  {
827  entry = data[i] / max_norm;
828  fnorm2 += entry * entry;
829  }
830 
831  scale_factor = max_norm;
832  scaled_fnorm2 = fnorm2;
833 }
834 
836 {
837 #ifdef MFEM_USE_LAPACK
838  ev.SetSize(a.Width());
839 
840  char JOBZ = 'N';
841  char RANGE = 'A';
842  char UPLO = 'U';
843  int N = a.Width();
844  double *A = new double[N*N];
845  int LDA = N;
846  double VL = 0.0;
847  double VU = 1.0;
848  int IL = 0;
849  int IU = 1;
850  double ABSTOL = 0.0;
851  int M;
852  double *W = ev.GetData();
853  double *Z = NULL;
854  int LDZ = 1;
855  int *ISUPPZ = new int[2*N];
856  int LWORK = -1; // query optimal (double) workspace size
857  double QWORK;
858  double *WORK = NULL;
859  int LIWORK = -1; // query optimal (int) workspace size
860  int QIWORK;
861  int *IWORK = NULL;
862  int INFO;
863 
864  if (evect) // Compute eigenvectors too
865  {
866  evect->SetSize(N);
867 
868  JOBZ = 'V';
869  Z = evect->Data();
870  LDZ = N;
871  }
872 
873  int hw = a.Height() * a.Width();
874  double *data = a.Data();
875 
876  for (int i = 0; i < hw; i++)
877  {
878  A[i] = data[i];
879  }
880 
881  dsyevr_( &JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU,
882  &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, &QWORK, &LWORK,
883  &QIWORK, &LIWORK, &INFO );
884 
885  LWORK = (int) QWORK;
886  LIWORK = QIWORK;
887 
888  WORK = new double[LWORK];
889  IWORK = new int[LIWORK];
890 
891  dsyevr_( &JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU,
892  &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, WORK, &LWORK,
893  IWORK, &LIWORK, &INFO );
894 
895  if (INFO != 0)
896  {
897  mfem::err << "dsyevr_Eigensystem(...): DSYEVR error code: "
898  << INFO << endl;
899  mfem_error();
900  }
901 
902 #ifdef MFEM_DEBUG
903  if (M < N)
904  {
905  mfem::err << "dsyevr_Eigensystem(...):\n"
906  << " DSYEVR did not find all eigenvalues "
907  << M << "/" << N << endl;
908  mfem_error();
909  }
910  if (CheckFinite(W, N) > 0)
911  {
912  mfem_error("dsyevr_Eigensystem(...): inf/nan values in W");
913  }
914  if (CheckFinite(Z, N*N) > 0)
915  {
916  mfem_error("dsyevr_Eigensystem(...): inf/nan values in Z");
917  }
918  VU = 0.0;
919  for (IL = 0; IL < N; IL++)
920  for (IU = 0; IU <= IL; IU++)
921  {
922  VL = 0.0;
923  for (M = 0; M < N; M++)
924  {
925  VL += Z[M+IL*N] * Z[M+IU*N];
926  }
927  if (IU < IL)
928  {
929  VL = fabs(VL);
930  }
931  else
932  {
933  VL = fabs(VL-1.0);
934  }
935  if (VL > VU)
936  {
937  VU = VL;
938  }
939  if (VU > 0.5)
940  {
941  mfem::err << "dsyevr_Eigensystem(...):"
942  << " Z^t Z - I deviation = " << VU
943  << "\n W[max] = " << W[N-1] << ", W[min] = "
944  << W[0] << ", N = " << N << endl;
945  mfem_error();
946  }
947  }
948  if (VU > 1e-9)
949  {
950  mfem::err << "dsyevr_Eigensystem(...):"
951  << " Z^t Z - I deviation = " << VU
952  << "\n W[max] = " << W[N-1] << ", W[min] = "
953  << W[0] << ", N = " << N << endl;
954  }
955  if (VU > 1e-5)
956  {
957  mfem_error("dsyevr_Eigensystem(...): ERROR: ...");
958  }
959  VU = 0.0;
960  for (IL = 0; IL < N; IL++)
961  for (IU = 0; IU < N; IU++)
962  {
963  VL = 0.0;
964  for (M = 0; M < N; M++)
965  {
966  VL += Z[IL+M*N] * W[M] * Z[IU+M*N];
967  }
968  VL = fabs(VL-data[IL+N*IU]);
969  if (VL > VU)
970  {
971  VU = VL;
972  }
973  }
974  if (VU > 1e-9)
975  {
976  mfem::err << "dsyevr_Eigensystem(...):"
977  << " max matrix deviation = " << VU
978  << "\n W[max] = " << W[N-1] << ", W[min] = "
979  << W[0] << ", N = " << N << endl;
980  }
981  if (VU > 1e-5)
982  {
983  mfem_error("dsyevr_Eigensystem(...): ERROR: ...");
984  }
985 #endif
986 
987  delete [] IWORK;
988  delete [] WORK;
989  delete [] ISUPPZ;
990  delete [] A;
991 #else
992  MFEM_CONTRACT_VAR(a);
993  MFEM_CONTRACT_VAR(ev);
994  MFEM_CONTRACT_VAR(evect);
995 #endif
996 }
997 
999 {
1000 #ifdef MFEM_USE_LAPACK
1001  int N = a.Width();
1002  char JOBZ = 'N';
1003  char UPLO = 'U';
1004  int LDA = N;
1005  int LWORK = -1; /* query optimal workspace size */
1006  int INFO;
1007 
1008  ev.SetSize(N);
1009 
1010  double *A = NULL;
1011  double *W = ev.GetData();
1012  double *WORK = NULL;
1013  double QWORK;
1014 
1015  if (evect)
1016  {
1017  JOBZ = 'V';
1018  evect->SetSize(N);
1019  A = evect->Data();
1020  }
1021  else
1022  {
1023  A = new double[N*N];
1024  }
1025 
1026  int hw = a.Height() * a.Width();
1027  double *data = a.Data();
1028  for (int i = 0; i < hw; i++)
1029  {
1030  A[i] = data[i];
1031  }
1032 
1033  dsyev_(&JOBZ, &UPLO, &N, A, &LDA, W, &QWORK, &LWORK, &INFO);
1034 
1035  LWORK = (int) QWORK;
1036  WORK = new double[LWORK];
1037 
1038  dsyev_(&JOBZ, &UPLO, &N, A, &LDA, W, WORK, &LWORK, &INFO);
1039 
1040  if (INFO != 0)
1041  {
1042  mfem::err << "dsyev_Eigensystem: DSYEV error code: " << INFO << endl;
1043  mfem_error();
1044  }
1045 
1046  delete [] WORK;
1047  if (evect == NULL) { delete [] A; }
1048 #else
1049  MFEM_CONTRACT_VAR(a);
1050  MFEM_CONTRACT_VAR(ev);
1051  MFEM_CONTRACT_VAR(evect);
1052 #endif
1053 }
1054 
1055 void DenseMatrix::Eigensystem(Vector &ev, DenseMatrix *evect)
1056 {
1057 #ifdef MFEM_USE_LAPACK
1058 
1059  // dsyevr_Eigensystem(*this, ev, evect);
1060 
1061  dsyev_Eigensystem(*this, ev, evect);
1062 
1063 #else
1064 
1065  MFEM_CONTRACT_VAR(ev);
1066  MFEM_CONTRACT_VAR(evect);
1067  mfem_error("DenseMatrix::Eigensystem: Compiled without LAPACK");
1068 
1069 #endif
1070 }
1071 
1073  DenseMatrix *evect)
1074 {
1075 #ifdef MFEM_USE_LAPACK
1076  int N = a.Width();
1077  int ITYPE = 1;
1078  char JOBZ = 'N';
1079  char UPLO = 'U';
1080  int LDA = N;
1081  int LDB = N;
1082  int LWORK = -1; /* query optimal workspace size */
1083  int INFO;
1084 
1085  ev.SetSize(N);
1086 
1087  double *A = NULL;
1088  double *B = new double[N*N];
1089  double *W = ev.GetData();
1090  double *WORK = NULL;
1091  double QWORK;
1092 
1093  if (evect)
1094  {
1095  JOBZ = 'V';
1096  evect->SetSize(N);
1097  A = evect->Data();
1098  }
1099  else
1100  {
1101  A = new double[N*N];
1102  }
1103 
1104  int hw = a.Height() * a.Width();
1105  double *a_data = a.Data();
1106  double *b_data = b.Data();
1107  for (int i = 0; i < hw; i++)
1108  {
1109  A[i] = a_data[i];
1110  B[i] = b_data[i];
1111  }
1112 
1113  dsygv_(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, &QWORK, &LWORK, &INFO);
1114 
1115  LWORK = (int) QWORK;
1116  WORK = new double[LWORK];
1117 
1118  dsygv_(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, WORK, &LWORK, &INFO);
1119 
1120  if (INFO != 0)
1121  {
1122  mfem::err << "dsygv_Eigensystem: DSYGV error code: " << INFO << endl;
1123  mfem_error();
1124  }
1125 
1126  delete [] WORK;
1127  delete [] B;
1128  if (evect == NULL) { delete [] A; }
1129 #else
1130  MFEM_CONTRACT_VAR(a);
1131  MFEM_CONTRACT_VAR(b);
1132  MFEM_CONTRACT_VAR(ev);
1133  MFEM_CONTRACT_VAR(evect);
1134 #endif
1135 }
1136 
1137 void DenseMatrix::Eigensystem(DenseMatrix &b, Vector &ev,
1138  DenseMatrix *evect)
1139 {
1140 #ifdef MFEM_USE_LAPACK
1141 
1142  dsygv_Eigensystem(*this, b, ev, evect);
1143 
1144 #else
1145  MFEM_CONTRACT_VAR(b);
1146  MFEM_CONTRACT_VAR(ev);
1147  MFEM_CONTRACT_VAR(evect);
1148  mfem_error("DenseMatrix::Eigensystem(generalized): Compiled without LAPACK");
1149 #endif
1150 }
1151 
1153 {
1154 #ifdef MFEM_USE_LAPACK
1155  DenseMatrix copy_of_this = *this;
1156  char jobu = 'N';
1157  char jobvt = 'N';
1158  int m = Height();
1159  int n = Width();
1160  double *a = copy_of_this.data;
1161  sv.SetSize(min(m, n));
1162  double *s = sv;
1163  double *u = NULL;
1164  double *vt = NULL;
1165  double *work = NULL;
1166  int lwork = -1;
1167  int info;
1168  double qwork;
1169 
1170  dgesvd_(&jobu, &jobvt, &m, &n, a, &m,
1171  s, u, &m, vt, &n, &qwork, &lwork, &info);
1172 
1173  lwork = (int) qwork;
1174  work = new double[lwork];
1175 
1176  dgesvd_(&jobu, &jobvt, &m, &n, a, &m,
1177  s, u, &m, vt, &n, work, &lwork, &info);
1178 
1179  delete [] work;
1180  if (info)
1181  {
1182  mfem::err << "DenseMatrix::SingularValues : info = " << info << endl;
1183  mfem_error();
1184  }
1185 #else
1186  MFEM_CONTRACT_VAR(sv);
1187  // compiling without lapack
1188  mfem_error("DenseMatrix::SingularValues: Compiled without LAPACK");
1189 #endif
1190 }
1191 
1192 int DenseMatrix::Rank(double tol) const
1193 {
1194  int rank=0;
1195  Vector sv(min(Height(), Width()));
1196  SingularValues(sv);
1197 
1198  for (int i=0; i < sv.Size(); ++i)
1199  if (sv(i) >= tol)
1200  {
1201  ++rank;
1202  }
1203 
1204  return rank;
1205 }
1206 
1207 double DenseMatrix::CalcSingularvalue(const int i) const
1208 {
1209  MFEM_ASSERT(Height() == Width() && Height() > 0 && Height() < 4,
1210  "The matrix must be square and sized 1, 2, or 3 to compute the"
1211  " singular values."
1212  << " Height() = " << Height()
1213  << ", Width() = " << Width());
1214 
1215  const int n = Height();
1216  const double *d = data;
1217 
1218  if (n == 1)
1219  {
1220  return d[0];
1221  }
1222  else if (n == 2)
1223  {
1224  return kernels::CalcSingularvalue<2>(d,i);
1225  }
1226  else
1227  {
1228  return kernels::CalcSingularvalue<3>(d,i);
1229  }
1230 }
1231 
1232 void DenseMatrix::CalcEigenvalues(double *lambda, double *vec) const
1233 {
1234 #ifdef MFEM_DEBUG
1235  if (Height() != Width() || Height() < 2 || Height() > 3)
1236  {
1237  mfem_error("DenseMatrix::CalcEigenvalues");
1238  }
1239 #endif
1240 
1241  const int n = Height();
1242  const double *d = data;
1243 
1244  if (n == 2)
1245  {
1246  kernels::CalcEigenvalues<2>(d, lambda, vec);
1247  }
1248  else
1249  {
1250  kernels::CalcEigenvalues<3>(d, lambda, vec);
1251  }
1252 }
1253 
1254 void DenseMatrix::GetRow(int r, Vector &row) const
1255 {
1256  int m = Height();
1257  int n = Width();
1258  row.SetSize(n);
1259 
1260  const double* rp = data + r;
1261  double* vp = row.GetData();
1262 
1263  for (int i = 0; i < n; i++)
1264  {
1265  vp[i] = *rp;
1266  rp += m;
1267  }
1268 }
1269 
1270 void DenseMatrix::GetColumn(int c, Vector &col) const
1271 {
1272  int m = Height();
1273  col.SetSize(m);
1274 
1275  double *cp = Data() + c * m;
1276  double *vp = col.GetData();
1277 
1278  for (int i = 0; i < m; i++)
1279  {
1280  vp[i] = cp[i];
1281  }
1282 }
1283 
1285 {
1286  if (height != width)
1287  {
1288  mfem_error("DenseMatrix::GetDiag\n");
1289  }
1290  d.SetSize(height);
1291 
1292  for (int i = 0; i < height; ++i)
1293  {
1294  d(i) = (*this)(i,i);
1295  }
1296 }
1297 
1299 {
1300  if (height != width)
1301  {
1302  mfem_error("DenseMatrix::Getl1Diag\n");
1303  }
1304  l.SetSize(height);
1305 
1306  l = 0.0;
1307 
1308  for (int j = 0; j < width; ++j)
1309  for (int i = 0; i < height; ++i)
1310  {
1311  l(i) += fabs((*this)(i,j));
1312  }
1313 }
1314 
1316 {
1317  l.SetSize(height);
1318  for (int i = 0; i < height; i++)
1319  {
1320  double d = 0.0;
1321  for (int j = 0; j < width; j++)
1322  {
1323  d += operator()(i, j);
1324  }
1325  l(i) = d;
1326  }
1327 }
1328 
1329 void DenseMatrix::Diag(double c, int n)
1330 {
1331  SetSize(n);
1332 
1333  const int N = n*n;
1334  for (int i = 0; i < N; i++)
1335  {
1336  data[i] = 0.0;
1337  }
1338  for (int i = 0; i < n; i++)
1339  {
1340  data[i*(n+1)] = c;
1341  }
1342 }
1343 
1344 void DenseMatrix::Diag(double *diag, int n)
1345 {
1346  SetSize(n);
1347 
1348  int i, N = n*n;
1349  for (i = 0; i < N; i++)
1350  {
1351  data[i] = 0.0;
1352  }
1353  for (i = 0; i < n; i++)
1354  {
1355  data[i*(n+1)] = diag[i];
1356  }
1357 }
1358 
1360 {
1361  int i, j;
1362  double t;
1363 
1364  if (Width() == Height())
1365  {
1366  for (i = 0; i < Height(); i++)
1367  for (j = i+1; j < Width(); j++)
1368  {
1369  t = (*this)(i,j);
1370  (*this)(i,j) = (*this)(j,i);
1371  (*this)(j,i) = t;
1372  }
1373  }
1374  else
1375  {
1376  DenseMatrix T(*this,'t');
1377  (*this) = T;
1378  }
1379 }
1380 
1382 {
1383  SetSize(A.Width(),A.Height());
1384 
1385  for (int i = 0; i < Height(); i++)
1386  for (int j = 0; j < Width(); j++)
1387  {
1388  (*this)(i,j) = A(j,i);
1389  }
1390 }
1391 
1393 {
1394 #ifdef MFEM_DEBUG
1395  if (Width() != Height())
1396  {
1397  mfem_error("DenseMatrix::Symmetrize() : not a square matrix!");
1398  }
1399 #endif
1401 }
1402 
1404 {
1405  for (int i = 0; i < Height(); i++)
1406  {
1407  double L = 0.0;
1408  for (int j = 0; j < Width(); j++)
1409  {
1410  L += (*this)(i, j);
1411  (*this)(i, j) = 0.0;
1412  }
1413  (*this)(i, i) = L;
1414  }
1415 }
1416 
1418 {
1419  int n = Height();
1420 
1421 #ifdef MFEM_DEBUG
1422  if ((Width() != 2 || curl.Width() != 1 || 2*n != curl.Height()) &&
1423  (Width() != 3 || curl.Width() != 3 || 3*n != curl.Height()))
1424  {
1425  mfem_error("DenseMatrix::GradToCurl(...): dimension mismatch");
1426  }
1427 #endif
1428 
1429  if (Width() == 2)
1430  {
1431  for (int i = 0; i < n; i++)
1432  {
1433  // (x,y) is grad of Ui
1434  double x = (*this)(i,0);
1435  double y = (*this)(i,1);
1436 
1437  int j = i+n;
1438 
1439  // curl of (Ui,0)
1440  curl(i,0) = -y;
1441 
1442  // curl of (0,Ui)
1443  curl(j,0) = x;
1444  }
1445  }
1446  else
1447  {
1448  for (int i = 0; i < n; i++)
1449  {
1450  // (x,y,z) is grad of Ui
1451  double x = (*this)(i,0);
1452  double y = (*this)(i,1);
1453  double z = (*this)(i,2);
1454 
1455  int j = i+n;
1456  int k = j+n;
1457 
1458  // curl of (Ui,0,0)
1459  curl(i,0) = 0.;
1460  curl(i,1) = z;
1461  curl(i,2) = -y;
1462 
1463  // curl of (0,Ui,0)
1464  curl(j,0) = -z;
1465  curl(j,1) = 0.;
1466  curl(j,2) = x;
1467 
1468  // curl of (0,0,Ui)
1469  curl(k,0) = y;
1470  curl(k,1) = -x;
1471  curl(k,2) = 0.;
1472  }
1473  }
1474 }
1475 
1477 {
1478  MFEM_ASSERT(Width()*Height() == div.Size(), "incompatible Vector 'div'!");
1479 
1480  // div(dof*j+i) <-- (*this)(i,j)
1481 
1482  const int n = height * width;
1483  double *ddata = div.GetData();
1484 
1485  for (int i = 0; i < n; i++)
1486  {
1487  ddata[i] = data[i];
1488  }
1489 }
1490 
1491 void DenseMatrix::CopyRows(const DenseMatrix &A, int row1, int row2)
1492 {
1493  SetSize(row2 - row1 + 1, A.Width());
1494 
1495  for (int j = 0; j < Width(); j++)
1496  {
1497  for (int i = row1; i <= row2; i++)
1498  {
1499  (*this)(i-row1,j) = A(i,j);
1500  }
1501  }
1502 }
1503 
1504 void DenseMatrix::CopyCols(const DenseMatrix &A, int col1, int col2)
1505 {
1506  SetSize(A.Height(), col2 - col1 + 1);
1507 
1508  for (int j = col1; j <= col2; j++)
1509  {
1510  for (int i = 0; i < Height(); i++)
1511  {
1512  (*this)(i,j-col1) = A(i,j);
1513  }
1514  }
1515 }
1516 
1517 void DenseMatrix::CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
1518 {
1519  SetSize(m,n);
1520 
1521  for (int j = 0; j < n; j++)
1522  {
1523  for (int i = 0; i < m; i++)
1524  {
1525  (*this)(i,j) = A(Aro+i,Aco+j);
1526  }
1527  }
1528 }
1529 
1530 void DenseMatrix::CopyMN(const DenseMatrix &A, int row_offset, int col_offset)
1531 {
1532  double *v = A.Data();
1533 
1534  for (int j = 0; j < A.Width(); j++)
1535  {
1536  for (int i = 0; i < A.Height(); i++)
1537  {
1538  (*this)(row_offset+i,col_offset+j) = *(v++);
1539  }
1540  }
1541 }
1542 
1543 void DenseMatrix::CopyMNt(const DenseMatrix &A, int row_offset, int col_offset)
1544 {
1545  double *v = A.Data();
1546 
1547  for (int i = 0; i < A.Width(); i++)
1548  {
1549  for (int j = 0; j < A.Height(); j++)
1550  {
1551  (*this)(row_offset+i,col_offset+j) = *(v++);
1552  }
1553  }
1554 }
1555 
1556 void DenseMatrix::CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco,
1557  int row_offset, int col_offset)
1558 {
1559  MFEM_VERIFY(row_offset+m <= this->Height() && col_offset+n <= this->Width(),
1560  "this DenseMatrix is too small to accomodate the submatrix. "
1561  << "row_offset = " << row_offset
1562  << ", m = " << m
1563  << ", this->Height() = " << this->Height()
1564  << ", col_offset = " << col_offset
1565  << ", n = " << n
1566  << ", this->Width() = " << this->Width()
1567  );
1568  MFEM_VERIFY(Aro+m <= A.Height() && Aco+n <= A.Width(),
1569  "The A DenseMatrix is too small to accomodate the submatrix. "
1570  << "Aro = " << Aro
1571  << ", m = " << m
1572  << ", A.Height() = " << A.Height()
1573  << ", Aco = " << Aco
1574  << ", n = " << n
1575  << ", A.Width() = " << A.Width()
1576  );
1577 
1578  for (int j = 0; j < n; j++)
1579  {
1580  for (int i = 0; i < m; i++)
1581  {
1582  (*this)(row_offset+i,col_offset+j) = A(Aro+i,Aco+j);
1583  }
1584  }
1585 }
1586 
1587 void DenseMatrix::CopyMNDiag(double c, int n, int row_offset, int col_offset)
1588 {
1589  for (int i = 0; i < n; i++)
1590  {
1591  for (int j = i+1; j < n; j++)
1592  {
1593  (*this)(row_offset+i,col_offset+j) =
1594  (*this)(row_offset+j,col_offset+i) = 0.0;
1595  }
1596  }
1597 
1598  for (int i = 0; i < n; i++)
1599  {
1600  (*this)(row_offset+i,col_offset+i) = c;
1601  }
1602 }
1603 
1604 void DenseMatrix::CopyMNDiag(double *diag, int n, int row_offset,
1605  int col_offset)
1606 {
1607  for (int i = 0; i < n; i++)
1608  {
1609  for (int j = i+1; j < n; j++)
1610  {
1611  (*this)(row_offset+i,col_offset+j) =
1612  (*this)(row_offset+j,col_offset+i) = 0.0;
1613  }
1614  }
1615 
1616  for (int i = 0; i < n; i++)
1617  {
1618  (*this)(row_offset+i,col_offset+i) = diag[i];
1619  }
1620 }
1621 
1622 void DenseMatrix::CopyExceptMN(const DenseMatrix &A, int m, int n)
1623 {
1624  SetSize(A.Width()-1,A.Height()-1);
1625 
1626  int i, j, i_off = 0, j_off = 0;
1627 
1628  for (j = 0; j < A.Width(); j++)
1629  {
1630  if ( j == n )
1631  {
1632  j_off = 1;
1633  continue;
1634  }
1635  for (i = 0; i < A.Height(); i++)
1636  {
1637  if ( i == m )
1638  {
1639  i_off = 1;
1640  continue;
1641  }
1642  (*this)(i-i_off,j-j_off) = A(i,j);
1643  }
1644  i_off = 0;
1645  }
1646 }
1647 
1648 void DenseMatrix::AddMatrix(DenseMatrix &A, int ro, int co)
1649 {
1650  int h, ah, aw;
1651  double *p, *ap;
1652 
1653  h = Height();
1654  ah = A.Height();
1655  aw = A.Width();
1656 
1657 #ifdef MFEM_DEBUG
1658  if (co+aw > Width() || ro+ah > h)
1659  {
1660  mfem_error("DenseMatrix::AddMatrix(...) 1 : dimension mismatch");
1661  }
1662 #endif
1663 
1664  p = data + ro + co * h;
1665  ap = A.data;
1666 
1667  for (int c = 0; c < aw; c++)
1668  {
1669  for (int r = 0; r < ah; r++)
1670  {
1671  p[r] += ap[r];
1672  }
1673  p += h;
1674  ap += ah;
1675  }
1676 }
1677 
1678 void DenseMatrix::AddMatrix(double a, const DenseMatrix &A, int ro, int co)
1679 {
1680  int h, ah, aw;
1681  double *p, *ap;
1682 
1683  h = Height();
1684  ah = A.Height();
1685  aw = A.Width();
1686 
1687 #ifdef MFEM_DEBUG
1688  if (co+aw > Width() || ro+ah > h)
1689  {
1690  mfem_error("DenseMatrix::AddMatrix(...) 2 : dimension mismatch");
1691  }
1692 #endif
1693 
1694  p = data + ro + co * h;
1695  ap = A.Data();
1696 
1697  for (int c = 0; c < aw; c++)
1698  {
1699  for (int r = 0; r < ah; r++)
1700  {
1701  p[r] += a * ap[r];
1702  }
1703  p += h;
1704  ap += ah;
1705  }
1706 }
1707 
1708 void DenseMatrix::AddToVector(int offset, Vector &v) const
1709 {
1710  const int n = height * width;
1711  double *vdata = v.GetData() + offset;
1712 
1713  for (int i = 0; i < n; i++)
1714  {
1715  vdata[i] += data[i];
1716  }
1717 }
1718 
1719 void DenseMatrix::GetFromVector(int offset, const Vector &v)
1720 {
1721  const int n = height * width;
1722  const double *vdata = v.GetData() + offset;
1723 
1724  for (int i = 0; i < n; i++)
1725  {
1726  data[i] = vdata[i];
1727  }
1728 }
1729 
1731 {
1732  const int n = Height();
1733 
1734 #ifdef MFEM_DEBUG
1735  if (dofs.Size() != n || Width() != n)
1736  {
1737  mfem_error("DenseMatrix::AdjustDofDirection(...): dimension mismatch");
1738  }
1739 #endif
1740 
1741  int *dof = dofs;
1742  for (int i = 0; i < n-1; i++)
1743  {
1744  const int s = (dof[i] < 0) ? (-1) : (1);
1745  for (int j = i+1; j < n; j++)
1746  {
1747  const int t = (dof[j] < 0) ? (-s) : (s);
1748  if (t < 0)
1749  {
1750  (*this)(i,j) = -(*this)(i,j);
1751  (*this)(j,i) = -(*this)(j,i);
1752  }
1753  }
1754  }
1755 }
1756 
1757 void DenseMatrix::SetRow(int row, double value)
1758 {
1759  for (int j = 0; j < Width(); j++)
1760  {
1761  (*this)(row, j) = value;
1762  }
1763 }
1764 
1765 void DenseMatrix::SetCol(int col, double value)
1766 {
1767  for (int i = 0; i < Height(); i++)
1768  {
1769  (*this)(i, col) = value;
1770  }
1771 }
1772 
1773 void DenseMatrix::SetRow(int r, const double* row)
1774 {
1775  MFEM_ASSERT(row != nullptr, "supplied row pointer is null");
1776  for (int j = 0; j < Width(); j++)
1777  {
1778  (*this)(r, j) = row[j];
1779  }
1780 }
1781 
1782 void DenseMatrix::SetRow(int r, const Vector &row)
1783 {
1784  MFEM_ASSERT(Width() == row.Size(), "");
1785  SetRow(r, row.GetData());
1786 }
1787 
1788 void DenseMatrix::SetCol(int c, const double* col)
1789 {
1790  MFEM_ASSERT(col != nullptr, "supplied column pointer is null");
1791  for (int i = 0; i < Height(); i++)
1792  {
1793  (*this)(i, c) = col[i];
1794  }
1795 }
1796 
1797 void DenseMatrix::SetCol(int c, const Vector &col)
1798 {
1799  MFEM_ASSERT(Height() == col.Size(), "");
1800  SetCol(c, col.GetData());
1801 }
1802 
1803 void DenseMatrix::Threshold(double eps)
1804 {
1805  for (int col = 0; col < Width(); col++)
1806  {
1807  for (int row = 0; row < Height(); row++)
1808  {
1809  if (std::abs(operator()(row,col)) <= eps)
1810  {
1811  operator()(row,col) = 0.0;
1812  }
1813  }
1814  }
1815 }
1816 
1817 void DenseMatrix::Print(std::ostream &os, int width_) const
1818 {
1819  // save current output flags
1820  ios::fmtflags old_flags = os.flags();
1821  // output flags = scientific + show sign
1822  os << setiosflags(ios::scientific | ios::showpos);
1823  for (int i = 0; i < height; i++)
1824  {
1825  os << "[row " << i << "]\n";
1826  for (int j = 0; j < width; j++)
1827  {
1828  os << (*this)(i,j);
1829  if (j+1 == width || (j+1) % width_ == 0)
1830  {
1831  os << '\n';
1832  }
1833  else
1834  {
1835  os << ' ';
1836  }
1837  }
1838  }
1839  // reset output flags to original values
1840  os.flags(old_flags);
1841 }
1842 
1843 void DenseMatrix::PrintMatlab(std::ostream &os) const
1844 {
1845  // save current output flags
1846  ios::fmtflags old_flags = os.flags();
1847  // output flags = scientific + show sign
1848  os << setiosflags(ios::scientific | ios::showpos);
1849  for (int i = 0; i < height; i++)
1850  {
1851  for (int j = 0; j < width; j++)
1852  {
1853  os << (*this)(i,j);
1854  os << ' ';
1855  }
1856  os << "\n";
1857  }
1858  // reset output flags to original values
1859  os.flags(old_flags);
1860 }
1861 
1862 void DenseMatrix::PrintT(std::ostream &os, int width_) const
1863 {
1864  // save current output flags
1865  ios::fmtflags old_flags = os.flags();
1866  // output flags = scientific + show sign
1867  os << setiosflags(ios::scientific | ios::showpos);
1868  for (int j = 0; j < width; j++)
1869  {
1870  os << "[col " << j << "]\n";
1871  for (int i = 0; i < height; i++)
1872  {
1873  os << (*this)(i,j);
1874  if (i+1 == height || (i+1) % width_ == 0)
1875  {
1876  os << '\n';
1877  }
1878  else
1879  {
1880  os << ' ';
1881  }
1882  }
1883  }
1884  // reset output flags to original values
1885  os.flags(old_flags);
1886 }
1887 
1889 {
1890  DenseMatrix copy(*this), C(width);
1891  Invert();
1892  mfem::Mult(*this, copy, C);
1893 
1894  for (int i = 0; i < width; i++)
1895  {
1896  C(i,i) -= 1.0;
1897  }
1898  mfem::out << "size = " << width << ", i_max = " << C.MaxMaxNorm()
1899  << ", cond_F = " << FNorm()*copy.FNorm() << endl;
1900 }
1901 
1903 {
1904  mfem::Swap(width, other.width);
1905  mfem::Swap(height, other.height);
1906  mfem::Swap(data, other.data);
1907 }
1908 
1910 {
1911  data.Delete();
1912 }
1913 
1914 
1915 
1916 void Add(const DenseMatrix &A, const DenseMatrix &B,
1917  double alpha, DenseMatrix &C)
1918 {
1919  kernels::Add(C.Height(), C.Width(), alpha, A.Data(), B.Data(), C.Data());
1920 }
1921 
1922 void Add(double alpha, const double *A,
1923  double beta, const double *B, DenseMatrix &C)
1924 {
1925  kernels::Add(C.Height(), C.Width(), alpha, A, beta, B, C.Data());
1926 }
1927 
1928 void Add(double alpha, const DenseMatrix &A,
1929  double beta, const DenseMatrix &B, DenseMatrix &C)
1930 {
1931  MFEM_ASSERT(A.Height() == C.Height(), "");
1932  MFEM_ASSERT(B.Height() == C.Height(), "");
1933  MFEM_ASSERT(A.Width() == C.Width(), "");
1934  MFEM_ASSERT(B.Width() == C.Width(), "");
1935  Add(alpha, A.GetData(), beta, B.GetData(), C);
1936 }
1937 
1938 bool LinearSolve(DenseMatrix& A, double* X, double TOL)
1939 {
1940  MFEM_VERIFY(A.IsSquare(), "A must be a square matrix!");
1941  MFEM_ASSERT(A.NumCols() > 0, "supplied matrix, A, is empty!");
1942  MFEM_ASSERT(X != nullptr, "supplied vector, X, is null!");
1943 
1944  int N = A.NumCols();
1945 
1946  switch (N)
1947  {
1948  case 1:
1949  {
1950  double det = A(0,0);
1951  if (std::abs(det) <= TOL) { return false; } // singular
1952 
1953  X[0] /= det;
1954  break;
1955  }
1956  case 2:
1957  {
1958  double det = A.Det();
1959  if (std::abs(det) <= TOL) { return false; } // singular
1960 
1961  double invdet = 1. / det;
1962 
1963  double b0 = X[0];
1964  double b1 = X[1];
1965 
1966  X[0] = ( A(1,1)*b0 - A(0,1)*b1) * invdet;
1967  X[1] = (-A(1,0)*b0 + A(0,0)*b1) * invdet;
1968  break;
1969  }
1970  default:
1971  {
1972  // default to LU factorization for the general case
1973  Array<int> ipiv(N);
1974  LUFactors lu(A.Data(), ipiv);
1975 
1976  if (!lu.Factor(N,TOL)) { return false; } // singular
1977 
1978  lu.Solve(N, 1, X);
1979  }
1980 
1981  } // END switch
1982 
1983  return true;
1984 }
1985 
1986 void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
1987 {
1988  MFEM_ASSERT(a.Height() == b.Height() && a.Width() == c.Width() &&
1989  b.Width() == c.Height(), "incompatible dimensions");
1990 
1991 #ifdef MFEM_USE_LAPACK
1992  static char transa = 'N', transb = 'N';
1993  static double alpha = 1.0, beta = 0.0;
1994  int m = b.Height(), n = c.Width(), k = b.Width();
1995 
1996  dgemm_(&transa, &transb, &m, &n, &k, &alpha, b.Data(), &m,
1997  c.Data(), &k, &beta, a.Data(), &m);
1998 #else
1999  const int ah = a.Height();
2000  const int aw = a.Width();
2001  const int bw = b.Width();
2002  double *ad = a.Data();
2003  const double *bd = b.Data();
2004  const double *cd = c.Data();
2005  kernels::Mult(ah,aw,bw,bd,cd,ad);
2006 #endif
2007 }
2008 
2009 void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c,
2010  DenseMatrix &a)
2011 {
2012  MFEM_ASSERT(a.Height() == b.Height() && a.Width() == c.Width() &&
2013  b.Width() == c.Height(), "incompatible dimensions");
2014 
2015 #ifdef MFEM_USE_LAPACK
2016  static char transa = 'N', transb = 'N';
2017  static double beta = 1.0;
2018  int m = b.Height(), n = c.Width(), k = b.Width();
2019 
2020  dgemm_(&transa, &transb, &m, &n, &k, &alpha, b.Data(), &m,
2021  c.Data(), &k, &beta, a.Data(), &m);
2022 #else
2023  const int ah = a.Height();
2024  const int aw = a.Width();
2025  const int bw = b.Width();
2026  double *ad = a.Data();
2027  const double *bd = b.Data();
2028  const double *cd = c.Data();
2029  for (int j = 0; j < aw; j++)
2030  {
2031  for (int k = 0; k < bw; k++)
2032  {
2033  for (int i = 0; i < ah; i++)
2034  {
2035  ad[i+j*ah] += alpha * bd[i+k*ah] * cd[k+j*bw];
2036  }
2037  }
2038  }
2039 #endif
2040 }
2041 
2042 void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
2043 {
2044  MFEM_ASSERT(a.Height() == b.Height() && a.Width() == c.Width() &&
2045  b.Width() == c.Height(), "incompatible dimensions");
2046 
2047 #ifdef MFEM_USE_LAPACK
2048  static char transa = 'N', transb = 'N';
2049  static double alpha = 1.0, beta = 1.0;
2050  int m = b.Height(), n = c.Width(), k = b.Width();
2051 
2052  dgemm_(&transa, &transb, &m, &n, &k, &alpha, b.Data(), &m,
2053  c.Data(), &k, &beta, a.Data(), &m);
2054 #else
2055  const int ah = a.Height();
2056  const int aw = a.Width();
2057  const int bw = b.Width();
2058  double *ad = a.Data();
2059  const double *bd = b.Data();
2060  const double *cd = c.Data();
2061  for (int j = 0; j < aw; j++)
2062  {
2063  for (int k = 0; k < bw; k++)
2064  {
2065  for (int i = 0; i < ah; i++)
2066  {
2067  ad[i+j*ah] += bd[i+k*ah] * cd[k+j*bw];
2068  }
2069  }
2070  }
2071 #endif
2072 }
2073 
2075 {
2076 #ifdef MFEM_DEBUG
2077  if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
2078  {
2079  mfem_error("CalcAdjugate(...): unsupported dimensions");
2080  }
2081  if (a.Width() != adja.Height() || a.Height() != adja.Width())
2082  {
2083  mfem_error("CalcAdjugate(...): dimension mismatch");
2084  }
2085 #endif
2086 
2087  if (a.Width() < a.Height())
2088  {
2089  const double *d = a.Data();
2090  double *ad = adja.Data();
2091  if (a.Width() == 1)
2092  {
2093  // N x 1, N = 2,3
2094  ad[0] = d[0];
2095  ad[1] = d[1];
2096  if (a.Height() == 3)
2097  {
2098  ad[2] = d[2];
2099  }
2100  }
2101  else
2102  {
2103  // 3 x 2
2104  double e, g, f;
2105  e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
2106  g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
2107  f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
2108 
2109  ad[0] = d[0]*g - d[3]*f;
2110  ad[1] = d[3]*e - d[0]*f;
2111  ad[2] = d[1]*g - d[4]*f;
2112  ad[3] = d[4]*e - d[1]*f;
2113  ad[4] = d[2]*g - d[5]*f;
2114  ad[5] = d[5]*e - d[2]*f;
2115  }
2116  return;
2117  }
2118 
2119  if (a.Width() == 1)
2120  {
2121  adja(0,0) = 1.0;
2122  }
2123  else if (a.Width() == 2)
2124  {
2125  adja(0,0) = a(1,1);
2126  adja(0,1) = -a(0,1);
2127  adja(1,0) = -a(1,0);
2128  adja(1,1) = a(0,0);
2129  }
2130  else
2131  {
2132  adja(0,0) = a(1,1)*a(2,2)-a(1,2)*a(2,1);
2133  adja(0,1) = a(0,2)*a(2,1)-a(0,1)*a(2,2);
2134  adja(0,2) = a(0,1)*a(1,2)-a(0,2)*a(1,1);
2135 
2136  adja(1,0) = a(1,2)*a(2,0)-a(1,0)*a(2,2);
2137  adja(1,1) = a(0,0)*a(2,2)-a(0,2)*a(2,0);
2138  adja(1,2) = a(0,2)*a(1,0)-a(0,0)*a(1,2);
2139 
2140  adja(2,0) = a(1,0)*a(2,1)-a(1,1)*a(2,0);
2141  adja(2,1) = a(0,1)*a(2,0)-a(0,0)*a(2,1);
2142  adja(2,2) = a(0,0)*a(1,1)-a(0,1)*a(1,0);
2143  }
2144 }
2145 
2147 {
2148 #ifdef MFEM_DEBUG
2149  if (a.Height() != a.Width() || adjat.Height() != adjat.Width() ||
2150  a.Width() != adjat.Width() || a.Width() < 1 || a.Width() > 3)
2151  {
2152  mfem_error("CalcAdjugateTranspose(...): dimension mismatch");
2153  }
2154 #endif
2155  if (a.Width() == 1)
2156  {
2157  adjat(0,0) = 1.0;
2158  }
2159  else if (a.Width() == 2)
2160  {
2161  adjat(0,0) = a(1,1);
2162  adjat(1,0) = -a(0,1);
2163  adjat(0,1) = -a(1,0);
2164  adjat(1,1) = a(0,0);
2165  }
2166  else
2167  {
2168  adjat(0,0) = a(1,1)*a(2,2)-a(1,2)*a(2,1);
2169  adjat(1,0) = a(0,2)*a(2,1)-a(0,1)*a(2,2);
2170  adjat(2,0) = a(0,1)*a(1,2)-a(0,2)*a(1,1);
2171 
2172  adjat(0,1) = a(1,2)*a(2,0)-a(1,0)*a(2,2);
2173  adjat(1,1) = a(0,0)*a(2,2)-a(0,2)*a(2,0);
2174  adjat(2,1) = a(0,2)*a(1,0)-a(0,0)*a(1,2);
2175 
2176  adjat(0,2) = a(1,0)*a(2,1)-a(1,1)*a(2,0);
2177  adjat(1,2) = a(0,1)*a(2,0)-a(0,0)*a(2,1);
2178  adjat(2,2) = a(0,0)*a(1,1)-a(0,1)*a(1,0);
2179  }
2180 }
2181 
2183 {
2184  MFEM_ASSERT(a.Width() <= a.Height() && a.Width() >= 1 && a.Height() <= 3, "");
2185  MFEM_ASSERT(inva.Height() == a.Width(), "incorrect dimensions");
2186  MFEM_ASSERT(inva.Width() == a.Height(), "incorrect dimensions");
2187 
2188  double t;
2189 
2190  if (a.Width() < a.Height())
2191  {
2192  const double *d = a.Data();
2193  double *id = inva.Data();
2194  if (a.Height() == 2)
2195  {
2196  t = 1.0 / (d[0]*d[0] + d[1]*d[1]);
2197  id[0] = d[0] * t;
2198  id[1] = d[1] * t;
2199  }
2200  else
2201  {
2202  if (a.Width() == 1)
2203  {
2204  t = 1.0 / (d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
2205  id[0] = d[0] * t;
2206  id[1] = d[1] * t;
2207  id[2] = d[2] * t;
2208  }
2209  else
2210  {
2211  double e, g, f;
2212  e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
2213  g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
2214  f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
2215  t = 1.0 / (e*g - f*f);
2216  e *= t; g *= t; f *= t;
2217 
2218  id[0] = d[0]*g - d[3]*f;
2219  id[1] = d[3]*e - d[0]*f;
2220  id[2] = d[1]*g - d[4]*f;
2221  id[3] = d[4]*e - d[1]*f;
2222  id[4] = d[2]*g - d[5]*f;
2223  id[5] = d[5]*e - d[2]*f;
2224  }
2225  }
2226  return;
2227  }
2228 
2229 #ifdef MFEM_DEBUG
2230  t = a.Det();
2231  MFEM_ASSERT(std::abs(t) > 1.0e-14 * pow(a.FNorm()/a.Width(), a.Width()),
2232  "singular matrix!");
2233 #endif
2234 
2235  switch (a.Height())
2236  {
2237  case 1:
2238  inva(0,0) = 1.0 / a.Det();
2239  break;
2240  case 2:
2241  kernels::CalcInverse<2>(a.Data(), inva.Data());
2242  break;
2243  case 3:
2244  kernels::CalcInverse<3>(a.Data(), inva.Data());
2245  break;
2246  }
2247 }
2248 
2250 {
2251 #ifdef MFEM_DEBUG
2252  if ( (a.Width() != a.Height()) || ( (a.Height()!= 1) && (a.Height()!= 2)
2253  && (a.Height()!= 3) ) )
2254  {
2255  mfem_error("CalcInverseTranspose(...): dimension mismatch");
2256  }
2257 #endif
2258 
2259  double t = 1. / a.Det() ;
2260 
2261  switch (a.Height())
2262  {
2263  case 1:
2264  inva(0,0) = 1.0 / a(0,0);
2265  break;
2266  case 2:
2267  inva(0,0) = a(1,1) * t ;
2268  inva(1,0) = -a(0,1) * t ;
2269  inva(0,1) = -a(1,0) * t ;
2270  inva(1,1) = a(0,0) * t ;
2271  break;
2272  case 3:
2273  inva(0,0) = (a(1,1)*a(2,2)-a(1,2)*a(2,1))*t;
2274  inva(1,0) = (a(0,2)*a(2,1)-a(0,1)*a(2,2))*t;
2275  inva(2,0) = (a(0,1)*a(1,2)-a(0,2)*a(1,1))*t;
2276 
2277  inva(0,1) = (a(1,2)*a(2,0)-a(1,0)*a(2,2))*t;
2278  inva(1,1) = (a(0,0)*a(2,2)-a(0,2)*a(2,0))*t;
2279  inva(2,1) = (a(0,2)*a(1,0)-a(0,0)*a(1,2))*t;
2280 
2281  inva(0,2) = (a(1,0)*a(2,1)-a(1,1)*a(2,0))*t;
2282  inva(1,2) = (a(0,1)*a(2,0)-a(0,0)*a(2,1))*t;
2283  inva(2,2) = (a(0,0)*a(1,1)-a(0,1)*a(1,0))*t;
2284  break;
2285  }
2286 }
2287 
2288 void CalcOrtho(const DenseMatrix &J, Vector &n)
2289 {
2290  MFEM_ASSERT( ((J.Height() == 2 && J.Width() == 1)
2291  || (J.Height() == 3 && J.Width() == 2))
2292  && (J.Height() == n.Size()),
2293  "Matrix must be 3x2 or 2x1, "
2294  << "and the Vector must be sized with the rows. "
2295  << " J.Height() = " << J.Height()
2296  << ", J.Width() = " << J.Width()
2297  << ", n.Size() = " << n.Size()
2298  );
2299 
2300  const double *d = J.Data();
2301  if (J.Height() == 2)
2302  {
2303  n(0) = d[1];
2304  n(1) = -d[0];
2305  }
2306  else
2307  {
2308  n(0) = d[1]*d[5] - d[2]*d[4];
2309  n(1) = d[2]*d[3] - d[0]*d[5];
2310  n(2) = d[0]*d[4] - d[1]*d[3];
2311  }
2312 }
2313 
2314 void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
2315 {
2316  const int height = a.Height();
2317  const int width = a.Width();
2318  for (int i = 0; i < height; i++)
2319  {
2320  for (int j = 0; j <= i; j++)
2321  {
2322  double temp = 0.;
2323  for (int k = 0; k < width; k++)
2324  {
2325  temp += a(i,k) * a(j,k);
2326  }
2327  aat(j,i) = aat(i,j) = temp;
2328  }
2329  }
2330 }
2331 
2332 void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
2333 {
2334  for (int i = 0; i < A.Height(); i++)
2335  {
2336  for (int j = 0; j < i; j++)
2337  {
2338  double t = 0.;
2339  for (int k = 0; k < A.Width(); k++)
2340  {
2341  t += D(k) * A(i, k) * A(j, k);
2342  }
2343  ADAt(i, j) += t;
2344  ADAt(j, i) += t;
2345  }
2346  }
2347 
2348  // process diagonal
2349  for (int i = 0; i < A.Height(); i++)
2350  {
2351  double t = 0.;
2352  for (int k = 0; k < A.Width(); k++)
2353  {
2354  t += D(k) * A(i, k) * A(i, k);
2355  }
2356  ADAt(i, i) += t;
2357  }
2358 }
2359 
2360 void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
2361 {
2362  for (int i = 0; i < A.Height(); i++)
2363  {
2364  for (int j = 0; j <= i; j++)
2365  {
2366  double t = 0.;
2367  for (int k = 0; k < A.Width(); k++)
2368  {
2369  t += D(k) * A(i, k) * A(j, k);
2370  }
2371  ADAt(j, i) = ADAt(i, j) = t;
2372  }
2373  }
2374 }
2375 
2376 void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
2377 {
2378 #ifdef MFEM_DEBUG
2379  if (A.Height() != ABt.Height() || B.Height() != ABt.Width() ||
2380  A.Width() != B.Width())
2381  {
2382  mfem_error("MultABt(...): dimension mismatch");
2383  }
2384 #endif
2385 
2386 #ifdef MFEM_USE_LAPACK
2387  static char transa = 'N', transb = 'T';
2388  static double alpha = 1.0, beta = 0.0;
2389  int m = A.Height(), n = B.Height(), k = A.Width();
2390 
2391  dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &m,
2392  B.Data(), &n, &beta, ABt.Data(), &m);
2393 #elif 1
2394  const int ah = A.Height();
2395  const int bh = B.Height();
2396  const int aw = A.Width();
2397  const double *ad = A.Data();
2398  const double *bd = B.Data();
2399  double *cd = ABt.Data();
2400 
2401  kernels::MultABt(ah, aw, bh, ad, bd, cd);
2402 #elif 1
2403  const int ah = A.Height();
2404  const int bh = B.Height();
2405  const int aw = A.Width();
2406  const double *ad = A.Data();
2407  const double *bd = B.Data();
2408  double *cd = ABt.Data();
2409 
2410  for (int j = 0; j < bh; j++)
2411  for (int i = 0; i < ah; i++)
2412  {
2413  double d = 0.0;
2414  const double *ap = ad + i;
2415  const double *bp = bd + j;
2416  for (int k = 0; k < aw; k++)
2417  {
2418  d += (*ap) * (*bp);
2419  ap += ah;
2420  bp += bh;
2421  }
2422  *(cd++) = d;
2423  }
2424 #else
2425  int i, j, k;
2426  double d;
2427 
2428  for (i = 0; i < A.Height(); i++)
2429  for (j = 0; j < B.Height(); j++)
2430  {
2431  d = 0.0;
2432  for (k = 0; k < A.Width(); k++)
2433  {
2434  d += A(i, k) * B(j, k);
2435  }
2436  ABt(i, j) = d;
2437  }
2438 #endif
2439 }
2440 
2441 void MultADBt(const DenseMatrix &A, const Vector &D,
2442  const DenseMatrix &B, DenseMatrix &ADBt)
2443 {
2444 #ifdef MFEM_DEBUG
2445  if (A.Height() != ADBt.Height() || B.Height() != ADBt.Width() ||
2446  A.Width() != B.Width() || A.Width() != D.Size())
2447  {
2448  mfem_error("MultADBt(...): dimension mismatch");
2449  }
2450 #endif
2451 
2452  const int ah = A.Height();
2453  const int bh = B.Height();
2454  const int aw = A.Width();
2455  const double *ad = A.Data();
2456  const double *bd = B.Data();
2457  const double *dd = D.GetData();
2458  double *cd = ADBt.Data();
2459 
2460  for (int i = 0, s = ah*bh; i < s; i++)
2461  {
2462  cd[i] = 0.0;
2463  }
2464  for (int k = 0; k < aw; k++)
2465  {
2466  double *cp = cd;
2467  for (int j = 0; j < bh; j++)
2468  {
2469  const double dk_bjk = dd[k] * bd[j];
2470  for (int i = 0; i < ah; i++)
2471  {
2472  cp[i] += ad[i] * dk_bjk;
2473  }
2474  cp += ah;
2475  }
2476  ad += ah;
2477  bd += bh;
2478  }
2479 }
2480 
2481 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
2482 {
2483 #ifdef MFEM_DEBUG
2484  if (A.Height() != ABt.Height() || B.Height() != ABt.Width() ||
2485  A.Width() != B.Width())
2486  {
2487  mfem_error("AddMultABt(...): dimension mismatch");
2488  }
2489 #endif
2490 
2491 #ifdef MFEM_USE_LAPACK
2492  static char transa = 'N', transb = 'T';
2493  static double alpha = 1.0, beta = 1.0;
2494  int m = A.Height(), n = B.Height(), k = A.Width();
2495 
2496  dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &m,
2497  B.Data(), &n, &beta, ABt.Data(), &m);
2498 #elif 1
2499  const int ah = A.Height();
2500  const int bh = B.Height();
2501  const int aw = A.Width();
2502  const double *ad = A.Data();
2503  const double *bd = B.Data();
2504  double *cd = ABt.Data();
2505 
2506  for (int k = 0; k < aw; k++)
2507  {
2508  double *cp = cd;
2509  for (int j = 0; j < bh; j++)
2510  {
2511  const double bjk = bd[j];
2512  for (int i = 0; i < ah; i++)
2513  {
2514  cp[i] += ad[i] * bjk;
2515  }
2516  cp += ah;
2517  }
2518  ad += ah;
2519  bd += bh;
2520  }
2521 #else
2522  int i, j, k;
2523  double d;
2524 
2525  for (i = 0; i < A.Height(); i++)
2526  for (j = 0; j < B.Height(); j++)
2527  {
2528  d = 0.0;
2529  for (k = 0; k < A.Width(); k++)
2530  {
2531  d += A(i, k) * B(j, k);
2532  }
2533  ABt(i, j) += d;
2534  }
2535 #endif
2536 }
2537 
2538 void AddMultADBt(const DenseMatrix &A, const Vector &D,
2539  const DenseMatrix &B, DenseMatrix &ADBt)
2540 {
2541 #ifdef MFEM_DEBUG
2542  if (A.Height() != ADBt.Height() || B.Height() != ADBt.Width() ||
2543  A.Width() != B.Width() || A.Width() != D.Size())
2544  {
2545  mfem_error("AddMultADBt(...): dimension mismatch");
2546  }
2547 #endif
2548 
2549  const int ah = A.Height();
2550  const int bh = B.Height();
2551  const int aw = A.Width();
2552  const double *ad = A.Data();
2553  const double *bd = B.Data();
2554  const double *dd = D.GetData();
2555  double *cd = ADBt.Data();
2556 
2557  for (int k = 0; k < aw; k++)
2558  {
2559  double *cp = cd;
2560  for (int j = 0; j < bh; j++)
2561  {
2562  const double dk_bjk = dd[k] * bd[j];
2563  for (int i = 0; i < ah; i++)
2564  {
2565  cp[i] += ad[i] * dk_bjk;
2566  }
2567  cp += ah;
2568  }
2569  ad += ah;
2570  bd += bh;
2571  }
2572 }
2573 
2574 void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B,
2575  DenseMatrix &ABt)
2576 {
2577 #ifdef MFEM_DEBUG
2578  if (A.Height() != ABt.Height() || B.Height() != ABt.Width() ||
2579  A.Width() != B.Width())
2580  {
2581  mfem_error("AddMult_a_ABt(...): dimension mismatch");
2582  }
2583 #endif
2584 
2585 #ifdef MFEM_USE_LAPACK
2586  static char transa = 'N', transb = 'T';
2587  double alpha = a;
2588  static double beta = 1.0;
2589  int m = A.Height(), n = B.Height(), k = A.Width();
2590 
2591  dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &m,
2592  B.Data(), &n, &beta, ABt.Data(), &m);
2593 #elif 1
2594  const int ah = A.Height();
2595  const int bh = B.Height();
2596  const int aw = A.Width();
2597  const double *ad = A.Data();
2598  const double *bd = B.Data();
2599  double *cd = ABt.Data();
2600 
2601  for (int k = 0; k < aw; k++)
2602  {
2603  double *cp = cd;
2604  for (int j = 0; j < bh; j++)
2605  {
2606  const double bjk = a * bd[j];
2607  for (int i = 0; i < ah; i++)
2608  {
2609  cp[i] += ad[i] * bjk;
2610  }
2611  cp += ah;
2612  }
2613  ad += ah;
2614  bd += bh;
2615  }
2616 #else
2617  int i, j, k;
2618  double d;
2619 
2620  for (i = 0; i < A.Height(); i++)
2621  for (j = 0; j < B.Height(); j++)
2622  {
2623  d = 0.0;
2624  for (k = 0; k < A.Width(); k++)
2625  {
2626  d += A(i, k) * B(j, k);
2627  }
2628  ABt(i, j) += a * d;
2629  }
2630 #endif
2631 }
2632 
2633 void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
2634 {
2635 #ifdef MFEM_DEBUG
2636  if (A.Width() != AtB.Height() || B.Width() != AtB.Width() ||
2637  A.Height() != B.Height())
2638  {
2639  mfem_error("MultAtB(...): dimension mismatch");
2640  }
2641 #endif
2642 
2643 #ifdef MFEM_USE_LAPACK
2644  static char transa = 'T', transb = 'N';
2645  static double alpha = 1.0, beta = 0.0;
2646  int m = A.Width(), n = B.Width(), k = A.Height();
2647 
2648  dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &k,
2649  B.Data(), &k, &beta, AtB.Data(), &m);
2650 #elif 1
2651  const int ah = A.Height();
2652  const int aw = A.Width();
2653  const int bw = B.Width();
2654  const double *ad = A.Data();
2655  const double *bd = B.Data();
2656  double *cd = AtB.Data();
2657 
2658  for (int j = 0; j < bw; j++)
2659  {
2660  const double *ap = ad;
2661  for (int i = 0; i < aw; i++)
2662  {
2663  double d = 0.0;
2664  for (int k = 0; k < ah; k++)
2665  {
2666  d += ap[k] * bd[k];
2667  }
2668  *(cd++) = d;
2669  ap += ah;
2670  }
2671  bd += ah;
2672  }
2673 #else
2674  int i, j, k;
2675  double d;
2676 
2677  for (i = 0; i < A.Width(); i++)
2678  for (j = 0; j < B.Width(); j++)
2679  {
2680  d = 0.0;
2681  for (k = 0; k < A.Height(); k++)
2682  {
2683  d += A(k, i) * B(k, j);
2684  }
2685  AtB(i, j) = d;
2686  }
2687 #endif
2688 }
2689 
2690 void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
2691 {
2692  double d;
2693 
2694  for (int i = 0; i < A.Height(); i++)
2695  {
2696  for (int j = 0; j < i; j++)
2697  {
2698  d = 0.;
2699  for (int k = 0; k < A.Width(); k++)
2700  {
2701  d += A(i,k) * A(j,k);
2702  }
2703  AAt(i, j) += (d *= a);
2704  AAt(j, i) += d;
2705  }
2706  d = 0.;
2707  for (int k = 0; k < A.Width(); k++)
2708  {
2709  d += A(i,k) * A(i,k);
2710  }
2711  AAt(i, i) += a * d;
2712  }
2713 }
2714 
2715 void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
2716 {
2717  for (int i = 0; i < A.Height(); i++)
2718  {
2719  for (int j = 0; j <= i; j++)
2720  {
2721  double d = 0.;
2722  for (int k = 0; k < A.Width(); k++)
2723  {
2724  d += A(i,k) * A(j,k);
2725  }
2726  AAt(i, j) = AAt(j, i) = a * d;
2727  }
2728  }
2729 }
2730 
2731 void MultVVt(const Vector &v, DenseMatrix &vvt)
2732 {
2733  for (int i = 0; i < v.Size(); i++)
2734  {
2735  for (int j = 0; j <= i; j++)
2736  {
2737  vvt(i,j) = vvt(j,i) = v(i) * v(j);
2738  }
2739  }
2740 }
2741 
2742 void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
2743 {
2744 #ifdef MFEM_DEBUG
2745  if (v.Size() != VWt.Height() || w.Size() != VWt.Width())
2746  {
2747  mfem_error("MultVWt(...): dimension mismatch");
2748  }
2749 #endif
2750 
2751  for (int i = 0; i < v.Size(); i++)
2752  {
2753  const double vi = v(i);
2754  for (int j = 0; j < w.Size(); j++)
2755  {
2756  VWt(i, j) = vi * w(j);
2757  }
2758  }
2759 }
2760 
2761 void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
2762 {
2763  const int m = v.Size(), n = w.Size();
2764 
2765 #ifdef MFEM_DEBUG
2766  if (VWt.Height() != m || VWt.Width() != n)
2767  {
2768  mfem_error("AddMultVWt(...): dimension mismatch");
2769  }
2770 #endif
2771 
2772  for (int i = 0; i < m; i++)
2773  {
2774  const double vi = v(i);
2775  for (int j = 0; j < n; j++)
2776  {
2777  VWt(i, j) += vi * w(j);
2778  }
2779  }
2780 }
2781 
2782 void AddMultVVt(const Vector &v, DenseMatrix &VVt)
2783 {
2784  const int n = v.Size();
2785 
2786 #ifdef MFEM_DEBUG
2787  if (VVt.Height() != n || VVt.Width() != n)
2788  {
2789  mfem_error("AddMultVVt(...): dimension mismatch");
2790  }
2791 #endif
2792 
2793  for (int i = 0; i < n; i++)
2794  {
2795  const double vi = v(i);
2796  for (int j = 0; j < i; j++)
2797  {
2798  const double vivj = vi * v(j);
2799  VVt(i, j) += vivj;
2800  VVt(j, i) += vivj;
2801  }
2802  VVt(i, i) += vi * vi;
2803  }
2804 }
2805 
2806 void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
2807  DenseMatrix &VWt)
2808 {
2809  const int m = v.Size(), n = w.Size();
2810 
2811 #ifdef MFEM_DEBUG
2812  if (VWt.Height() != m || VWt.Width() != n)
2813  {
2814  mfem_error("AddMult_a_VWt(...): dimension mismatch");
2815  }
2816 #endif
2817 
2818  for (int j = 0; j < n; j++)
2819  {
2820  const double awj = a * w(j);
2821  for (int i = 0; i < m; i++)
2822  {
2823  VWt(i, j) += v(i) * awj;
2824  }
2825  }
2826 }
2827 
2828 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
2829 {
2830  MFEM_ASSERT(VVt.Height() == v.Size() && VVt.Width() == v.Size(),
2831  "incompatible dimensions!");
2832 
2833  const int n = v.Size();
2834  for (int i = 0; i < n; i++)
2835  {
2836  double avi = a * v(i);
2837  for (int j = 0; j < i; j++)
2838  {
2839  const double avivj = avi * v(j);
2840  VVt(i, j) += avivj;
2841  VVt(j, i) += avivj;
2842  }
2843  VVt(i, i) += avi * v(i);
2844  }
2845 }
2846 
2847 
2848 bool LUFactors::Factor(int m, double TOL)
2849 {
2850 #ifdef MFEM_USE_LAPACK
2851  int info = 0;
2852  if (m) { dgetrf_(&m, &m, data, &m, ipiv, &info); }
2853  return info == 0;
2854 #else
2855  // compiling without LAPACK
2856  double *data_ptr = this->data;
2857  for (int i = 0; i < m; i++)
2858  {
2859  // pivoting
2860  {
2861  int piv = i;
2862  double a = std::abs(data_ptr[piv+i*m]);
2863  for (int j = i+1; j < m; j++)
2864  {
2865  const double b = std::abs(data_ptr[j+i*m]);
2866  if (b > a)
2867  {
2868  a = b;
2869  piv = j;
2870  }
2871  }
2872  ipiv[i] = piv;
2873  if (piv != i)
2874  {
2875  // swap rows i and piv in both L and U parts
2876  for (int j = 0; j < m; j++)
2877  {
2878  mfem::Swap<double>(data_ptr[i+j*m], data_ptr[piv+j*m]);
2879  }
2880  }
2881  }
2882 
2883  if (abs(data_ptr[i + i*m]) <= TOL)
2884  {
2885  return false; // failed
2886  }
2887 
2888  const double a_ii_inv = 1.0 / data_ptr[i+i*m];
2889  for (int j = i+1; j < m; j++)
2890  {
2891  data_ptr[j+i*m] *= a_ii_inv;
2892  }
2893  for (int k = i+1; k < m; k++)
2894  {
2895  const double a_ik = data_ptr[i+k*m];
2896  for (int j = i+1; j < m; j++)
2897  {
2898  data_ptr[j+k*m] -= a_ik * data_ptr[j+i*m];
2899  }
2900  }
2901  }
2902 #endif
2903 
2904  return true; // success
2905 }
2906 
2907 double LUFactors::Det(int m) const
2908 {
2909  double det = 1.0;
2910  for (int i=0; i<m; i++)
2911  {
2912  if (ipiv[i] != i-ipiv_base)
2913  {
2914  det *= -data[m * i + i];
2915  }
2916  else
2917  {
2918  det *= data[m * i + i];
2919  }
2920  }
2921  return det;
2922 }
2923 
2924 void LUFactors::Mult(int m, int n, double *X) const
2925 {
2926  double *x = X;
2927  for (int k = 0; k < n; k++)
2928  {
2929  // X <- U X
2930  for (int i = 0; i < m; i++)
2931  {
2932  double x_i = x[i] * data[i+i*m];
2933  for (int j = i+1; j < m; j++)
2934  {
2935  x_i += x[j] * data[i+j*m];
2936  }
2937  x[i] = x_i;
2938  }
2939  // X <- L X
2940  for (int i = m-1; i >= 0; i--)
2941  {
2942  double x_i = x[i];
2943  for (int j = 0; j < i; j++)
2944  {
2945  x_i += x[j] * data[i+j*m];
2946  }
2947  x[i] = x_i;
2948  }
2949  // X <- P^{-1} X
2950  for (int i = m-1; i >= 0; i--)
2951  {
2952  mfem::Swap<double>(x[i], x[ipiv[i]-ipiv_base]);
2953  }
2954  x += m;
2955  }
2956 }
2957 
2958 void LUFactors::LSolve(int m, int n, double *X) const
2959 {
2960  double *x = X;
2961  for (int k = 0; k < n; k++)
2962  {
2963  // X <- P X
2964  for (int i = 0; i < m; i++)
2965  {
2966  mfem::Swap<double>(x[i], x[ipiv[i]-ipiv_base]);
2967  }
2968  // X <- L^{-1} X
2969  for (int j = 0; j < m; j++)
2970  {
2971  const double x_j = x[j];
2972  for (int i = j+1; i < m; i++)
2973  {
2974  x[i] -= data[i+j*m] * x_j;
2975  }
2976  }
2977  x += m;
2978  }
2979 }
2980 
2981 void LUFactors::USolve(int m, int n, double *X) const
2982 {
2983  double *x = X;
2984  // X <- U^{-1} X
2985  for (int k = 0; k < n; k++)
2986  {
2987  for (int j = m-1; j >= 0; j--)
2988  {
2989  const double x_j = ( x[j] /= data[j+j*m] );
2990  for (int i = 0; i < j; i++)
2991  {
2992  x[i] -= data[i+j*m] * x_j;
2993  }
2994  }
2995  x += m;
2996  }
2997 }
2998 
2999 void LUFactors::Solve(int m, int n, double *X) const
3000 {
3001 #ifdef MFEM_USE_LAPACK
3002  char trans = 'N';
3003  int info = 0;
3004  if (m > 0 && n > 0) { dgetrs_(&trans, &m, &n, data, &m, ipiv, X, &m, &info); }
3005  MFEM_VERIFY(!info, "LAPACK: error in DGETRS");
3006 #else
3007  // compiling without LAPACK
3008  LSolve(m, n, X);
3009  USolve(m, n, X);
3010 #endif
3011 }
3012 
3013 void LUFactors::RightSolve(int m, int n, double *X) const
3014 {
3015  double *x;
3016 #ifdef MFEM_USE_LAPACK
3017  char n_ch = 'N', side = 'R', u_ch = 'U', l_ch = 'L';
3018  double alpha = 1.0;
3019  if (m > 0 && n > 0)
3020  {
3021  dtrsm_(&side,&u_ch,&n_ch,&n_ch,&n,&m,&alpha,data,&m,X,&n);
3022  dtrsm_(&side,&l_ch,&n_ch,&u_ch,&n,&m,&alpha,data,&m,X,&n);
3023  }
3024 #else
3025  // compiling without LAPACK
3026  // X <- X U^{-1}
3027  x = X;
3028  for (int k = 0; k < n; k++)
3029  {
3030  for (int j = 0; j < m; j++)
3031  {
3032  const double x_j = ( x[j*n] /= data[j+j*m]);
3033  for (int i = j+1; i < m; i++)
3034  {
3035  x[i*n] -= data[j + i*m] * x_j;
3036  }
3037  }
3038  ++x;
3039  }
3040 
3041  // X <- X L^{-1}
3042  x = X;
3043  for (int k = 0; k < n; k++)
3044  {
3045  for (int j = m-1; j >= 0; j--)
3046  {
3047  const double x_j = x[j*n];
3048  for (int i = 0; i < j; i++)
3049  {
3050  x[i*n] -= data[j + i*m] * x_j;
3051  }
3052  }
3053  ++x;
3054  }
3055 #endif
3056  // X <- X P
3057  x = X;
3058  for (int k = 0; k < n; k++)
3059  {
3060  for (int i = m-1; i >= 0; --i)
3061  {
3062  mfem::Swap<double>(x[i*n], x[(ipiv[i]-ipiv_base)*n]);
3063  }
3064  ++x;
3065  }
3066 }
3067 
3068 void LUFactors::GetInverseMatrix(int m, double *X) const
3069 {
3070  // A^{-1} = U^{-1} L^{-1} P
3071  // X <- U^{-1} (set only the upper triangular part of X)
3072  double *x = X;
3073  for (int k = 0; k < m; k++)
3074  {
3075  const double minus_x_k = -( x[k] = 1.0/data[k+k*m] );
3076  for (int i = 0; i < k; i++)
3077  {
3078  x[i] = data[i+k*m] * minus_x_k;
3079  }
3080  for (int j = k-1; j >= 0; j--)
3081  {
3082  const double x_j = ( x[j] /= data[j+j*m] );
3083  for (int i = 0; i < j; i++)
3084  {
3085  x[i] -= data[i+j*m] * x_j;
3086  }
3087  }
3088  x += m;
3089  }
3090  // X <- X L^{-1} (use input only from the upper triangular part of X)
3091  {
3092  int k = m-1;
3093  for (int j = 0; j < k; j++)
3094  {
3095  const double minus_L_kj = -data[k+j*m];
3096  for (int i = 0; i <= j; i++)
3097  {
3098  X[i+j*m] += X[i+k*m] * minus_L_kj;
3099  }
3100  for (int i = j+1; i < m; i++)
3101  {
3102  X[i+j*m] = X[i+k*m] * minus_L_kj;
3103  }
3104  }
3105  }
3106  for (int k = m-2; k >= 0; k--)
3107  {
3108  for (int j = 0; j < k; j++)
3109  {
3110  const double L_kj = data[k+j*m];
3111  for (int i = 0; i < m; i++)
3112  {
3113  X[i+j*m] -= X[i+k*m] * L_kj;
3114  }
3115  }
3116  }
3117  // X <- X P
3118  for (int k = m-1; k >= 0; k--)
3119  {
3120  const int piv_k = ipiv[k]-ipiv_base;
3121  if (k != piv_k)
3122  {
3123  for (int i = 0; i < m; i++)
3124  {
3125  Swap<double>(X[i+k*m], X[i+piv_k*m]);
3126  }
3127  }
3128  }
3129 }
3130 
3131 void LUFactors::SubMult(int m, int n, int r, const double *A21,
3132  const double *X1, double *X2)
3133 {
3134  // X2 <- X2 - A21 X1
3135  for (int k = 0; k < r; k++)
3136  {
3137  for (int j = 0; j < m; j++)
3138  {
3139  const double x1_jk = X1[j+k*m];
3140  for (int i = 0; i < n; i++)
3141  {
3142  X2[i+k*n] -= A21[i+j*n] * x1_jk;
3143  }
3144  }
3145  }
3146 }
3147 
3149  int m, int n, double *A12, double *A21, double *A22) const
3150 {
3151  // A12 <- L^{-1} P A12
3152  LSolve(m, n, A12);
3153  // A21 <- A21 U^{-1}
3154  for (int j = 0; j < m; j++)
3155  {
3156  const double u_jj_inv = 1.0/data[j+j*m];
3157  for (int i = 0; i < n; i++)
3158  {
3159  A21[i+j*n] *= u_jj_inv;
3160  }
3161  for (int k = j+1; k < m; k++)
3162  {
3163  const double u_jk = data[j+k*m];
3164  for (int i = 0; i < n; i++)
3165  {
3166  A21[i+k*n] -= A21[i+j*n] * u_jk;
3167  }
3168  }
3169  }
3170  // A22 <- A22 - A21 A12
3171  SubMult(m, n, n, A21, A12, A22);
3172 }
3173 
3174 void LUFactors::BlockForwSolve(int m, int n, int r, const double *L21,
3175  double *B1, double *B2) const
3176 {
3177  // B1 <- L^{-1} P B1
3178  LSolve(m, r, B1);
3179  // B2 <- B2 - L21 B1
3180  SubMult(m, n, r, L21, B1, B2);
3181 }
3182 
3183 void LUFactors::BlockBackSolve(int m, int n, int r, const double *U12,
3184  const double *X2, double *Y1) const
3185 {
3186  // Y1 <- Y1 - U12 X2
3187  SubMult(n, m, r, U12, X2, Y1);
3188  // Y1 <- U^{-1} Y1
3189  USolve(m, r, Y1);
3190 }
3191 
3192 
3194  : MatrixInverse(mat)
3195 {
3196  MFEM_ASSERT(height == width, "not a square matrix");
3197  a = &mat;
3198  lu.data = new double[width*width];
3199  lu.ipiv = new int[width];
3200  Factor();
3201 }
3202 
3204  : MatrixInverse(*mat)
3205 {
3206  MFEM_ASSERT(height == width, "not a square matrix");
3207  a = mat;
3208  lu.data = new double[width*width];
3209  lu.ipiv = new int[width];
3210 }
3211 
3213 {
3214  MFEM_ASSERT(a, "DenseMatrix is not given");
3215  const double *adata = a->data;
3216  const int s = width*width;
3217  for (int i = 0; i < s; i++)
3218  {
3219  lu.data[i] = adata[i];
3220  }
3221  lu.Factor(width);
3222 }
3223 
3225 {
3226  Ainv.SetSize(width);
3227  lu.GetInverseMatrix(width, Ainv.Data());
3228 }
3229 
3231 {
3232  MFEM_VERIFY(mat.height == mat.width, "DenseMatrix is not square!");
3233  if (width != mat.width)
3234  {
3235  height = width = mat.width;
3236  delete [] lu.data;
3237  lu.data = new double[width*width];
3238  delete [] lu.ipiv;
3239  lu.ipiv = new int[width];
3240  }
3241  a = &mat;
3242  Factor();
3243 }
3244 
3246 {
3247  const DenseMatrix *p = dynamic_cast<const DenseMatrix*>(&op);
3248  MFEM_VERIFY(p != NULL, "Operator is not a DenseMatrix!");
3249  Factor(*p);
3250 }
3251 
3252 void DenseMatrixInverse::Mult(const double *x, double *y) const
3253 {
3254  for (int row = 0; row < height; row++)
3255  {
3256  y[row] = x[row];
3257  }
3258  lu.Solve(width, 1, y);
3259 }
3260 
3261 void DenseMatrixInverse::Mult(const Vector &x, Vector &y) const
3262 {
3263  y = x;
3264  lu.Solve(width, 1, y.GetData());
3265 }
3266 
3268 {
3269  X = B;
3270  lu.Solve(width, X.Width(), X.Data());
3271 }
3272 
3274 {
3275  DenseMatrix C(width);
3276  Mult(*a, C);
3277  for (int i = 0; i < width; i++)
3278  {
3279  C(i,i) -= 1.0;
3280  }
3281  mfem::out << "size = " << width << ", i_max = " << C.MaxMaxNorm() << endl;
3282 }
3283 
3285 {
3286  delete [] lu.data;
3287  delete [] lu.ipiv;
3288 }
3289 
3290 #ifdef MFEM_USE_LAPACK
3291 
3293  : mat(m)
3294 {
3295  n = mat.Width();
3296  EVal.SetSize(n);
3297  EVect.SetSize(n);
3298  ev.SetDataAndSize(NULL, n);
3299 
3300  jobz = 'V';
3301  uplo = 'U';
3302  lwork = -1;
3303  double qwork;
3304  dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
3305  &qwork, &lwork, &info);
3306 
3307  lwork = (int) qwork;
3308  work = new double[lwork];
3309 }
3310 
3312  const DenseMatrixEigensystem &other)
3313  : mat(other.mat), EVal(other.EVal), EVect(other.EVect), ev(NULL, other.n),
3314  n(other.n)
3315 {
3316  jobz = other.jobz;
3317  uplo = other.uplo;
3318  lwork = other.lwork;
3319 
3320  work = new double[lwork];
3321 }
3322 
3324 {
3325 #ifdef MFEM_DEBUG
3326  if (mat.Width() != n)
3327  {
3328  mfem_error("DenseMatrixEigensystem::Eval(): dimension mismatch");
3329  }
3330 #endif
3331 
3332  EVect = mat;
3333  dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
3334  work, &lwork, &info);
3335 
3336  if (info != 0)
3337  {
3338  mfem::err << "DenseMatrixEigensystem::Eval(): DSYEV error code: "
3339  << info << endl;
3340  mfem_error();
3341  }
3342 }
3343 
3345 {
3346  delete [] work;
3347 }
3348 
3349 
3352  bool left_eigen_vectors,
3353  bool right_eigen_vectors)
3354  : A(a), B(b)
3355 {
3356  MFEM_VERIFY(A.Height() == A.Width(), "A has to be a square matrix");
3357  MFEM_VERIFY(B.Height() == B.Width(), "B has to be a square matrix");
3358  n = A.Width();
3359  MFEM_VERIFY(B.Height() == n, "A and B dimension mismatch");
3360 
3361  jobvl = 'N';
3362  jobvr = 'N';
3363  A_copy.SetSize(n);
3364  B_copy.SetSize(n);
3365  if (left_eigen_vectors)
3366  {
3367  jobvl = 'V';
3368  Vl.SetSize(n);
3369  }
3370  if (right_eigen_vectors)
3371  {
3372  jobvr = 'V';
3373  Vr.SetSize(n);
3374  }
3375 
3376  lwork = -1;
3377  double qwork;
3378 
3379  alphar = new double[n];
3380  alphai = new double[n];
3381  beta = new double[n];
3382 
3383  int nl = max(1,Vl.Height());
3384  int nr = max(1,Vr.Height());
3385 
3386  dggev_(&jobvl,&jobvr,&n,A_copy.Data(),&n,B_copy.Data(),&n,alphar,
3387  alphai, beta, Vl.Data(), &nl, Vr.Data(), &nr,
3388  &qwork, &lwork, &info);
3389 
3390  lwork = (int) qwork;
3391  work = new double[lwork];
3392 
3393 }
3394 
3396 {
3397  int nl = max(1,Vl.Height());
3398  int nr = max(1,Vr.Height());
3399 
3400  A_copy = A;
3401  B_copy = B;
3402  dggev_(&jobvl,&jobvr,&n,A_copy.Data(),&n,B_copy.Data(),&n,alphar,
3403  alphai, beta, Vl.Data(), &nl, Vr.Data(), &nr,
3404  work, &lwork, &info);
3405 
3406  if (info != 0)
3407  {
3408  mfem::err << "DenseMatrixGeneralizedEigensystem::Eval(): DGGEV error code: "
3409  << info << endl;
3410  mfem_error();
3411  }
3412  evalues_r.SetSize(n);
3413  evalues_i.SetSize(n);
3414  for (int i = 0; i<n; i++)
3415  {
3416  if (beta[i] != 0.)
3417  {
3418  evalues_r(i) = alphar[i]/beta[i];
3419  evalues_i(i) = alphai[i]/beta[i];
3420  }
3421  else
3422  {
3423  evalues_r(i) = infinity();
3424  evalues_i(i) = infinity();
3425  }
3426  }
3427 }
3428 
3430 {
3431  delete [] alphar;
3432  delete [] alphai;
3433  delete [] beta;
3434  delete [] work;
3435 }
3436 
3438  bool left_singular_vectors,
3439  bool right_singular_vectors)
3440 {
3441  m = M.Height();
3442  n = M.Width();
3443  jobu = (left_singular_vectors)? 'S' : 'N';
3444  jobvt = (right_singular_vectors)? 'S' : 'N';
3445  Init();
3446 }
3447 
3449  bool left_singular_vectors,
3450  bool right_singular_vectors)
3451 {
3452  m = h;
3453  n = w;
3454  jobu = (left_singular_vectors)? 'S' : 'N';
3455  jobvt = (right_singular_vectors)? 'S' : 'N';
3456  Init();
3457 }
3458 
3459 void DenseMatrixSVD::Init()
3460 {
3461  sv.SetSize(min(m, n));
3462  double qwork;
3463  lwork = -1;
3464  dgesvd_(&jobu, &jobvt, &m, &n, NULL, &m, sv.GetData(), NULL, &m,
3465  NULL, &n, &qwork, &lwork, &info);
3466 
3467  lwork = (int) qwork;
3468  work = new double[lwork];
3469 }
3470 
3472 {
3473 #ifdef MFEM_DEBUG
3474  if (M.Height() != m || M.Width() != n)
3475  {
3476  mfem_error("DenseMatrixSVD::Eval()");
3477  }
3478 #endif
3479  double * datau = nullptr;
3480  double * datavt = nullptr;
3481  if (jobu == 'S')
3482  {
3483  U.SetSize(m,min(m,n));
3484  datau = U.Data();
3485  }
3486  if (jobvt == 'S')
3487  {
3488  Vt.SetSize(min(m,n),n);
3489  datavt = Vt.Data();
3490  }
3491  Mc = M;
3492  dgesvd_(&jobu, &jobvt, &m, &n, Mc.Data(), &m, sv.GetData(), datau, &m,
3493  datavt, &n, work, &lwork, &info);
3494 
3495  if (info)
3496  {
3497  mfem::err << "DenseMatrixSVD::Eval() : info = " << info << endl;
3498  mfem_error();
3499  }
3500 }
3501 
3503 {
3504  delete [] work;
3505 }
3506 
3507 #endif // if MFEM_USE_LAPACK
3508 
3509 
3510 void DenseTensor::AddMult(const Table &elem_dof, const Vector &x, Vector &y)
3511 const
3512 {
3513  int n = SizeI(), ne = SizeK();
3514  const int *I = elem_dof.GetI(), *J = elem_dof.GetJ(), *dofs;
3515  const double *d_col = tdata;
3516  double *yp = y.HostReadWrite();
3517  double x_col;
3518  const double *xp = x;
3519  // the '4' here can be tuned for given platform and compiler
3520  if (n <= 4)
3521  {
3522  for (int i = 0; i < ne; i++)
3523  {
3524  dofs = J + I[i];
3525  for (int col = 0; col < n; col++)
3526  {
3527  x_col = xp[dofs[col]];
3528  for (int row = 0; row < n; row++)
3529  {
3530  yp[dofs[row]] += x_col*d_col[row];
3531  }
3532  d_col += n;
3533  }
3534  }
3535  }
3536  else
3537  {
3538  Vector ye(n);
3539  for (int i = 0; i < ne; i++)
3540  {
3541  dofs = J + I[i];
3542  x_col = xp[dofs[0]];
3543  for (int row = 0; row < n; row++)
3544  {
3545  ye(row) = x_col*d_col[row];
3546  }
3547  d_col += n;
3548  for (int col = 1; col < n; col++)
3549  {
3550  x_col = xp[dofs[col]];
3551  for (int row = 0; row < n; row++)
3552  {
3553  ye(row) += x_col*d_col[row];
3554  }
3555  d_col += n;
3556  }
3557  for (int row = 0; row < n; row++)
3558  {
3559  yp[dofs[row]] += ye(row);
3560  }
3561  }
3562  }
3563 }
3564 
3566 {
3567  int s = SizeI() * SizeJ() * SizeK();
3568  for (int i=0; i<s; i++)
3569  {
3570  tdata[i] = c;
3571  }
3572  return *this;
3573 }
3574 
3576 {
3577  DenseTensor new_tensor(other);
3578  Swap(new_tensor);
3579  return *this;
3580 }
3581 
3582 void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const double TOL)
3583 {
3584  const int m = Mlu.SizeI();
3585  const int NE = Mlu.SizeK();
3586  P.SetSize(m*NE);
3587 
3588  auto data_all = mfem::Reshape(Mlu.ReadWrite(), m, m, NE);
3589  auto ipiv_all = mfem::Reshape(P.Write(), m, NE);
3590  Array<bool> pivot_flag(1);
3591  pivot_flag[0] = true;
3592  bool *d_pivot_flag = pivot_flag.ReadWrite();
3593 
3594  MFEM_FORALL(e, NE,
3595  {
3596  for (int i = 0; i < m; i++)
3597  {
3598  // pivoting
3599  {
3600  int piv = i;
3601  double a = fabs(data_all(piv,i,e));
3602  for (int j = i+1; j < m; j++)
3603  {
3604  const double b = fabs(data_all(j,i,e));
3605  if (b > a)
3606  {
3607  a = b;
3608  piv = j;
3609  }
3610  }
3611  ipiv_all(i,e) = piv;
3612  if (piv != i)
3613  {
3614  // swap rows i and piv in both L and U parts
3615  for (int j = 0; j < m; j++)
3616  {
3617  mfem::kernels::internal::Swap<double>(data_all(i,j,e), data_all(piv,j,e));
3618  }
3619  }
3620  } // pivot end
3621 
3622  if (abs(data_all(i,i,e)) <= TOL)
3623  {
3624  d_pivot_flag[0] = false;
3625  }
3626 
3627  const double a_ii_inv = 1.0 / data_all(i,i,e);
3628  for (int j = i+1; j < m; j++)
3629  {
3630  data_all(j,i,e) *= a_ii_inv;
3631  }
3632 
3633  for (int k = i+1; k < m; k++)
3634  {
3635  const double a_ik = data_all(i,k,e);
3636  for (int j = i+1; j < m; j++)
3637  {
3638  data_all(j,k,e) -= a_ik * data_all(j,i,e);
3639  }
3640  }
3641 
3642  } // m loop
3643 
3644  });
3645 
3646  MFEM_ASSERT(pivot_flag.HostRead()[0], "Batch LU factorization failed \n");
3647 }
3648 
3649 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X)
3650 {
3651 
3652  const int m = Mlu.SizeI();
3653  const int NE = Mlu.SizeK();
3654 
3655  auto data_all = mfem::Reshape(Mlu.Read(), m, m, NE);
3656  auto piv_all = mfem::Reshape(P.Read(), m, NE);
3657  auto x_all = mfem::Reshape(X.ReadWrite(), m, NE);
3658 
3659  MFEM_FORALL(e, NE,
3660  {
3661  kernels::LUSolve(&data_all(0, 0,e), m, &piv_all(0, e), &x_all(0,e));
3662  });
3663 
3664 }
3665 
3666 } // namespace mfem
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: densemat.cpp:1392
void dsyevr_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
Definition: densemat.cpp:835
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:2376
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
DenseMatrix & operator-=(const DenseMatrix &m)
Definition: densemat.cpp:586
void SymmetricScaling(const Vector &s)
SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
Definition: densemat.cpp:357
void SquareRootInverse()
Replaces the current matrix with its square root inverse.
Definition: densemat.cpp:728
int CheckFinite(const double *v, const int n)
Definition: vector.hpp:501
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
void BatchLUFactor(DenseTensor &Mlu, Array< int > &P, const double TOL)
Compute the LU factorization of a batch of matrices.
Definition: densemat.cpp:3582
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:2761
int * GetJ()
Definition: table.hpp:114
DenseMatrix & operator*=(double c)
Definition: densemat.cpp:599
void GetDiag(Vector &d) const
Returns the diagonal of the matrix.
Definition: densemat.cpp:1284
void SetCol(int c, const double *col)
Definition: densemat.cpp:1788
DenseTensor & operator=(double c)
Sets the tensor elements equal to constant c.
Definition: densemat.cpp:3565
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:2742
void SetRow(int r, const double *row)
Definition: densemat.cpp:1773
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition: densemat.cpp:343
MFEM_HOST_DEVICE void MultABt(const int Aheight, const int Awidth, const int Bheight, const TA *Adata, const TB *Bdata, TC *ABtdata)
Multiply a matrix of size Aheight x Awidth and data Adata with the transpose of a matrix of size Bhei...
Definition: kernels.hpp:363
void SingularValues(Vector &sv) const
Definition: densemat.cpp:1152
void dsyev_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
Definition: densemat.cpp:998
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
void Delete()
Delete the owned pointers and reset the Memory object.
double Det() const
Definition: densemat.cpp:436
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int SizeK() const
Definition: densemat.hpp:834
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
Definition: densemat.cpp:3148
MFEM_HOST_DEVICE void Add(const int height, const int width, const TALPHA alpha, const TA *Adata, const TB *Bdata, TC *Cdata)
Compute C = A + alpha*B, where the matrices A, B and C are of size height x width with data Adata...
Definition: kernels.hpp:266
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:3183
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:284
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2074
DenseMatrixSVD(DenseMatrix &M, bool left_singular_vectors=false, bool right_singlular_vectors=false)
Definition: densemat.cpp:3437
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:909
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:3510
void TestInversion()
Invert and print the numerical conditioning of the inversion.
Definition: densemat.cpp:1888
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
void CopyRows(const DenseMatrix &A, int row1, int row2)
Copy rows row1 through row2 from A to *this.
Definition: densemat.cpp:1491
double Det(int m) const
Definition: densemat.cpp:2907
void Swap(DenseTensor &t)
Definition: densemat.hpp:932
void Eval(DenseMatrix &M)
Definition: densemat.cpp:3471
Abstract data type for matrix inverse.
Definition: matrix.hpp:62
void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
Definition: densemat.cpp:2574
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3224
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3212
MFEM_HOST_DEVICE double CalcSingularvalue< 3 >(const double *data, const int i)
Return the i&#39;th singular value of the matrix of size 3 with given data.
Definition: kernels.hpp:1353
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:3068
bool IsSquare() const
Returns whether the matrix is a square matrix.
Definition: matrix.hpp:39
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2288
DenseMatrix & operator=(double c)
Sets the matrix elements equal to constant c.
Definition: densemat.cpp:540
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:520
int Capacity() const
Return the size of the allocated memory.
void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt = a * A * A^t.
Definition: densemat.cpp:2715
void dtrsm_(char *side, char *uplo, char *transa, char *diag, int *m, int *n, double *alpha, double *a, int *lda, double *b, int *ldb)
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:1072
static void SubMult(int m, int n, int r, const double *A21, const double *X1, double *X2)
Definition: densemat.cpp:3131
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: densemat.cpp:1817
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
void dggev_(char *jobvl, char *jobvr, int *n, double *a, int *lda, double *B, int *ldb, double *alphar, double *alphai, double *beta, double *vl, int *ldvl, double *vr, int *ldvr, double *work, int *lwork, int *info)
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1916
double & operator()(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.hpp:967
double Weight() const
Definition: densemat.cpp:493
void USolve(int m, int n, double *X) const
Definition: densemat.cpp:2981
double FNorm() const
Compute the Frobenius norm of the matrix.
Definition: densemat.hpp:239
MFEM_HOST_DEVICE void CalcEigenvalues< 2 >(const double *data, double *lambda, double *vec)
Definition: kernels.hpp:1072
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:188
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2146
friend class DenseMatrixInverse
Definition: densemat.hpp:26
void RightSolve(int m, int n, double *X) const
Definition: densemat.cpp:3013
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2042
double operator*(const DenseMatrix &m) const
Matrix inner product: tr(A^t B)
Definition: densemat.cpp:173
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:529
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:154
void AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
Definition: densemat.cpp:2806
double b
Definition: lissajous.cpp:42
void InvSymmetricScaling(const Vector &s)
InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
Definition: densemat.cpp:385
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1254
void BlockForwSolve(int m, int n, int r, const double *L21, double *B1, double *B2) const
Definition: densemat.cpp:3174
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:773
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:2441
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:618
bool LinearSolve(DenseMatrix &A, double *X, double TOL)
Solves the dense linear system, A * X = B for X
Definition: densemat.cpp:1938
virtual ~DenseMatrixInverse()
Destroys dense inverse matrix.
Definition: densemat.cpp:3284
void LSolve(int m, int n, double *X) const
Definition: densemat.cpp:2958
DenseMatrixGeneralizedEigensystem(DenseMatrix &a, DenseMatrix &b, bool left_eigen_vectors=false, bool right_eigen_vectors=false)
Definition: densemat.cpp:3350
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
Definition: densemat.cpp:302
double Det() const
Compute the determinant of the original DenseMatrix using the LU factors.
Definition: densemat.hpp:679
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
Definition: densemat.cpp:2782
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:304
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:1587
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:2828
void Neg()
(*this) = -(*this)
Definition: densemat.cpp:609
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: densemat.cpp:3245
MFEM_HOST_DEVICE void LUSolve(const double *data, const int m, const int *ipiv, double *x)
Assuming L.U = P.A for a factored matrix (m x m),.
Definition: kernels.hpp:1594
void Solve(int m, int n, double *X) const
Definition: densemat.cpp:2999
double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:925
void Getl1Diag(Vector &l) const
Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|.
Definition: densemat.cpp:1298
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:1708
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1270
void AddMult(const Vector &x, Vector &y) const
y += A.x
Definition: densemat.cpp:211
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:1803
int SizeI() const
Definition: densemat.hpp:832
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2182
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
void TestInversion()
Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
Definition: densemat.cpp:3273
MFEM_HOST_DEVICE void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
Matrix vector multiplication: y = A x, where the matrix A is of size height x width with given data...
Definition: kernels.hpp:163
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:786
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:630
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1359
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:2731
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:412
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2481
void dgetrs_(char *, int *, int *, double *, int *, int *, double *, int *, int *)
void BatchLUSolve(const DenseTensor &Mlu, const Array< int > &P, Vector &X)
Solve batch linear systems.
Definition: densemat.cpp:3649
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2314
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
void Swap(DenseMatrix &other)
Definition: densemat.cpp:1902
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:1648
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3252
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2249
A class to initialize the size of a Tensor.
Definition: dtensor.hpp:54
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
MFEM_HOST_DEVICE double CalcSingularvalue< 2 >(const double *data, const int i)
Return the i&#39;th singular value of the matrix of size 2 with given data.
Definition: kernels.hpp:1305
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:2360
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
virtual void PrintT(std::ostream &out=mfem::out, int width_=4) const
Prints the transpose matrix to stream out.
Definition: densemat.cpp:1862
void CopyCols(const DenseMatrix &A, int col1, int col2)
Copy columns col1 through col2 from A to *this.
Definition: densemat.cpp:1504
int SizeJ() const
Definition: densemat.hpp:833
virtual MatrixInverse * Inverse() const
Returns a pointer to the inverse matrix.
Definition: densemat.cpp:431
double a
Definition: lissajous.cpp:41
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:454
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:2538
virtual ~DenseMatrix()
Destroys dense matrix.
Definition: densemat.cpp:1909
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:1543
void AddMultTranspose(const Vector &x, Vector &y) const
y += A^t x
Definition: densemat.cpp:229
void CopyExceptMN(const DenseMatrix &A, int m, int n)
Copy All rows and columns except m and n from A.
Definition: densemat.cpp:1622
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
void Diag(double c, int n)
Creates n x n diagonal matrix with diagonal elements c.
Definition: densemat.cpp:1329
void Mult(int m, int n, double *X) const
Definition: densemat.cpp:2924
MFEM_HOST_DEVICE void CalcEigenvalues< 3 >(const double *data, double *lambda, double *vec)
Definition: kernels.hpp:1103
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Definition: densemat.cpp:2848
static const int ipiv_base
Definition: densemat.hpp:535
void GradToCurl(DenseMatrix &curl)
Definition: densemat.cpp:1417
DenseMatrixInverse()
Default constructor.
Definition: densemat.hpp:643
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1207
void GetRowSums(Vector &l) const
Compute the row sums of the DenseMatrix.
Definition: densemat.cpp:1315
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
void CalcEigenvalues(double *lambda, double *vec) const
Definition: densemat.cpp:1232
RefCoord t[3]
DenseMatrixEigensystem(DenseMatrix &m)
Definition: densemat.cpp:3292
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)
T * Write(bool on_dev=true)
Shortcut for mfem::Write(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:312
int Rank(double tol) const
Definition: densemat.cpp:1192
const double alpha
Definition: ex15.cpp:369
void AddMult_a(double a, const Vector &x, Vector &y) const
y += a * A.x
Definition: densemat.cpp:247
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
Definition: densemat.cpp:328
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:2633
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:160
void AddMultTranspose_a(double a, const Vector &x, Vector &y) const
y += a * A^t x
Definition: densemat.cpp:265
void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt += A D A^t, where D is diagonal.
Definition: densemat.cpp:2332
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:1719
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:1517
RefCoord s[3]
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
Definition: densemat.cpp:315
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
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:66
Abstract operator.
Definition: operator.hpp:24
void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
Definition: densemat.cpp:2009
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:786
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.cpp:150
void AdjustDofDirection(Array< int > &dofs)
Definition: densemat.cpp:1730
MFEM_HOST_DEVICE void Symmetrize(const int size, T *data)
Symmetrize a square matrix with given size and data: A -&gt; (A+A^T)/2.
Definition: kernels.hpp:223
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:458
void GradToDiv(Vector &div)
Definition: densemat.cpp:1476
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:2690
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
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
Prints operator in Matlab format.
Definition: densemat.cpp:1843
double * data
Definition: densemat.hpp:532
DenseMatrix & operator+=(const double *m)
Definition: densemat.cpp:573