 MFEM  v4.4.0 Finite element discretization library
operator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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
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::FormLinearSystem(const Array<int> &ess_tdof_list,
52  Vector &x, Vector &b,
53  Operator* &Aout, Vector &X, Vector &B,
54  int copy_interior)
55 {
56  const Operator *P = this->GetProlongation();
57  const Operator *R = this->GetRestriction();
58  InitTVectors(P, R, P, x, b, X, B);
59
60  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
61
62  ConstrainedOperator *constrainedA;
63  FormConstrainedSystemOperator(ess_tdof_list, constrainedA);
64  constrainedA->EliminateRHS(X, B);
65  Aout = constrainedA;
66 }
67
69  const Array<int> &trial_tdof_list,
70  const Array<int> &test_tdof_list, Vector &x, Vector &b,
71  Operator* &Aout, Vector &X, Vector &B)
72 {
73  const Operator *Pi = this->GetProlongation();
74  const Operator *Po = this->GetOutputProlongation();
75  const Operator *Ri = this->GetRestriction();
76  InitTVectors(Po, Ri, Pi, x, b, X, B);
77
78  RectangularConstrainedOperator *constrainedA;
79  FormRectangularConstrainedSystemOperator(trial_tdof_list, test_tdof_list,
80  constrainedA);
81  constrainedA->EliminateRHS(X, B);
82  Aout = constrainedA;
83 }
84
86 {
87  // Same for Rectangular and Square operators
88  const Operator *P = this->GetProlongation();
89  if (!IsIdentityProlongation(P))
90  {
91  // Apply conforming prolongation
92  x.SetSize(P->Height());
93  P->Mult(X, x);
94  }
95  else
96  {
97  // X and x point to the same data
98
99  // If the validity flags of X's Memory were changed (e.g. if it was moved
100  // to device memory) then we need to tell x about that.
101  x.SyncMemory(X);
102  }
103 }
104
106 {
107  Operator *rap;
108  if (!IsIdentityProlongation(Pi))
109  {
110  if (!IsIdentityProlongation(Po))
111  {
112  rap = new RAPOperator(*Po, *this, *Pi);
113  }
114  else
115  {
116  rap = new ProductOperator(this, Pi, false,false);
117  }
118  }
119  else
120  {
121  if (!IsIdentityProlongation(Po))
122  {
123  TransposeOperator * PoT = new TransposeOperator(Po);
124  rap = new ProductOperator(PoT, this, true,false);
125  }
126  else
127  {
128  rap = this;
129  }
130  }
131  return rap;
132 }
133
135  const Array<int> &ess_tdof_list, ConstrainedOperator* &Aout)
136 {
137  const Operator *P = this->GetProlongation();
138  Operator *rap = SetupRAP(P, P);
139
140  // Impose the boundary conditions through a ConstrainedOperator, which owns
141  // the rap operator when P and R are non-trivial
142  ConstrainedOperator *A = new ConstrainedOperator(rap, ess_tdof_list,
143  rap != this);
144  Aout = A;
145 }
146
148  const Array<int> &trial_tdof_list, const Array<int> &test_tdof_list,
150 {
151  const Operator *Pi = this->GetProlongation();
152  const Operator *Po = this->GetOutputProlongation();
153  Operator *rap = SetupRAP(Pi, Po);
154
155  // Impose the boundary conditions through a RectangularConstrainedOperator,
156  // which owns the rap operator when P and R are non-trivial
159  trial_tdof_list, test_tdof_list,
160  rap != this);
161  Aout = A;
162 }
163
164 void Operator::FormSystemOperator(const Array<int> &ess_tdof_list,
165  Operator* &Aout)
166 {
168  FormConstrainedSystemOperator(ess_tdof_list, A);
169  Aout = A;
170 }
171
173  const Array<int> &test_tdof_list,
174  Operator* &Aout)
175 {
177  FormRectangularConstrainedSystemOperator(trial_tdof_list, test_tdof_list, A);
178  Aout = A;
179 }
180
182 {
183  const Operator *Pin = this->GetProlongation();
184  const Operator *Rout = this->GetOutputRestriction();
185  Aout = new TripleProductOperator(Rout, this, Pin,false, false, false);
186 }
187
188 void Operator::PrintMatlab(std::ostream & os, int n, int m) const
189 {
190  using namespace std;
191  if (n == 0) { n = width; }
192  if (m == 0) { m = height; }
193
194  Vector x(n), y(m);
195  x = 0.0;
196
197  os << setiosflags(ios::scientific | ios::showpos);
198  for (int i = 0; i < n; i++)
199  {
200  x(i) = 1.0;
201  Mult(x, y);
202  for (int j = 0; j < m; j++)
203  {
204  if (y(j))
205  {
206  os << j+1 << " " << i+1 << " " << y(j) << '\n';
207  }
208  }
209  x(i) = 0.0;
210  }
211 }
212
213 void Operator::PrintMatlab(std::ostream &os) const
214 {
215  PrintMatlab(os, width, height);
216 }
217
218
220 {
221  mfem_error("TimeDependentOperator::ExplicitMult() is not overridden!");
222 }
223
225  Vector &) const
226 {
227  mfem_error("TimeDependentOperator::ImplicitMult() is not overridden!");
228 }
229
231 {
232  mfem_error("TimeDependentOperator::Mult() is not overridden!");
233 }
234
235 void TimeDependentOperator::ImplicitSolve(const double, const Vector &,
236  Vector &)
237 {
238  mfem_error("TimeDependentOperator::ImplicitSolve() is not overridden!");
239 }
240
242  const Vector &, const Vector &, double) const
243 {
245  "not overridden!");
246  return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
247 }
248
250 {
252  "not overridden!");
253  return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
254 }
255
257  const Vector &,
258  int, int *, double)
259 {
260  mfem_error("TimeDependentOperator::SUNImplicitSetup() is not overridden!");
261  return (-1);
262 }
263
265 {
266  mfem_error("TimeDependentOperator::SUNImplicitSolve() is not overridden!");
267  return (-1);
268 }
269
271 {
272  mfem_error("TimeDependentOperator::SUNMassSetup() is not overridden!");
273  return (-1);
274 }
275
277 {
278  mfem_error("TimeDependentOperator::SUNMassSolve() is not overridden!");
279  return (-1);
280 }
281
283 {
284  mfem_error("TimeDependentOperator::SUNMassMult() is not overridden!");
285  return (-1);
286 }
287
288
290  const Vector &dxdt,
291  Vector &y) const
292 {
293  mfem_error("SecondOrderTimeDependentOperator::Mult() is not overridden!");
294 }
295
297  const double dt1,
298  const Vector &x,
299  const Vector &dxdt,
300  Vector &k)
301 {
302  mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
303 }
304
305
307  bool ownA, bool ownB)
308  : Operator(A->Height(), B->Width()),
309  A(A), B(B), ownA(ownA), ownB(ownB), z(A->Width())
310 {
311  MFEM_VERIFY(A->Width() == B->Height(),
312  "incompatible Operators: A->Width() = " << A->Width()
313  << ", B->Height() = " << B->Height());
314
315  {
316  const Solver* SolverB = dynamic_cast<const Solver*>(B);
317  if (SolverB)
318  {
319  MFEM_VERIFY(!(SolverB->iterative_mode),
320  "Operator B of a ProductOperator should not be in iterative mode");
321  }
322  }
323 }
324
326 {
327  if (ownA) { delete A; }
328  if (ownB) { delete B; }
329 }
330
331
333  const Operator &P_)
334  : Operator(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_)
335 {
336  MFEM_VERIFY(Rt.Height() == A.Height(),
337  "incompatible Operators: Rt.Height() = " << Rt.Height()
338  << ", A.Height() = " << A.Height());
339  MFEM_VERIFY(A.Width() == P.Height(),
340  "incompatible Operators: A.Width() = " << A.Width()
341  << ", P.Height() = " << P.Height());
342
343  {
344  const Solver* SolverA = dynamic_cast<const Solver*>(&A);
345  if (SolverA)
346  {
347  MFEM_VERIFY(!(SolverA->iterative_mode),
348  "Operator A of an RAPOperator should not be in iterative mode");
349  }
350
351  const Solver* SolverP = dynamic_cast<const Solver*>(&P);
352  if (SolverP)
353  {
354  MFEM_VERIFY(!(SolverP->iterative_mode),
355  "Operator P of an RAPOperator should not be in iterative mode");
356  }
357  }
358
359  mem_class = Rt.GetMemoryClass()*P.GetMemoryClass();
360  MemoryType mem_type = GetMemoryType(A.GetMemoryClass()*mem_class);
361  Px.SetSize(P.Height(), mem_type);
362  APx.SetSize(A.Height(), mem_type);
363 }
364
365
367  const Operator *A, const Operator *B, const Operator *C,
368  bool ownA, bool ownB, bool ownC)
369  : Operator(A->Height(), C->Width())
370  , A(A), B(B), C(C)
371  , ownA(ownA), ownB(ownB), ownC(ownC)
372 {
373  MFEM_VERIFY(A->Width() == B->Height(),
374  "incompatible Operators: A->Width() = " << A->Width()
375  << ", B->Height() = " << B->Height());
376  MFEM_VERIFY(B->Width() == C->Height(),
377  "incompatible Operators: B->Width() = " << B->Width()
378  << ", C->Height() = " << C->Height());
379
380  {
381  const Solver* SolverB = dynamic_cast<const Solver*>(B);
382  if (SolverB)
383  {
384  MFEM_VERIFY(!(SolverB->iterative_mode),
385  "Operator B of a TripleProductOperator should not be in iterative mode");
386  }
387
388  const Solver* SolverC = dynamic_cast<const Solver*>(C);
389  if (SolverC)
390  {
391  MFEM_VERIFY(!(SolverC->iterative_mode),
392  "Operator C of a TripleProductOperator should not be in iterative mode");
393  }
394  }
395
396  mem_class = A->GetMemoryClass()*C->GetMemoryClass();
397  MemoryType mem_type = GetMemoryType(mem_class*B->GetMemoryClass());
398  t1.SetSize(C->Height(), mem_type);
399  t2.SetSize(B->Height(), mem_type);
400 }
401
403 {
404  if (ownA) { delete A; }
405  if (ownB) { delete B; }
406  if (ownC) { delete C; }
407 }
408
409
411  bool own_A_,
412  DiagonalPolicy diag_policy_)
413  : Operator(A->Height(), A->Width()), A(A), own_A(own_A_),
414  diag_policy(diag_policy_)
415 {
416  // 'mem_class' should work with A->Mult() and MFEM_FORALL():
418  MemoryType mem_type = GetMemoryType(mem_class);
419  list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
420  constraint_list.MakeRef(list);
421  // typically z and w are large vectors, so store them on the device
422  z.SetSize(height, mem_type); z.UseDevice(true);
423  w.SetSize(height, mem_type); w.UseDevice(true);
424 }
425
427 {
428  A->AssembleDiagonal(diag);
429
430  if (diag_policy == DIAG_KEEP) { return; }
431
432  const int csz = constraint_list.Size();
435  switch (diag_policy)
436  {
437  case DIAG_ONE:
438  MFEM_FORALL(i, csz,
439  {
440  const int id = idx[i];
441  d_diag[id] = 1.0;
442  });
443  break;
444  case DIAG_ZERO:
445  MFEM_FORALL(i, csz,
446  {
447  const int id = idx[i];
448  d_diag[id] = 0.0;
449  });
450  break;
451  default:
452  MFEM_ABORT("unknown diagonal policy");
453  break;
454  }
455 }
456
458 {
459  w = 0.0;
460  const int csz = constraint_list.Size();
463  // Use read+write access - we are modifying sub-vector of w
465  MFEM_FORALL(i, csz,
466  {
467  const int id = idx[i];
468  d_w[id] = d_x[id];
469  });
470
471  // A.AddMult(w, b, -1.0); // if available to all Operators
472  A->Mult(w, z);
473  b -= z;
474
475  // Use read+write access - we are modifying sub-vector of b
477  MFEM_FORALL(i, csz,
478  {
479  const int id = idx[i];
480  d_b[id] = d_x[id];
481  });
482 }
483
484 void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
485 {
486  const int csz = constraint_list.Size();
487  if (csz == 0)
488  {
489  A->Mult(x, y);
490  return;
491  }
492
493  z = x;
494
496  // Use read+write access - we are modifying sub-vector of z
498  MFEM_FORALL(i, csz, d_z[idx[i]] = 0.0;);
499
500  A->Mult(z, y);
501
503  // Use read+write access - we are modifying sub-vector of y
505  switch (diag_policy)
506  {
507  case DIAG_ONE:
508  MFEM_FORALL(i, csz,
509  {
510  const int id = idx[i];
511  d_y[id] = d_x[id];
512  });
513  break;
514  case DIAG_ZERO:
515  MFEM_FORALL(i, csz,
516  {
517  const int id = idx[i];
518  d_y[id] = 0.0;
519  });
520  break;
521  case DIAG_KEEP:
522  // Needs action of the operator diagonal on vector
523  mfem_error("ConstrainedOperator::Mult #1");
524  break;
525  default:
526  mfem_error("ConstrainedOperator::Mult #2");
527  break;
528  }
529 }
530
532  Operator *A,
533  const Array<int> &trial_list,
534  const Array<int> &test_list,
535  bool own_A_)
536  : Operator(A->Height(), A->Width()), A(A), own_A(own_A_)
537 {
538  // 'mem_class' should work with A->Mult() and MFEM_FORALL():
540  MemoryType mem_type = GetMemoryType(mem_class);
541  trial_list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
542  test_list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
543  trial_constraints.MakeRef(trial_list);
544  test_constraints.MakeRef(test_list);
545  // typically z and w are large vectors, so store them on the device
546  z.SetSize(height, mem_type); z.UseDevice(true);
547  w.SetSize(width, mem_type); w.UseDevice(true);
548 }
549
551  Vector &b) const
552 {
553  w = 0.0;
554  const int trial_csz = trial_constraints.Size();
557  // Use read+write access - we are modifying sub-vector of w
559  MFEM_FORALL(i, trial_csz,
560  {
561  const int id = trial_idx[i];
562  d_w[id] = d_x[id];
563  });
564
565  // A.AddMult(w, b, -1.0); // if available to all Operators
566  A->Mult(w, z);
567  b -= z;
568
569  const int test_csz = test_constraints.Size();
572  MFEM_FORALL(i, test_csz, d_b[test_idx[i]] = 0.0;);
573 }
574
576 {
577  const int trial_csz = trial_constraints.Size();
578  const int test_csz = test_constraints.Size();
579  if (trial_csz == 0)
580  {
581  A->Mult(x, y);
582  }
583  else
584  {
585  w = x;
586
588  // Use read+write access - we are modifying sub-vector of w
590  MFEM_FORALL(i, trial_csz, d_w[idx[i]] = 0.0;);
591
592  A->Mult(w, y);
593  }
594
595  if (test_csz != 0)
596  {
599  MFEM_FORALL(i, test_csz, d_y[idx[i]] = 0.0;);
600  }
601 }
602
604  Vector &y) const
605 {
606  const int trial_csz = trial_constraints.Size();
607  const int test_csz = test_constraints.Size();
608  if (test_csz == 0)
609  {
610  A->MultTranspose(x, y);
611  }
612  else
613  {
614  z = x;
615
617  // Use read+write access - we are modifying sub-vector of z
619  MFEM_FORALL(i, test_csz, d_z[idx[i]] = 0.0;);
620
621  A->MultTranspose(z, y);
622  }
623
624  if (trial_csz != 0)
625  {
628  MFEM_FORALL(i, trial_csz, d_y[idx[i]] = 0.0;);
629  }
630 }
631
633  int numSteps, double tolerance, int seed)
634 {
635  v1.SetSize(v0.Size());
636  if (seed != 0)
637  {
638  v0.Randomize(seed);
639  }
640
641  double eigenvalue = 1.0;
642
643  for (int iter = 0; iter < numSteps; ++iter)
644  {
645  double normV0;
646
647 #ifdef MFEM_USE_MPI
648  if (comm != MPI_COMM_NULL)
649  {
650  normV0 = InnerProduct(comm, v0, v0);
651  }
652  else
653  {
654  normV0 = InnerProduct(v0, v0);
655  }
656 #else
657  normV0 = InnerProduct(v0, v0);
658 #endif
659
660  v0 /= sqrt(normV0);
661  opr.Mult(v0, v1);
662
663  double eigenvalueNew;
664 #ifdef MFEM_USE_MPI
665  if (comm != MPI_COMM_NULL)
666  {
667  eigenvalueNew = InnerProduct(comm, v0, v1);
668  }
669  else
670  {
671  eigenvalueNew = InnerProduct(v0, v1);
672  }
673 #else
674  eigenvalueNew = InnerProduct(v0, v1);
675 #endif
676  double diff = std::abs((eigenvalueNew - eigenvalue) / eigenvalue);
677
678  eigenvalue = eigenvalueNew;
679  std::swap(v0, v1);
680
681  if (diff < tolerance)
682  {
683  break;
684  }
685  }
686
687  return eigenvalue;
688 }
689
690 }
static MemoryClass GetMemoryClass()
(DEPRECATED) Equivalent to GetDeviceMemoryClass().
Definition: device.hpp:285
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:289
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
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:172
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:410
virtual ~ProductOperator()
Definition: operator.cpp:325
virtual void AssembleDiagonal(Vector &diag) const
Diagonal of A, modified according to the used DiagonalPolicy.
Definition: operator.cpp:426
Array< int > constraint_list
List of constrained indices/dofs.
Definition: operator.hpp:844
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:296
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
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:235
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual void Mult(const Vector &x, Vector &y) const
Constrained operator action.
Definition: operator.cpp:484
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:249
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
virtual const Operator * GetOutputProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to output vectors for the operator...
Definition: operator.hpp:121
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate columns corresponding to &quot;essential boundary condition&quot; values specified in x from the give...
Definition: operator.cpp:550
void FormDiscreteOperator(Operator *&A)
Return in A a parallel (on truedofs) version of this rectangular operator.
Definition: operator.cpp:181
ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB)
Definition: operator.cpp:306
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:655
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:772
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
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 FormRectangularConstrainedSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, RectangularConstrainedOperator *&Aout)
see FormRectangularSystemOperator()
Definition: operator.cpp:147
virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol)
Solve the ODE linear system as setup by the method SUNImplicitSetup().
Definition: operator.cpp:264
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:219
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:68
virtual const Operator * GetProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to input vectors for the operator...
Definition: operator.hpp:115
General product operator: x -&gt; (A*B)(x) = A(B(x)).
Definition: operator.hpp:743
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:86
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
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
void FormConstrainedSystemOperator(const Array< int > &ess_tdof_list, ConstrainedOperator *&Aout)
see FormSystemOperator()
Definition: operator.cpp:134
double b
Definition: lissajous.cpp:42
virtual int SUNMassSetup()
Setup the mass matrix in the ODE system .
Definition: operator.cpp:270
const T * Read(bool on_dev=true) const
Definition: array.hpp:304
The operator x -&gt; R*A*P*x constructed through the actions of R^T, A and P.
Definition: operator.hpp:763
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:603
virtual void Mult(const Vector &x, Vector &y) const
Rectangular-constrained operator action.
Definition: operator.cpp:575
Set the diagonal value to one.
Definition: operator.hpp:49
RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
Construct the RAP operator given R^T, A and P.
Definition: operator.cpp:332
virtual const Operator * GetRestriction() const
Restriction operator from input vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:118
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
Definition: vector.hpp:241
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:695
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:689
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:276
Vector w
Auxiliary vectors.
Definition: operator.hpp:847
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
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 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:241
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:51
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition: operator.hpp:718
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:85
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
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:224
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:426
Definition: vector.hpp:454
Rectangular Operator for imposing essential boundary conditions on the input space using only the act...
Definition: operator.hpp:905
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Definition: operator.hpp:107
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate &quot;essential boundary condition&quot; values specified in x from the given right-hand side b...
Definition: operator.cpp:457
Keep the diagonal value.
Definition: operator.hpp:50
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
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:164
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:46
General triple product operator x -&gt; A*B*C*x, with ownership of the factors.
Definition: operator.hpp:807
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:55
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:256
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:230
TripleProductOperator(const Operator *A, const Operator *B, const Operator *C, bool ownA, bool ownB, bool ownC)
Definition: operator.cpp:366
virtual int SUNMassMult(const Vector &x, Vector &v)
Compute the mass matrix-vector product .
Definition: operator.cpp:282
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:188
Operator * SetupRAP(const Operator *Pi, const Operator *Po)
Returns RAP Operator of this, using input/output Prolongation matrices Pi corresponds to &quot;P&quot;...
Definition: operator.cpp:105
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:864
Vector data type.
Definition: vector.hpp:60
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:585
DiagonalPolicy diag_policy
Diagonal policy for constrained dofs.
Definition: operator.hpp:849
Base class for solvers.
Definition: operator.hpp:651
Operator * A
The unconstrained Operator.
Definition: operator.hpp:845
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:841
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:632
Abstract operator.
Definition: operator.hpp:24
virtual const double * Read(bool on_dev=true) const