1 // Copyright (c) 2010-2021, 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
InitTVectors(const Operator * Po,const Operator * Ri,const Operator * Pi,Vector & x,Vector & b,Vector & X,Vector & B) const22 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
FormLinearSystem(const Array<int> & ess_tdof_list,Vector & x,Vector & b,Operator * & Aout,Vector & X,Vector & B,int copy_interior)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
FormRectangularLinearSystem(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,Vector & x,Vector & b,Operator * & Aout,Vector & X,Vector & B)68 void Operator::FormRectangularLinearSystem(
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
RecoverFEMSolution(const Vector & X,const Vector & b,Vector & x)85 void Operator::RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
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
SetupRAP(const Operator * Pi,const Operator * Po)105 Operator * Operator::SetupRAP(const Operator *Pi, const Operator *Po)
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
FormConstrainedSystemOperator(const Array<int> & ess_tdof_list,ConstrainedOperator * & Aout)134 void Operator::FormConstrainedSystemOperator(
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
FormRectangularConstrainedSystemOperator(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,RectangularConstrainedOperator * & Aout)147 void Operator::FormRectangularConstrainedSystemOperator(
148 const Array<int> &trial_tdof_list, const Array<int> &test_tdof_list,
149 RectangularConstrainedOperator* &Aout)
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
157 RectangularConstrainedOperator *A
158 = new RectangularConstrainedOperator(rap,
159 trial_tdof_list, test_tdof_list,
160 rap != this);
161 Aout = A;
162 }
163
FormSystemOperator(const Array<int> & ess_tdof_list,Operator * & Aout)164 void Operator::FormSystemOperator(const Array<int> &ess_tdof_list,
165 Operator* &Aout)
166 {
167 ConstrainedOperator *A;
168 FormConstrainedSystemOperator(ess_tdof_list, A);
169 Aout = A;
170 }
171
FormRectangularSystemOperator(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,Operator * & Aout)172 void Operator::FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
173 const Array<int> &test_tdof_list,
174 Operator* &Aout)
175 {
176 RectangularConstrainedOperator *A;
177 FormRectangularConstrainedSystemOperator(trial_tdof_list, test_tdof_list, A);
178 Aout = A;
179 }
180
FormDiscreteOperator(Operator * & Aout)181 void Operator::FormDiscreteOperator(Operator* &Aout)
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
PrintMatlab(std::ostream & out,int n,int m) const188 void Operator::PrintMatlab(std::ostream & out, 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 out << 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 out << j+1 << " " << i+1 << " " << y(j) << '\n';
207 }
208 }
209 x(i) = 0.0;
210 }
211 }
212
213
ExplicitMult(const Vector &,Vector &) const214 void TimeDependentOperator::ExplicitMult(const Vector &, Vector &) const
215 {
216 mfem_error("TimeDependentOperator::ExplicitMult() is not overridden!");
217 }
218
ImplicitMult(const Vector &,const Vector &,Vector &) const219 void TimeDependentOperator::ImplicitMult(const Vector &, const Vector &,
220 Vector &) const
221 {
222 mfem_error("TimeDependentOperator::ImplicitMult() is not overridden!");
223 }
224
Mult(const Vector &,Vector &) const225 void TimeDependentOperator::Mult(const Vector &, Vector &) const
226 {
227 mfem_error("TimeDependentOperator::Mult() is not overridden!");
228 }
229
ImplicitSolve(const double,const Vector &,Vector &)230 void TimeDependentOperator::ImplicitSolve(const double, const Vector &,
231 Vector &)
232 {
233 mfem_error("TimeDependentOperator::ImplicitSolve() is not overridden!");
234 }
235
GetImplicitGradient(const Vector &,const Vector &,double) const236 Operator &TimeDependentOperator::GetImplicitGradient(
237 const Vector &, const Vector &, double) const
238 {
239 mfem_error("TimeDependentOperator::GetImplicitGradient() is "
240 "not overridden!");
241 return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
242 }
243
GetExplicitGradient(const Vector &) const244 Operator &TimeDependentOperator::GetExplicitGradient(const Vector &) const
245 {
246 mfem_error("TimeDependentOperator::GetExplicitGradient() is "
247 "not overridden!");
248 return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
249 }
250
SUNImplicitSetup(const Vector &,const Vector &,int,int *,double)251 int TimeDependentOperator::SUNImplicitSetup(const Vector &,
252 const Vector &,
253 int, int *, double)
254 {
255 mfem_error("TimeDependentOperator::SUNImplicitSetup() is not overridden!");
256 return (-1);
257 }
258
SUNImplicitSolve(const Vector &,Vector &,double)259 int TimeDependentOperator::SUNImplicitSolve(const Vector &, Vector &, double)
260 {
261 mfem_error("TimeDependentOperator::SUNImplicitSolve() is not overridden!");
262 return (-1);
263 }
264
SUNMassSetup()265 int TimeDependentOperator::SUNMassSetup()
266 {
267 mfem_error("TimeDependentOperator::SUNMassSetup() is not overridden!");
268 return (-1);
269 }
270
SUNMassSolve(const Vector &,Vector &,double)271 int TimeDependentOperator::SUNMassSolve(const Vector &, Vector &, double)
272 {
273 mfem_error("TimeDependentOperator::SUNMassSolve() is not overridden!");
274 return (-1);
275 }
276
SUNMassMult(const Vector &,Vector &)277 int TimeDependentOperator::SUNMassMult(const Vector &, Vector &)
278 {
279 mfem_error("TimeDependentOperator::SUNMassMult() is not overridden!");
280 return (-1);
281 }
282
283
Mult(const Vector & x,const Vector & dxdt,Vector & y) const284 void SecondOrderTimeDependentOperator::Mult(const Vector &x,
285 const Vector &dxdt,
286 Vector &y) const
287 {
288 mfem_error("SecondOrderTimeDependentOperator::Mult() is not overridden!");
289 }
290
ImplicitSolve(const double dt0,const double dt1,const Vector & x,const Vector & dxdt,Vector & k)291 void SecondOrderTimeDependentOperator::ImplicitSolve(const double dt0,
292 const double dt1,
293 const Vector &x,
294 const Vector &dxdt,
295 Vector &k)
296 {
297 mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
298 }
299
300
ProductOperator(const Operator * A,const Operator * B,bool ownA,bool ownB)301 ProductOperator::ProductOperator(const Operator *A, const Operator *B,
302 bool ownA, bool ownB)
303 : Operator(A->Height(), B->Width()),
304 A(A), B(B), ownA(ownA), ownB(ownB), z(A->Width())
305 {
306 MFEM_VERIFY(A->Width() == B->Height(),
307 "incompatible Operators: A->Width() = " << A->Width()
308 << ", B->Height() = " << B->Height());
309
310 {
311 const Solver* SolverB = dynamic_cast<const Solver*>(B);
312 if (SolverB)
313 {
314 MFEM_VERIFY(!(SolverB->iterative_mode),
315 "Operator B of a ProductOperator should not be in iterative mode");
316 }
317 }
318 }
319
~ProductOperator()320 ProductOperator::~ProductOperator()
321 {
322 if (ownA) { delete A; }
323 if (ownB) { delete B; }
324 }
325
326
RAPOperator(const Operator & Rt_,const Operator & A_,const Operator & P_)327 RAPOperator::RAPOperator(const Operator &Rt_, const Operator &A_,
328 const Operator &P_)
329 : Operator(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_)
330 {
331 MFEM_VERIFY(Rt.Height() == A.Height(),
332 "incompatible Operators: Rt.Height() = " << Rt.Height()
333 << ", A.Height() = " << A.Height());
334 MFEM_VERIFY(A.Width() == P.Height(),
335 "incompatible Operators: A.Width() = " << A.Width()
336 << ", P.Height() = " << P.Height());
337
338 {
339 const Solver* SolverA = dynamic_cast<const Solver*>(&A);
340 if (SolverA)
341 {
342 MFEM_VERIFY(!(SolverA->iterative_mode),
343 "Operator A of an RAPOperator should not be in iterative mode");
344 }
345
346 const Solver* SolverP = dynamic_cast<const Solver*>(&P);
347 if (SolverP)
348 {
349 MFEM_VERIFY(!(SolverP->iterative_mode),
350 "Operator P of an RAPOperator should not be in iterative mode");
351 }
352 }
353
354 mem_class = Rt.GetMemoryClass()*P.GetMemoryClass();
355 MemoryType mem_type = GetMemoryType(A.GetMemoryClass()*mem_class);
356 Px.SetSize(P.Height(), mem_type);
357 APx.SetSize(A.Height(), mem_type);
358 }
359
360
TripleProductOperator(const Operator * A,const Operator * B,const Operator * C,bool ownA,bool ownB,bool ownC)361 TripleProductOperator::TripleProductOperator(
362 const Operator *A, const Operator *B, const Operator *C,
363 bool ownA, bool ownB, bool ownC)
364 : Operator(A->Height(), C->Width())
365 , A(A), B(B), C(C)
366 , ownA(ownA), ownB(ownB), ownC(ownC)
367 {
368 MFEM_VERIFY(A->Width() == B->Height(),
369 "incompatible Operators: A->Width() = " << A->Width()
370 << ", B->Height() = " << B->Height());
371 MFEM_VERIFY(B->Width() == C->Height(),
372 "incompatible Operators: B->Width() = " << B->Width()
373 << ", C->Height() = " << C->Height());
374
375 {
376 const Solver* SolverB = dynamic_cast<const Solver*>(B);
377 if (SolverB)
378 {
379 MFEM_VERIFY(!(SolverB->iterative_mode),
380 "Operator B of a TripleProductOperator should not be in iterative mode");
381 }
382
383 const Solver* SolverC = dynamic_cast<const Solver*>(C);
384 if (SolverC)
385 {
386 MFEM_VERIFY(!(SolverC->iterative_mode),
387 "Operator C of a TripleProductOperator should not be in iterative mode");
388 }
389 }
390
391 mem_class = A->GetMemoryClass()*C->GetMemoryClass();
392 MemoryType mem_type = GetMemoryType(mem_class*B->GetMemoryClass());
393 t1.SetSize(C->Height(), mem_type);
394 t2.SetSize(B->Height(), mem_type);
395 }
396
~TripleProductOperator()397 TripleProductOperator::~TripleProductOperator()
398 {
399 if (ownA) { delete A; }
400 if (ownB) { delete B; }
401 if (ownC) { delete C; }
402 }
403
404
ConstrainedOperator(Operator * A,const Array<int> & list,bool own_A_,DiagonalPolicy diag_policy_)405 ConstrainedOperator::ConstrainedOperator(Operator *A, const Array<int> &list,
406 bool own_A_,
407 DiagonalPolicy diag_policy_)
408 : Operator(A->Height(), A->Width()), A(A), own_A(own_A_),
409 diag_policy(diag_policy_)
410 {
411 // 'mem_class' should work with A->Mult() and MFEM_FORALL():
412 mem_class = A->GetMemoryClass()*Device::GetDeviceMemoryClass();
413 MemoryType mem_type = GetMemoryType(mem_class);
414 list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
415 constraint_list.MakeRef(list);
416 // typically z and w are large vectors, so store them on the device
417 z.SetSize(height, mem_type); z.UseDevice(true);
418 w.SetSize(height, mem_type); w.UseDevice(true);
419 }
420
AssembleDiagonal(Vector & diag) const421 void ConstrainedOperator::AssembleDiagonal(Vector &diag) const
422 {
423 A->AssembleDiagonal(diag);
424
425 if (diag_policy == DIAG_KEEP) { return; }
426
427 const int csz = constraint_list.Size();
428 auto d_diag = diag.ReadWrite();
429 auto idx = constraint_list.Read();
430 switch (diag_policy)
431 {
432 case DIAG_ONE:
433 MFEM_FORALL(i, csz,
434 {
435 const int id = idx[i];
436 d_diag[id] = 1.0;
437 });
438 break;
439 case DIAG_ZERO:
440 MFEM_FORALL(i, csz,
441 {
442 const int id = idx[i];
443 d_diag[id] = 0.0;
444 });
445 break;
446 default:
447 MFEM_ABORT("unknown diagonal policy");
448 break;
449 }
450 }
451
EliminateRHS(const Vector & x,Vector & b) const452 void ConstrainedOperator::EliminateRHS(const Vector &x, Vector &b) const
453 {
454 w = 0.0;
455 const int csz = constraint_list.Size();
456 auto idx = constraint_list.Read();
457 auto d_x = x.Read();
458 // Use read+write access - we are modifying sub-vector of w
459 auto d_w = w.ReadWrite();
460 MFEM_FORALL(i, csz,
461 {
462 const int id = idx[i];
463 d_w[id] = d_x[id];
464 });
465
466 // A.AddMult(w, b, -1.0); // if available to all Operators
467 A->Mult(w, z);
468 b -= z;
469
470 // Use read+write access - we are modifying sub-vector of b
471 auto d_b = b.ReadWrite();
472 MFEM_FORALL(i, csz,
473 {
474 const int id = idx[i];
475 d_b[id] = d_x[id];
476 });
477 }
478
Mult(const Vector & x,Vector & y) const479 void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
480 {
481 const int csz = constraint_list.Size();
482 if (csz == 0)
483 {
484 A->Mult(x, y);
485 return;
486 }
487
488 z = x;
489
490 auto idx = constraint_list.Read();
491 // Use read+write access - we are modifying sub-vector of z
492 auto d_z = z.ReadWrite();
493 MFEM_FORALL(i, csz, d_z[idx[i]] = 0.0;);
494
495 A->Mult(z, y);
496
497 auto d_x = x.Read();
498 // Use read+write access - we are modifying sub-vector of y
499 auto d_y = y.ReadWrite();
500 switch (diag_policy)
501 {
502 case DIAG_ONE:
503 MFEM_FORALL(i, csz,
504 {
505 const int id = idx[i];
506 d_y[id] = d_x[id];
507 });
508 break;
509 case DIAG_ZERO:
510 MFEM_FORALL(i, csz,
511 {
512 const int id = idx[i];
513 d_y[id] = 0.0;
514 });
515 break;
516 case DIAG_KEEP:
517 // Needs action of the operator diagonal on vector
518 mfem_error("ConstrainedOperator::Mult #1");
519 break;
520 default:
521 mfem_error("ConstrainedOperator::Mult #2");
522 break;
523 }
524 }
525
RectangularConstrainedOperator(Operator * A,const Array<int> & trial_list,const Array<int> & test_list,bool own_A_)526 RectangularConstrainedOperator::RectangularConstrainedOperator(
527 Operator *A,
528 const Array<int> &trial_list,
529 const Array<int> &test_list,
530 bool own_A_)
531 : Operator(A->Height(), A->Width()), A(A), own_A(own_A_)
532 {
533 // 'mem_class' should work with A->Mult() and MFEM_FORALL():
534 mem_class = A->GetMemoryClass()*Device::GetMemoryClass();
535 MemoryType mem_type = GetMemoryType(mem_class);
536 trial_list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
537 test_list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
538 trial_constraints.MakeRef(trial_list);
539 test_constraints.MakeRef(test_list);
540 // typically z and w are large vectors, so store them on the device
541 z.SetSize(height, mem_type); z.UseDevice(true);
542 w.SetSize(width, mem_type); w.UseDevice(true);
543 }
544
EliminateRHS(const Vector & x,Vector & b) const545 void RectangularConstrainedOperator::EliminateRHS(const Vector &x,
546 Vector &b) const
547 {
548 w = 0.0;
549 const int trial_csz = trial_constraints.Size();
550 auto trial_idx = trial_constraints.Read();
551 auto d_x = x.Read();
552 // Use read+write access - we are modifying sub-vector of w
553 auto d_w = w.ReadWrite();
554 MFEM_FORALL(i, trial_csz,
555 {
556 const int id = trial_idx[i];
557 d_w[id] = d_x[id];
558 });
559
560 // A.AddMult(w, b, -1.0); // if available to all Operators
561 A->Mult(w, z);
562 b -= z;
563
564 const int test_csz = test_constraints.Size();
565 auto test_idx = test_constraints.Read();
566 auto d_b = b.ReadWrite();
567 MFEM_FORALL(i, test_csz, d_b[test_idx[i]] = 0.0;);
568 }
569
Mult(const Vector & x,Vector & y) const570 void RectangularConstrainedOperator::Mult(const Vector &x, Vector &y) const
571 {
572 const int trial_csz = trial_constraints.Size();
573 const int test_csz = test_constraints.Size();
574 if (trial_csz == 0)
575 {
576 A->Mult(x, y);
577 }
578 else
579 {
580 w = x;
581
582 auto idx = trial_constraints.Read();
583 // Use read+write access - we are modifying sub-vector of w
584 auto d_w = w.ReadWrite();
585 MFEM_FORALL(i, trial_csz, d_w[idx[i]] = 0.0;);
586
587 A->Mult(w, y);
588 }
589
590 if (test_csz != 0)
591 {
592 auto idx = test_constraints.Read();
593 auto d_y = y.ReadWrite();
594 MFEM_FORALL(i, test_csz, d_y[idx[i]] = 0.0;);
595 }
596 }
597
MultTranspose(const Vector & x,Vector & y) const598 void RectangularConstrainedOperator::MultTranspose(const Vector &x,
599 Vector &y) const
600 {
601 const int trial_csz = trial_constraints.Size();
602 const int test_csz = test_constraints.Size();
603 if (test_csz == 0)
604 {
605 A->MultTranspose(x, y);
606 }
607 else
608 {
609 z = x;
610
611 auto idx = test_constraints.Read();
612 // Use read+write access - we are modifying sub-vector of z
613 auto d_z = z.ReadWrite();
614 MFEM_FORALL(i, test_csz, d_z[idx[i]] = 0.0;);
615
616 A->MultTranspose(z, y);
617 }
618
619 if (trial_csz != 0)
620 {
621 auto idx = trial_constraints.Read();
622 auto d_y = y.ReadWrite();
623 MFEM_FORALL(i, trial_csz, d_y[idx[i]] = 0.0;);
624 }
625 }
626
EstimateLargestEigenvalue(Operator & opr,Vector & v0,int numSteps,double tolerance,int seed)627 double PowerMethod::EstimateLargestEigenvalue(Operator& opr, Vector& v0,
628 int numSteps, double tolerance, int seed)
629 {
630 v1.SetSize(v0.Size());
631 if (seed != 0)
632 {
633 v0.Randomize(seed);
634 }
635
636 double eigenvalue = 1.0;
637
638 for (int iter = 0; iter < numSteps; ++iter)
639 {
640 double normV0;
641
642 #ifdef MFEM_USE_MPI
643 if (comm != MPI_COMM_NULL)
644 {
645 normV0 = InnerProduct(comm, v0, v0);
646 }
647 else
648 {
649 normV0 = InnerProduct(v0, v0);
650 }
651 #else
652 normV0 = InnerProduct(v0, v0);
653 #endif
654
655 v0 /= sqrt(normV0);
656 opr.Mult(v0, v1);
657
658 double eigenvalueNew;
659 #ifdef MFEM_USE_MPI
660 if (comm != MPI_COMM_NULL)
661 {
662 eigenvalueNew = InnerProduct(comm, v0, v1);
663 }
664 else
665 {
666 eigenvalueNew = InnerProduct(v0, v1);
667 }
668 #else
669 eigenvalueNew = InnerProduct(v0, v1);
670 #endif
671 double diff = std::abs((eigenvalueNew - eigenvalue) / eigenvalue);
672
673 eigenvalue = eigenvalueNew;
674 std::swap(v0, v1);
675
676 if (diff < tolerance)
677 {
678 break;
679 }
680 }
681
682 return eigenvalue;
683 }
684
685 }
686