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