MFEM  v4.6.0
Finite element discretization library
operator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #include "vector.hpp"
13 #include "operator.hpp"
14 #include "../general/forall.hpp"
15 
16 #include <iostream>
17 #include <iomanip>
18 
19 namespace mfem
20 {
21 
22 void Operator::InitTVectors(const Operator *Po, const Operator *Ri,
23  const Operator *Pi,
24  Vector &x, Vector &b,
25  Vector &X, Vector &B) const
26 {
27  if (!IsIdentityProlongation(Po))
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  }
38  if (!IsIdentityProlongation(Pi))
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 
51 void Operator::AddMult(const Vector &x, Vector &y, const double a) const
52 {
53  mfem::Vector z(y.Size());
54  Mult(x, z);
55  y.Add(a, z);
56 }
57 
59  const double a) const
60 {
61  mfem::Vector z(y.Size());
62  MultTranspose(x, z);
63  y.Add(a, 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 double 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!");
98  AddMult(*X[i], *Y[i], a);
99  }
100 }
101 
103  Array<Vector *> &Y, const double 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!");
110  AddMultTranspose(*X[i], *Y[i], a);
111  }
112 }
113 
114 void 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();
152  if (!IsIdentityProlongation(P))
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 
227 void 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 
251 void 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 
276 void Operator::PrintMatlab(std::ostream &os) const
277 {
278  PrintMatlab(os, width, height);
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 
298 void TimeDependentOperator::ImplicitSolve(const double, const Vector &,
299  Vector &)
300 {
301  mfem_error("TimeDependentOperator::ImplicitSolve() is not overridden!");
302 }
303 
305  const Vector &, const Vector &, double) const
306 {
307  mfem_error("TimeDependentOperator::GetImplicitGradient() is "
308  "not overridden!");
309  return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
310 }
311 
313 {
314  mfem_error("TimeDependentOperator::GetExplicitGradient() is "
315  "not overridden!");
316  return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
317 }
318 
320  const Vector &,
321  int, int *, double)
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 double dt1,
361  const Vector &x,
362  const Vector &dxdt,
363  Vector &k)
364 {
365  mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
366 }
367 
368 
370  bool ownA, bool ownB)
371  : Operator(A->Height(), B->Width()),
372  A(A), B(B), ownA(ownA), ownB(ownB), z(A->Width())
373 {
374  MFEM_VERIFY(A->Width() == B->Height(),
375  "incompatible Operators: A->Width() = " << A->Width()
376  << ", B->Height() = " << B->Height());
377 
378  {
379  const Solver* SolverB = dynamic_cast<const Solver*>(B);
380  if (SolverB)
381  {
382  MFEM_VERIFY(!(SolverB->iterative_mode),
383  "Operator B of a ProductOperator should not be in iterative mode");
384  }
385  }
386 }
387 
389 {
390  if (ownA) { delete A; }
391  if (ownB) { delete B; }
392 }
393 
394 
396  const Operator &P_)
397  : Operator(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_)
398 {
399  MFEM_VERIFY(Rt.Height() == A.Height(),
400  "incompatible Operators: Rt.Height() = " << Rt.Height()
401  << ", A.Height() = " << A.Height());
402  MFEM_VERIFY(A.Width() == P.Height(),
403  "incompatible Operators: A.Width() = " << A.Width()
404  << ", P.Height() = " << P.Height());
405 
406  {
407  const Solver* SolverA = dynamic_cast<const Solver*>(&A);
408  if (SolverA)
409  {
410  MFEM_VERIFY(!(SolverA->iterative_mode),
411  "Operator A of an RAPOperator should not be in iterative mode");
412  }
413 
414  const Solver* SolverP = dynamic_cast<const Solver*>(&P);
415  if (SolverP)
416  {
417  MFEM_VERIFY(!(SolverP->iterative_mode),
418  "Operator P of an RAPOperator should not be in iterative mode");
419  }
420  }
421 
422  mem_class = Rt.GetMemoryClass()*P.GetMemoryClass();
423  MemoryType mem_type = GetMemoryType(A.GetMemoryClass()*mem_class);
424  Px.SetSize(P.Height(), mem_type);
425  APx.SetSize(A.Height(), mem_type);
426 }
427 
428 
430  const Operator *A, const Operator *B, const Operator *C,
431  bool ownA, bool ownB, bool ownC)
432  : Operator(A->Height(), C->Width())
433  , A(A), B(B), C(C)
434  , ownA(ownA), ownB(ownB), ownC(ownC)
435 {
436  MFEM_VERIFY(A->Width() == B->Height(),
437  "incompatible Operators: A->Width() = " << A->Width()
438  << ", B->Height() = " << B->Height());
439  MFEM_VERIFY(B->Width() == C->Height(),
440  "incompatible Operators: B->Width() = " << B->Width()
441  << ", C->Height() = " << C->Height());
442 
443  {
444  const Solver* SolverB = dynamic_cast<const Solver*>(B);
445  if (SolverB)
446  {
447  MFEM_VERIFY(!(SolverB->iterative_mode),
448  "Operator B of a TripleProductOperator should not be in iterative mode");
449  }
450 
451  const Solver* SolverC = dynamic_cast<const Solver*>(C);
452  if (SolverC)
453  {
454  MFEM_VERIFY(!(SolverC->iterative_mode),
455  "Operator C of a TripleProductOperator should not be in iterative mode");
456  }
457  }
458 
459  mem_class = A->GetMemoryClass()*C->GetMemoryClass();
460  MemoryType mem_type = GetMemoryType(mem_class*B->GetMemoryClass());
461  t1.SetSize(C->Height(), mem_type);
462  t2.SetSize(B->Height(), mem_type);
463 }
464 
466 {
467  if (ownA) { delete A; }
468  if (ownB) { delete B; }
469  if (ownC) { delete C; }
470 }
471 
472 
474  bool own_A_,
475  DiagonalPolicy diag_policy_)
476  : Operator(A->Height(), A->Width()), A(A), own_A(own_A_),
477  diag_policy(diag_policy_)
478 {
479  // 'mem_class' should work with A->Mult() and mfem::forall():
481  MemoryType mem_type = GetMemoryType(mem_class);
482  list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
483  constraint_list.MakeRef(list);
484  // typically z and w are large vectors, so store them on the device
485  z.SetSize(height, mem_type); z.UseDevice(true);
486  w.SetSize(height, mem_type); w.UseDevice(true);
487 }
488 
490 {
491  A->AssembleDiagonal(diag);
492 
493  if (diag_policy == DIAG_KEEP) { return; }
494 
495  const int csz = constraint_list.Size();
496  auto d_diag = diag.ReadWrite();
497  auto idx = constraint_list.Read();
498  switch (diag_policy)
499  {
500  case DIAG_ONE:
501  mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
502  {
503  const int id = idx[i];
504  d_diag[id] = 1.0;
505  });
506  break;
507  case DIAG_ZERO:
508  mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
509  {
510  const int id = idx[i];
511  d_diag[id] = 0.0;
512  });
513  break;
514  default:
515  MFEM_ABORT("unknown diagonal policy");
516  break;
517  }
518 }
519 
521 {
522  w = 0.0;
523  const int csz = constraint_list.Size();
524  auto idx = constraint_list.Read();
525  auto d_x = x.Read();
526  // Use read+write access - we are modifying sub-vector of w
527  auto d_w = w.ReadWrite();
528  mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
529  {
530  const int id = idx[i];
531  d_w[id] = d_x[id];
532  });
533 
534  // A.AddMult(w, b, -1.0); // if available to all Operators
535  A->Mult(w, z);
536  b -= z;
537 
538  // Use read+write access - we are modifying sub-vector of b
539  auto d_b = b.ReadWrite();
540  mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
541  {
542  const int id = idx[i];
543  d_b[id] = d_x[id];
544  });
545 }
546 
547 void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
548 {
549  const int csz = constraint_list.Size();
550  if (csz == 0)
551  {
552  A->Mult(x, y);
553  return;
554  }
555 
556  z = x;
557 
558  auto idx = constraint_list.Read();
559  // Use read+write access - we are modifying sub-vector of z
560  auto d_z = z.ReadWrite();
561  mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i) { d_z[idx[i]] = 0.0; });
562 
563  A->Mult(z, y);
564 
565  auto d_x = x.Read();
566  // Use read+write access - we are modifying sub-vector of y
567  auto d_y = y.ReadWrite();
568  switch (diag_policy)
569  {
570  case DIAG_ONE:
571  mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
572  {
573  const int id = idx[i];
574  d_y[id] = d_x[id];
575  });
576  break;
577  case DIAG_ZERO:
578  mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
579  {
580  const int id = idx[i];
581  d_y[id] = 0.0;
582  });
583  break;
584  case DIAG_KEEP:
585  // Needs action of the operator diagonal on vector
586  mfem_error("ConstrainedOperator::Mult #1");
587  break;
588  default:
589  mfem_error("ConstrainedOperator::Mult #2");
590  break;
591  }
592 }
593 
595  Operator *A,
596  const Array<int> &trial_list,
597  const Array<int> &test_list,
598  bool own_A_)
599  : Operator(A->Height(), A->Width()), A(A), own_A(own_A_)
600 {
601  // 'mem_class' should work with A->Mult() and mfem::forall():
603  MemoryType mem_type = GetMemoryType(mem_class);
604  trial_list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
605  test_list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
606  trial_constraints.MakeRef(trial_list);
607  test_constraints.MakeRef(test_list);
608  // typically z and w are large vectors, so store them on the device
609  z.SetSize(height, mem_type); z.UseDevice(true);
610  w.SetSize(width, mem_type); w.UseDevice(true);
611 }
612 
614  Vector &b) const
615 {
616  w = 0.0;
617  const int trial_csz = trial_constraints.Size();
618  auto trial_idx = trial_constraints.Read();
619  auto d_x = x.Read();
620  // Use read+write access - we are modifying sub-vector of w
621  auto d_w = w.ReadWrite();
622  mfem::forall(trial_csz, [=] MFEM_HOST_DEVICE (int i)
623  {
624  const int id = trial_idx[i];
625  d_w[id] = d_x[id];
626  });
627 
628  // A.AddMult(w, b, -1.0); // if available to all Operators
629  A->Mult(w, z);
630  b -= z;
631 
632  const int test_csz = test_constraints.Size();
633  auto test_idx = test_constraints.Read();
634  auto d_b = b.ReadWrite();
635  mfem::forall(test_csz, [=] MFEM_HOST_DEVICE (int i)
636  {
637  d_b[test_idx[i]] = 0.0;
638  });
639 }
640 
642 {
643  const int trial_csz = trial_constraints.Size();
644  const int test_csz = test_constraints.Size();
645  if (trial_csz == 0)
646  {
647  A->Mult(x, y);
648  }
649  else
650  {
651  w = x;
652 
653  auto idx = trial_constraints.Read();
654  // Use read+write access - we are modifying sub-vector of w
655  auto d_w = w.ReadWrite();
656  mfem::forall(trial_csz, [=] MFEM_HOST_DEVICE (int i)
657  {
658  d_w[idx[i]] = 0.0;
659  });
660 
661  A->Mult(w, y);
662  }
663 
664  if (test_csz != 0)
665  {
666  auto idx = test_constraints.Read();
667  auto d_y = y.ReadWrite();
668  mfem::forall(test_csz, [=] MFEM_HOST_DEVICE (int i)
669  {
670  d_y[idx[i]] = 0.0;
671  });
672  }
673 }
674 
676  Vector &y) const
677 {
678  const int trial_csz = trial_constraints.Size();
679  const int test_csz = test_constraints.Size();
680  if (test_csz == 0)
681  {
682  A->MultTranspose(x, y);
683  }
684  else
685  {
686  z = x;
687 
688  auto idx = test_constraints.Read();
689  // Use read+write access - we are modifying sub-vector of z
690  auto d_z = z.ReadWrite();
691  mfem::forall(test_csz, [=] MFEM_HOST_DEVICE (int i)
692  {
693  d_z[idx[i]] = 0.0;
694  });
695 
696  A->MultTranspose(z, y);
697  }
698 
699  if (trial_csz != 0)
700  {
701  auto idx = trial_constraints.Read();
702  auto d_y = y.ReadWrite();
703  mfem::forall(trial_csz, [=] MFEM_HOST_DEVICE (int i)
704  {
705  d_y[idx[i]] = 0.0;
706  });
707  }
708 }
709 
711  int numSteps, double tolerance, int seed)
712 {
713  v1.SetSize(v0.Size());
714  if (seed != 0)
715  {
716  v0.Randomize(seed);
717  }
718 
719  double eigenvalue = 1.0;
720 
721  for (int iter = 0; iter < numSteps; ++iter)
722  {
723  double normV0;
724 
725 #ifdef MFEM_USE_MPI
726  if (comm != MPI_COMM_NULL)
727  {
728  normV0 = InnerProduct(comm, v0, v0);
729  }
730  else
731  {
732  normV0 = InnerProduct(v0, v0);
733  }
734 #else
735  normV0 = InnerProduct(v0, v0);
736 #endif
737 
738  v0 /= sqrt(normV0);
739  opr.Mult(v0, v1);
740 
741  double eigenvalueNew;
742 #ifdef MFEM_USE_MPI
743  if (comm != MPI_COMM_NULL)
744  {
745  eigenvalueNew = InnerProduct(comm, v0, v1);
746  }
747  else
748  {
749  eigenvalueNew = InnerProduct(v0, v1);
750  }
751 #else
752  eigenvalueNew = InnerProduct(v0, v1);
753 #endif
754  double diff = std::abs((eigenvalueNew - eigenvalue) / eigenvalue);
755 
756  eigenvalue = eigenvalueNew;
757  std::swap(v0, v1);
758 
759  if (diff < tolerance)
760  {
761  break;
762  }
763  }
764 
765  return eigenvalue;
766 }
767 
768 }
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
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
static MemoryClass GetMemoryClass()
(DEPRECATED) Equivalent to GetDeviceMemoryClass().
Definition: device.hpp:285
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
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:473
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
Operator transpose application: y+=A^t(x) (default) or y+=a*A^t(x).
Definition: operator.cpp:58
virtual Operator & GetImplicitGradient(const Vector &x, const Vector &k, double 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
virtual ~ProductOperator()
Definition: operator.cpp:388
Array< int > constraint_list
List of constrained indices/dofs.
Definition: operator.hpp:875
virtual void ImplicitSolve(const double fac0, const double 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
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:613
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
Definition: vector.hpp:235
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Operator application: y+=A(x) (default) or y+=a*A(x).
Definition: operator.cpp:51
virtual void ImplicitSolve(const double 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 void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
virtual const Operator * GetRestriction() const
Restriction operator from input vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:143
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
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 MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:86
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
void FormDiscreteOperator(Operator *&A)
Return in A a parallel (on truedofs) version of this rectangular operator.
Definition: operator.cpp:244
ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB)
Definition: operator.cpp:369
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
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:675
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:686
STL namespace.
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:817
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
void FormRectangularConstrainedSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, RectangularConstrainedOperator *&Aout)
see FormRectangularSystemOperator()
Definition: operator.cpp:210
virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol)
Solve the ODE linear system as setup by the method SUNImplicitSetup().
Definition: operator.cpp:327
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 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
General product operator: x -> (A*B)(x) = A(B(x)).
Definition: operator.hpp:774
virtual void ArrayAddMultTranspose(const Array< const Vector *> &X, Array< Vector *> &Y, const double 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
static MemoryClass GetDeviceMemoryClass()
Get the current Device MemoryClass. This is the MemoryClass used by most MFEM device kernels to acces...
Definition: device.hpp:281
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
virtual void ArrayMult(const Array< const Vector *> &X, Array< Vector *> &Y) const
Operator application on a matrix: Y=A(X).
Definition: operator.cpp:66
void FormConstrainedSystemOperator(const Array< int > &ess_tdof_list, ConstrainedOperator *&Aout)
see FormSystemOperator()
Definition: operator.cpp:197
double b
Definition: lissajous.cpp:42
virtual int SUNMassSetup()
Setup the mass matrix in the ODE system .
Definition: operator.cpp:333
The operator x -> R*A*P*x constructed through the actions of R^T, A and P.
Definition: operator.hpp:794
virtual void Mult(const Vector &x, Vector &y) const
Rectangular-constrained operator action.
Definition: operator.cpp:641
Set the diagonal value to one.
Definition: operator.hpp:50
RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
Construct the RAP operator given R^T, A and P.
Definition: operator.cpp:395
virtual const Operator * GetProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to input vectors for the operator...
Definition: operator.hpp:139
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the dofs Array to the given val.
Definition: vector.cpp:740
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:720
virtual const Operator * GetOutputRestriction() const
Restriction operator from output vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:160
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, double tol)
Solve the mass matrix linear system as setup by the method SUNMassSetup().
Definition: operator.cpp:339
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:520
void forall(int N, lambda &&body)
Definition: forall.hpp:742
Vector w
Auxiliary vectors.
Definition: operator.hpp:878
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
virtual void AssembleDiagonal(Vector &diag) const
Diagonal of A, modified according to the used DiagonalPolicy.
Definition: operator.cpp:489
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 Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition: operator.hpp:749
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
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
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
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
double a
Definition: lissajous.cpp:41
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:414
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:469
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
Rectangular Operator for imposing essential boundary conditions on the input space using only the act...
Definition: operator.hpp:938
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
Keep the diagonal value.
Definition: operator.hpp:51
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
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:47
General triple product operator x -> A*B*C*x, with ownership of the factors.
Definition: operator.hpp:838
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:60
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
virtual int SUNImplicitSetup(const Vector &x, const Vector &fx, int jok, int *jcur, double gamma)
Setup the ODE linear system or , where .
Definition: operator.cpp:319
TripleProductOperator(const Operator *A, const Operator *B, const Operator *C, bool ownA, bool ownB, bool ownC)
Definition: operator.cpp:429
virtual int SUNMassMult(const Vector &x, Vector &v)
Compute the mass matrix-vector product .
Definition: operator.cpp:345
virtual void ArrayAddMult(const Array< const Vector *> &X, Array< Vector *> &Y, const double a=1.0) const
Operator application on a matrix: Y+=A(X) (default) or Y+=a*A(X).
Definition: operator.cpp:90
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
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:872
Vector data type.
Definition: vector.hpp:58
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:581
virtual const Operator * GetOutputProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to output vectors for the operator...
Definition: operator.hpp:147
DiagonalPolicy diag_policy
Diagonal policy for constrained dofs.
Definition: operator.hpp:880
virtual void Mult(const Vector &x, Vector &y) const
Constrained operator action.
Definition: operator.cpp:547
Base class for solvers.
Definition: operator.hpp:682
Operator * A
The unconstrained Operator.
Definition: operator.hpp:876
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:872
double EstimateLargestEigenvalue(Operator &opr, Vector &v0, int numSteps=10, double 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:710
Abstract operator.
Definition: operator.hpp:24
Set the diagonal value to zero.
Definition: operator.hpp:49
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
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
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:594