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