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