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