MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
operator.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
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#include "vector.hpp"
13#include "operator.hpp"
14#include "../general/forall.hpp"
15
16#include <iostream>
17#include <iomanip>
18
19namespace mfem
20{
21
22void Operator::InitTVectors(const Operator *Po, const Operator *Ri,
23 const Operator *Pi,
24 Vector &x, Vector &b,
25 Vector &X, Vector &B) const
26{
28 {
29 // Variational restriction with Po
30 B.SetSize(Po->Width(), b);
31 Po->MultTranspose(b, B);
32 }
33 else
34 {
35 // B points to same data as b
36 B.MakeRef(b, 0, b.Size());
37 }
39 {
40 // Variational restriction with Ri
41 X.SetSize(Ri->Height(), x);
42 Ri->Mult(x, X);
43 }
44 else
45 {
46 // X points to same data as x
47 X.MakeRef(x, 0, x.Size());
48 }
49}
50
51void Operator::AddMult(const Vector &x, Vector &y, const real_t a) const
52{
53 mfem::Vector z(y.Size());
54 Mult(x, z);
56}
57
59 const real_t a) const
60{
61 mfem::Vector z(y.Size());
62 MultTranspose(x, z);
64}
65
67 Array<Vector *> &Y) const
68{
69 MFEM_ASSERT(X.Size() == Y.Size(),
70 "Number of columns mismatch in Operator::Mult!");
71 for (int i = 0; i < X.Size(); i++)
72 {
73 MFEM_ASSERT(X[i] && Y[i], "Missing Vector in Operator::Mult!");
74 Mult(*X[i], *Y[i]);
75 }
76}
77
79 Array<Vector *> &Y) const
80{
81 MFEM_ASSERT(X.Size() == Y.Size(),
82 "Number of columns mismatch in Operator::MultTranspose!");
83 for (int i = 0; i < X.Size(); i++)
84 {
85 MFEM_ASSERT(X[i] && Y[i], "Missing Vector in Operator::MultTranspose!");
86 MultTranspose(*X[i], *Y[i]);
87 }
88}
89
91 const real_t a) const
92{
93 MFEM_ASSERT(X.Size() == Y.Size(),
94 "Number of columns mismatch in Operator::AddMult!");
95 for (int i = 0; i < X.Size(); i++)
96 {
97 MFEM_ASSERT(X[i] && Y[i], "Missing Vector in Operator::AddMult!");
99 }
100}
101
103 Array<Vector *> &Y, const real_t a) const
104{
105 MFEM_ASSERT(X.Size() == Y.Size(),
106 "Number of columns mismatch in Operator::AddMultTranspose!");
107 for (int i = 0; i < X.Size(); i++)
108 {
109 MFEM_ASSERT(X[i] && Y[i], "Missing Vector in Operator::AddMultTranspose!");
111 }
112}
113
114void Operator::FormLinearSystem(const Array<int> &ess_tdof_list,
115 Vector &x, Vector &b,
116 Operator* &Aout, Vector &X, Vector &B,
117 int copy_interior)
118{
119 const Operator *P = this->GetProlongation();
120 const Operator *R = this->GetRestriction();
121 InitTVectors(P, R, P, x, b, X, B);
122
123 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
124
125 ConstrainedOperator *constrainedA;
126 FormConstrainedSystemOperator(ess_tdof_list, constrainedA);
127 constrainedA->EliminateRHS(X, B);
128 Aout = constrainedA;
129}
130
132 const Array<int> &trial_tdof_list,
133 const Array<int> &test_tdof_list, Vector &x, Vector &b,
134 Operator* &Aout, Vector &X, Vector &B)
135{
136 const Operator *Pi = this->GetProlongation();
137 const Operator *Po = this->GetOutputProlongation();
138 const Operator *Ri = this->GetRestriction();
139 InitTVectors(Po, Ri, Pi, x, b, X, B);
140
141 RectangularConstrainedOperator *constrainedA;
142 FormRectangularConstrainedSystemOperator(trial_tdof_list, test_tdof_list,
143 constrainedA);
144 constrainedA->EliminateRHS(X, B);
145 Aout = constrainedA;
146}
147
149{
150 // Same for Rectangular and Square operators
151 const Operator *P = this->GetProlongation();
153 {
154 // Apply conforming prolongation
155 x.SetSize(P->Height());
156 P->Mult(X, x);
157 }
158 else
159 {
160 // X and x point to the same data
161
162 // If the validity flags of X's Memory were changed (e.g. if it was moved
163 // to device memory) then we need to tell x about that.
164 x.SyncMemory(X);
165 }
166}
167
169{
170 Operator *rap;
171 if (!IsIdentityProlongation(Pi))
172 {
173 if (!IsIdentityProlongation(Po))
174 {
175 rap = new RAPOperator(*Po, *this, *Pi);
176 }
177 else
178 {
179 rap = new ProductOperator(this, Pi, false,false);
180 }
181 }
182 else
183 {
184 if (!IsIdentityProlongation(Po))
185 {
186 TransposeOperator * PoT = new TransposeOperator(Po);
187 rap = new ProductOperator(PoT, this, true,false);
188 }
189 else
190 {
191 rap = this;
192 }
193 }
194 return rap;
195}
196
198 const Array<int> &ess_tdof_list, ConstrainedOperator* &Aout)
199{
200 const Operator *P = this->GetProlongation();
201 Operator *rap = SetupRAP(P, P);
202
203 // Impose the boundary conditions through a ConstrainedOperator, which owns
204 // the rap operator when P and R are non-trivial
205 ConstrainedOperator *A = new ConstrainedOperator(rap, ess_tdof_list,
206 rap != this);
207 Aout = A;
208}
209
211 const Array<int> &trial_tdof_list, const Array<int> &test_tdof_list,
213{
214 const Operator *Pi = this->GetProlongation();
215 const Operator *Po = this->GetOutputProlongation();
216 Operator *rap = SetupRAP(Pi, Po);
217
218 // Impose the boundary conditions through a RectangularConstrainedOperator,
219 // which owns the rap operator when P and R are non-trivial
222 trial_tdof_list, test_tdof_list,
223 rap != this);
224 Aout = A;
225}
226
227void Operator::FormSystemOperator(const Array<int> &ess_tdof_list,
228 Operator* &Aout)
229{
231 FormConstrainedSystemOperator(ess_tdof_list, A);
232 Aout = A;
233}
234
236 const Array<int> &test_tdof_list,
237 Operator* &Aout)
238{
240 FormRectangularConstrainedSystemOperator(trial_tdof_list, test_tdof_list, A);
241 Aout = A;
242}
243
245{
246 const Operator *Pin = this->GetProlongation();
247 const Operator *Rout = this->GetOutputRestriction();
248 Aout = new TripleProductOperator(Rout, this, Pin,false, false, false);
249}
250
251void Operator::PrintMatlab(std::ostream & os, int n, int m) const
252{
253 using namespace std;
254 if (n == 0) { n = width; }
255 if (m == 0) { m = height; }
256
257 Vector x(n), y(m);
258 x = 0.0;
259
260 os << setiosflags(ios::scientific | ios::showpos);
261 for (int i = 0; i < n; i++)
262 {
263 x(i) = 1.0;
264 Mult(x, y);
265 for (int j = 0; j < m; j++)
266 {
267 if (y(j) != 0)
268 {
269 os << j+1 << " " << i+1 << " " << y(j) << '\n';
270 }
271 }
272 x(i) = 0.0;
273 }
274}
275
276void Operator::PrintMatlab(std::ostream &os) const
277{
279}
280
281
283{
284 mfem_error("TimeDependentOperator::ExplicitMult() is not overridden!");
285}
286
288 Vector &) const
289{
290 mfem_error("TimeDependentOperator::ImplicitMult() is not overridden!");
291}
292
294{
295 mfem_error("TimeDependentOperator::Mult() is not overridden!");
296}
297
299 Vector &)
300{
301 mfem_error("TimeDependentOperator::ImplicitSolve() is not overridden!");
302}
303
305 const Vector &, const Vector &, real_t) const
306{
308 "not overridden!");
309 return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
310}
311
313{
315 "not overridden!");
316 return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
317}
318
320 const Vector &,
321 int, int *, real_t)
322{
323 mfem_error("TimeDependentOperator::SUNImplicitSetup() is not overridden!");
324 return (-1);
325}
326
328{
329 mfem_error("TimeDependentOperator::SUNImplicitSolve() is not overridden!");
330 return (-1);
331}
332
334{
335 mfem_error("TimeDependentOperator::SUNMassSetup() is not overridden!");
336 return (-1);
337}
338
340{
341 mfem_error("TimeDependentOperator::SUNMassSolve() is not overridden!");
342 return (-1);
343}
344
346{
347 mfem_error("TimeDependentOperator::SUNMassMult() is not overridden!");
348 return (-1);
349}
350
351
353 const Vector &dxdt,
354 Vector &y) const
355{
356 mfem_error("SecondOrderTimeDependentOperator::Mult() is not overridden!");
357}
358
360 const real_t dt1,
361 const Vector &x,
362 const Vector &dxdt,
363 Vector &k)
364{
365 mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
366}
367
369 const Operator *B, const real_t beta,
370 bool ownA, bool ownB)
371 : Operator(A->Height(), A->Width()),
372 A(A), B(B), alpha(alpha), beta(beta), ownA(ownA), ownB(ownB),
373 z(A->Height())
374{
375 MFEM_VERIFY(A->Width() == B->Width(),
376 "incompatible Operators: different widths\n"
377 << "A->Width() = " << A->Width()
378 << ", B->Width() = " << B->Width() );
379 MFEM_VERIFY(A->Height() == B->Height(),
380 "incompatible Operators: different heights\n"
381 << "A->Height() = " << A->Height()
382 << ", B->Height() = " << B->Height() );
383
384 {
385 const Solver* SolverA = dynamic_cast<const Solver*>(A);
386 const Solver* SolverB = dynamic_cast<const Solver*>(B);
387 if (SolverA)
388 {
389 MFEM_VERIFY(!(SolverA->iterative_mode),
390 "Operator A of a SumOperator should not be in iterative mode");
391 }
392 if (SolverB)
393 {
394 MFEM_VERIFY(!(SolverB->iterative_mode),
395 "Operator B of a SumOperator should not be in iterative mode");
396 }
397 }
398
399}
400
402{
403 if (ownA) { delete A; }
404 if (ownB) { delete B; }
405}
406
408 bool ownA, bool ownB)
409 : Operator(A->Height(), B->Width()),
410 A(A), B(B), ownA(ownA), ownB(ownB), z(A->Width())
411{
412 MFEM_VERIFY(A->Width() == B->Height(),
413 "incompatible Operators: A->Width() = " << A->Width()
414 << ", B->Height() = " << B->Height());
415
416 {
417 const Solver* SolverB = dynamic_cast<const Solver*>(B);
418 if (SolverB)
419 {
420 MFEM_VERIFY(!(SolverB->iterative_mode),
421 "Operator B of a ProductOperator should not be in iterative mode");
422 }
423 }
424}
425
427{
428 if (ownA) { delete A; }
429 if (ownB) { delete B; }
430}
431
432
434 const Operator &P_)
435 : Operator(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_)
436{
437 MFEM_VERIFY(Rt.Height() == A.Height(),
438 "incompatible Operators: Rt.Height() = " << Rt.Height()
439 << ", A.Height() = " << A.Height());
440 MFEM_VERIFY(A.Width() == P.Height(),
441 "incompatible Operators: A.Width() = " << A.Width()
442 << ", P.Height() = " << P.Height());
443
444 {
445 const Solver* SolverA = dynamic_cast<const Solver*>(&A);
446 if (SolverA)
447 {
448 MFEM_VERIFY(!(SolverA->iterative_mode),
449 "Operator A of an RAPOperator should not be in iterative mode");
450 }
451
452 const Solver* SolverP = dynamic_cast<const Solver*>(&P);
453 if (SolverP)
454 {
455 MFEM_VERIFY(!(SolverP->iterative_mode),
456 "Operator P of an RAPOperator should not be in iterative mode");
457 }
458 }
459
460 mem_class = Rt.GetMemoryClass()*P.GetMemoryClass();
461 MemoryType mem_type = GetMemoryType(A.GetMemoryClass()*mem_class);
462 Px.SetSize(P.Height(), mem_type);
463 APx.SetSize(A.Height(), mem_type);
464}
465
466
468 const Operator *A, const Operator *B, const Operator *C,
469 bool ownA, bool ownB, bool ownC)
470 : Operator(A->Height(), C->Width())
471 , A(A), B(B), C(C)
472 , ownA(ownA), ownB(ownB), ownC(ownC)
473{
474 MFEM_VERIFY(A->Width() == B->Height(),
475 "incompatible Operators: A->Width() = " << A->Width()
476 << ", B->Height() = " << B->Height());
477 MFEM_VERIFY(B->Width() == C->Height(),
478 "incompatible Operators: B->Width() = " << B->Width()
479 << ", C->Height() = " << C->Height());
480
481 {
482 const Solver* SolverB = dynamic_cast<const Solver*>(B);
483 if (SolverB)
484 {
485 MFEM_VERIFY(!(SolverB->iterative_mode),
486 "Operator B of a TripleProductOperator should not be in iterative mode");
487 }
488
489 const Solver* SolverC = dynamic_cast<const Solver*>(C);
490 if (SolverC)
491 {
492 MFEM_VERIFY(!(SolverC->iterative_mode),
493 "Operator C of a TripleProductOperator should not be in iterative mode");
494 }
495 }
496
497 mem_class = A->GetMemoryClass()*C->GetMemoryClass();
498 MemoryType mem_type = GetMemoryType(mem_class*B->GetMemoryClass());
499 t1.SetSize(C->Height(), mem_type);
500 t2.SetSize(B->Height(), mem_type);
501}
502
504{
505 if (ownA) { delete A; }
506 if (ownB) { delete B; }
507 if (ownC) { delete C; }
508}
509
510
512 bool own_A_,
513 DiagonalPolicy diag_policy_)
514 : Operator(A->Height(), A->Width()), A(A), own_A(own_A_),
515 diag_policy(diag_policy_)
516{
517 // 'mem_class' should work with A->Mult() and mfem::forall():
520 list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
522 // typically z and w are large vectors, so use the device (GPU) to perform
523 // operations on them
524 z.SetSize(height, mem_type); z.UseDevice(true);
525 w.SetSize(height, mem_type); w.UseDevice(true);
526}
527
529{
530 A->AssembleDiagonal(diag);
531
532 if (diag_policy == DIAG_KEEP) { return; }
533
534 const int csz = constraint_list.Size();
537 switch (diag_policy)
538 {
539 case DIAG_ONE:
540 mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
541 {
542 const int id = idx[i];
543 d_diag[id] = 1.0;
544 });
545 break;
546 case DIAG_ZERO:
547 mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
548 {
549 const int id = idx[i];
550 d_diag[id] = 0.0;
551 });
552 break;
553 default:
554 MFEM_ABORT("unknown diagonal policy");
555 break;
556 }
557}
558
560{
561 w = 0.0;
562 const int csz = constraint_list.Size();
565 // Use read+write access - we are modifying sub-vector of w
567 mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
568 {
569 const int id = idx[i];
570 d_w[id] = d_x[id];
571 });
572
573 // A.AddMult(w, b, -1.0); // if available to all Operators
574 A->Mult(w, z);
575 b -= z;
576
577 // Use read+write access - we are modifying sub-vector of b
579 mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
580 {
581 const int id = idx[i];
582 d_b[id] = d_x[id];
583 });
584}
585
587 const bool transpose) const
588{
589 const int csz = constraint_list.Size();
590 if (csz == 0)
591 {
592 if (transpose)
593 {
594 A->MultTranspose(x, y);
595 }
596 else
597 {
598 A->Mult(x, y);
599 }
600 return;
601 }
602
603 z = x;
604
606 // Use read+write access - we are modifying sub-vector of z
608 mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i) { d_z[idx[i]] = 0.0; });
609
610 if (transpose)
611 {
612 A->MultTranspose(z, y);
613 }
614 else
615 {
616 A->Mult(z, y);
617 }
618
620 // Use read+write access - we are modifying sub-vector of y
622 switch (diag_policy)
623 {
624 case DIAG_ONE:
625 mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
626 {
627 const int id = idx[i];
628 d_y[id] = d_x[id];
629 });
630 break;
631 case DIAG_ZERO:
632 mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
633 {
634 const int id = idx[i];
635 d_y[id] = 0.0;
636 });
637 break;
638 case DIAG_KEEP:
639 // Needs action of the operator diagonal on vector
640 mfem_error("ConstrainedOperator::Mult #1");
641 break;
642 default:
643 mfem_error("ConstrainedOperator::Mult #2");
644 break;
645 }
646}
647
649{
650 constexpr bool transpose = false;
651 ConstrainedMult(x, y, transpose);
652}
653
655{
656 constexpr bool transpose = true;
657 ConstrainedMult(x, y, transpose);
658}
659
661 const real_t a) const
662{
663 Mult(x, w);
665}
666
668 Operator *A,
669 const Array<int> &trial_list,
670 const Array<int> &test_list,
671 bool own_A_)
672 : Operator(A->Height(), A->Width()), A(A), own_A(own_A_)
673{
674 // 'mem_class' should work with A->Mult() and mfem::forall():
677 trial_list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
678 test_list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
679 trial_constraints.MakeRef(trial_list);
680 test_constraints.MakeRef(test_list);
681 // typically z and w are large vectors, so store them on the device
682 z.SetSize(height, mem_type); z.UseDevice(true);
683 w.SetSize(width, mem_type); w.UseDevice(true);
684}
685
687 Vector &b) const
688{
689 w = 0.0;
690 const int trial_csz = trial_constraints.Size();
693 // Use read+write access - we are modifying sub-vector of w
695 mfem::forall(trial_csz, [=] MFEM_HOST_DEVICE (int i)
696 {
697 const int id = trial_idx[i];
698 d_w[id] = d_x[id];
699 });
700
702
703 const int test_csz = test_constraints.Size();
706 mfem::forall(test_csz, [=] MFEM_HOST_DEVICE (int i)
707 {
708 d_b[test_idx[i]] = 0.0;
709 });
710}
711
713{
714 const int trial_csz = trial_constraints.Size();
715 const int test_csz = test_constraints.Size();
716 if (trial_csz == 0)
717 {
718 A->Mult(x, y);
719 }
720 else
721 {
722 w = x;
723
725 // Use read+write access - we are modifying sub-vector of w
727 mfem::forall(trial_csz, [=] MFEM_HOST_DEVICE (int i)
728 {
729 d_w[idx[i]] = 0.0;
730 });
731
732 A->Mult(w, y);
733 }
734
735 if (test_csz != 0)
736 {
739 mfem::forall(test_csz, [=] MFEM_HOST_DEVICE (int i)
740 {
741 d_y[idx[i]] = 0.0;
742 });
743 }
744}
745
747 Vector &y) const
748{
749 const int trial_csz = trial_constraints.Size();
750 const int test_csz = test_constraints.Size();
751 if (test_csz == 0)
752 {
753 A->MultTranspose(x, y);
754 }
755 else
756 {
757 z = x;
758
760 // Use read+write access - we are modifying sub-vector of z
762 mfem::forall(test_csz, [=] MFEM_HOST_DEVICE (int i)
763 {
764 d_z[idx[i]] = 0.0;
765 });
766
767 A->MultTranspose(z, y);
768 }
769
770 if (trial_csz != 0)
771 {
774 mfem::forall(trial_csz, [=] MFEM_HOST_DEVICE (int i)
775 {
776 d_y[idx[i]] = 0.0;
777 });
778 }
779}
780
782 int numSteps, real_t tolerance, int seed)
783{
784 v1.SetSize(v0.Size());
785 if (seed != 0)
786 {
787 v0.Randomize(seed);
788 }
789
790 real_t eigenvalue = 1.0;
791
792 for (int iter = 0; iter < numSteps; ++iter)
793 {
794 real_t normV0;
795
796#ifdef MFEM_USE_MPI
797 if (comm != MPI_COMM_NULL)
798 {
799 normV0 = InnerProduct(comm, v0, v0);
800 }
801 else
802 {
803 normV0 = InnerProduct(v0, v0);
804 }
805#else
806 normV0 = InnerProduct(v0, v0);
807#endif
808
809 v0 /= sqrt(normV0);
810 opr.Mult(v0, v1);
811
812 real_t eigenvalueNew;
813#ifdef MFEM_USE_MPI
814 if (comm != MPI_COMM_NULL)
815 {
816 eigenvalueNew = InnerProduct(comm, v0, v1);
817 }
818 else
819 {
820 eigenvalueNew = InnerProduct(v0, v1);
821 }
822#else
823 eigenvalueNew = InnerProduct(v0, v1);
824#endif
825 real_t diff = std::abs((eigenvalueNew - eigenvalue) / eigenvalue);
826
827 eigenvalue = eigenvalueNew;
828 std::swap(v0, v1);
829
830 if (diff < tolerance)
831 {
832 break;
833 }
834 }
835
836 return eigenvalue;
837}
838
839}
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
Definition array.hpp:882
const T * Read(bool on_dev=true) const
Definition array.hpp:317
Square Operator for imposing essential boundary conditions using only the action, Mult(),...
Definition operator.hpp:895
Array< int > constraint_list
List of constrained indices/dofs.
Definition operator.hpp:897
void Mult(const Vector &x, Vector &y) const override
Constrained operator action.
Definition operator.cpp:648
Operator * A
The unconstrained Operator.
Definition operator.hpp:898
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate "essential boundary condition" values specified in x from the given right-hand side b.
Definition operator.cpp:559
ConstrainedOperator(Operator *A, const Array< int > &list, bool own_A=false, DiagonalPolicy diag_policy=DIAG_ONE)
Constructor from a general Operator and a list of essential indices/dofs.
Definition operator.cpp:511
void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const override
Operator application: y+=A(x) (default) or y+=a*A(x).
Definition operator.cpp:660
void AssembleDiagonal(Vector &diag) const override
Diagonal of A, modified according to the used DiagonalPolicy.
Definition operator.cpp:528
DiagonalPolicy diag_policy
Diagonal policy for constrained dofs.
Definition operator.hpp:902
void ConstrainedMult(const Vector &x, Vector &y, const bool transpose) const
Implementation of Mult or MultTranspose. TODO - Generalize to allow constraining rows and columns dif...
Definition operator.cpp:586
Vector w
Auxiliary vectors.
Definition operator.hpp:900
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition operator.cpp:654
static MemoryClass GetMemoryClass()
(DEPRECATED) Equivalent to GetDeviceMemoryClass().
Definition device.hpp:286
static MemoryClass GetDeviceMemoryClass()
Get the current Device MemoryClass. This is the MemoryClass used by most MFEM device kernels to acces...
Definition device.hpp:282
Abstract operator.
Definition operator.hpp:25
void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B)
Form a column-constrained linear system using a matrix-free approach.
Definition operator.cpp:131
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition operator.hpp:86
void FormConstrainedSystemOperator(const Array< int > &ess_tdof_list, ConstrainedOperator *&Aout)
see FormSystemOperator()
Definition operator.cpp:197
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B, int copy_interior=0)
Form a constrained linear system using a matrix-free approach.
Definition operator.cpp:114
int width
Dimension of the input / number of columns in the matrix.
Definition operator.hpp:28
void FormSystemOperator(const Array< int > &ess_tdof_list, Operator *&A)
Return in A a parallel (on truedofs) version of this square operator.
Definition operator.cpp:227
void FormDiscreteOperator(Operator *&A)
Return in A a parallel (on truedofs) version of this rectangular operator.
Definition operator.cpp:244
virtual void ArrayMultTranspose(const Array< const Vector * > &X, Array< Vector * > &Y) const
Action of the transpose operator on a matrix: Y=A^t(X).
Definition operator.cpp:78
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
virtual void ArrayAddMult(const Array< const Vector * > &X, Array< Vector * > &Y, const real_t a=1.0) const
Operator application on a matrix: Y+=A(X) (default) or Y+=a*A(X).
Definition operator.cpp:90
virtual const Operator * GetRestriction() const
Restriction operator from input vectors for the operator to linear algebra (linear system) vectors....
Definition operator.hpp:143
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition operator.hpp:48
@ DIAG_ONE
Set the diagonal value to one.
Definition operator.hpp:50
@ DIAG_KEEP
Keep the diagonal value.
Definition operator.hpp:51
@ DIAG_ZERO
Set the diagonal value to zero.
Definition operator.hpp:49
virtual const Operator * GetOutputRestriction() const
Restriction operator from output vectors for the operator to linear algebra (linear system) vectors....
Definition operator.hpp:160
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
virtual const Operator * GetOutputProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to output vectors for the operator...
Definition operator.hpp:147
virtual void ArrayAddMultTranspose(const Array< const Vector * > &X, Array< Vector * > &Y, const real_t a=1.0) const
Operator transpose application on a matrix: Y+=A^t(X) (default) or Y+=a*A^t(X).
Definition operator.cpp:102
virtual void ArrayMult(const Array< const Vector * > &X, Array< Vector * > &Y) const
Operator application on a matrix: Y=A(X).
Definition operator.cpp:66
virtual const Operator * GetProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to input vectors for the operator....
Definition operator.hpp:139
void FormRectangularSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Operator *&A)
Return in A a parallel (on truedofs) version of this rectangular operator (including constraints).
Definition operator.cpp:235
Operator * SetupRAP(const Operator *Pi, const Operator *Po)
Returns RAP Operator of this, using input/output Prolongation matrices Pi corresponds to "P",...
Definition operator.cpp:168
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition operator.cpp:148
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Definition operator.hpp:131
void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi, Vector &x, Vector &b, Vector &X, Vector &B) const
Initializes memory for true vectors of linear system.
Definition operator.cpp:22
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
Operator transpose application: y+=A^t(x) (default) or y+=a*A^t(x).
Definition operator.cpp:58
void FormRectangularConstrainedSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, RectangularConstrainedOperator *&Aout)
see FormRectangularSystemOperator()
Definition operator.cpp:210
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition operator.hpp:93
void PrintMatlab(std::ostream &out, int n, int m=0) const
Prints operator with input size n and output size m in Matlab format.
Definition operator.cpp:251
virtual void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const
Operator application: y+=A(x) (default) or y+=a*A(x).
Definition operator.cpp:51
real_t EstimateLargestEigenvalue(Operator &opr, Vector &v0, int numSteps=10, real_t tolerance=1e-8, int seed=12345)
Returns an estimate of the largest eigenvalue of the operator opr using the iterative power method.
Definition operator.cpp:781
General product operator: x -> (A*B)(x) = A(B(x)).
Definition operator.hpp:797
ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB)
Definition operator.cpp:407
virtual ~ProductOperator()
Definition operator.cpp:426
The operator x -> R*A*P*x constructed through the actions of R^T, A and P.
Definition operator.hpp:817
RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
Construct the RAP operator given R^T, A and P.
Definition operator.cpp:433
Rectangular Operator for imposing essential boundary conditions on the input space using only the act...
Definition operator.hpp:970
RectangularConstrainedOperator(Operator *A, const Array< int > &trial_list, const Array< int > &test_list, bool own_A=false)
Constructor from a general Operator and a list of essential indices/dofs.
Definition operator.cpp:667
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition operator.cpp:746
virtual void Mult(const Vector &x, Vector &y) const
Rectangular-constrained operator action.
Definition operator.cpp:712
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate columns corresponding to "essential boundary condition" values specified in x from the give...
Definition operator.cpp:686
virtual void Mult(const Vector &x, const Vector &dxdt, Vector &y) const
Perform the action of the operator: y = k = f(x,@ dxdt, t), where k solves the algebraic equation F(x...
Definition operator.cpp:352
virtual void ImplicitSolve(const real_t fac0, const real_t fac1, const Vector &x, const Vector &dxdt, Vector &k)
Solve the equation: k = f(x + fac0 k, dxdt + fac1 k, t), for the unknown k at the current time t.
Definition operator.cpp:359
Base class for solvers.
Definition operator.hpp:683
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
virtual ~SumOperator()
Definition operator.cpp:401
SumOperator(const Operator *A, const real_t alpha, const Operator *B, const real_t beta, bool ownA, bool ownB)
Definition operator.cpp:368
virtual Operator & GetExplicitGradient(const Vector &x) const
Return an Operator representing dG/dx at the given point x and the currently set time.
Definition operator.cpp:312
virtual int SUNImplicitSolve(const Vector &b, Vector &x, real_t tol)
Solve the ODE linear system as setup by the method SUNImplicitSetup().
Definition operator.cpp:327
virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const
Perform the action of the implicit part of the operator, F: y = F(x, k, t) where t is the current tim...
Definition operator.cpp:287
virtual int SUNMassMult(const Vector &x, Vector &v)
Compute the mass matrix-vector product .
Definition operator.cpp:345
virtual int SUNMassSetup()
Setup the mass matrix in the ODE system .
Definition operator.cpp:333
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x,...
Definition operator.cpp:293
virtual void ExplicitMult(const Vector &x, Vector &y) const
Perform the action of the explicit part of the operator, G: y = G(x, t) where t is the current time.
Definition operator.cpp:282
virtual int SUNMassSolve(const Vector &b, Vector &x, real_t tol)
Solve the mass matrix linear system as setup by the method SUNMassSetup().
Definition operator.cpp:339
virtual void ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Definition operator.cpp:298
virtual int SUNImplicitSetup(const Vector &x, const Vector &fx, int jok, int *jcur, real_t gamma)
Setup the ODE linear system or , where .
Definition operator.cpp:319
virtual Operator & GetImplicitGradient(const Vector &x, const Vector &k, real_t shift) const
Return an Operator representing (dF/dk shift + dF/dx) at the given x, k, and the currently set time.
Definition operator.cpp:304
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition operator.hpp:750
General triple product operator x -> A*B*C*x, with ownership of the factors.
Definition operator.hpp:861
TripleProductOperator(const Operator *A, const Operator *B, const Operator *C, bool ownA, bool ownB, bool ownC)
Definition operator.cpp:467
Vector data type.
Definition vector.hpp:80
void Randomize(int seed=0)
Set random values in the vector.
Definition vector.cpp:816
virtual const real_t * Read(bool on_dev=true) const
Definition vector.hpp:474
Definition vector.hpp:490
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
Definition vector.hpp:256
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:136
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
void SetSubVectorComplement(const Array< int > &dofs, const real_t val)
Set all vector entries NOT in the dofs Array to the given val.
Definition vector.cpp:739
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:247
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
Vector beta
const real_t alpha
Definition ex15.cpp:369
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
void mfem_error(const char *msg)
Definition error.cpp:154
real_t InnerProduct(HypreParVector *x, HypreParVector *y)
Definition hypre.cpp:439
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
bool IsIdentityProlongation(const Operator *P)
Definition operator.hpp:720
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.
void forall(int N, lambda &&body)
Definition forall.hpp:754