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