complex_smoother.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Include guards
27 #ifndef OOMPH_COMPLEX_SMOOTHER_HEADER
28 #define OOMPH_COMPLEX_SMOOTHER_HEADER
29 
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32 #include <oomph-lib-config.h>
33 #endif
34 
35 // Namespace extension
36 namespace oomph
37 {
38  //====================================================================
39  /// Helmholtz smoother class:
40  /// The smoother class is designed for the Helmholtz equation to be used
41  /// in conjunction with multigrid. The action of the smoother should
42  /// reduce the high frequency errors. These methods are inefficient as
43  /// stand-alone solvers
44  //====================================================================
46  {
47  public:
48  /// Empty constructor
50 
51  /// Virtual empty destructor
52  virtual ~HelmholtzSmoother(){};
53 
54  /// The smoother_solve function performs fixed number of iterations
55  /// on the system A*result=rhs. The number of (smoothing) iterations is
56  /// the same as the max. number of iterations in the underlying
57  /// IterativeLinearSolver class.
59  Vector<DoubleVector>& result) = 0;
60 
61 
62  /// Setup the smoother for the matrix specified by the pointer
63  virtual void complex_smoother_setup(Vector<CRDoubleMatrix*> matrix_pt) = 0;
64 
65  /// Helper function to calculate a complex matrix-vector product.
66  /// Assumes the matrix has been provided as a Vector of length two; the
67  /// first entry containing the real part of the system matrix and the
68  /// second entry containing the imaginary part
70  const Vector<DoubleVector>& x,
72  {
73 #ifdef PARANOID
74  // PARANOID check - Make sure the input matrix has the right size
75  if (matrices_pt.size() != 2)
76  {
77  // Create an output stream
78  std::ostringstream error_message_stream;
79 
80  // Create the error message
81  error_message_stream << "Can only deal with two matrices. You have "
82  << matrices_pt.size() << " matrices." << std::endl;
83 
84  // Throw an error
85  throw OomphLibError(error_message_stream.str(),
86  OOMPH_CURRENT_FUNCTION,
87  OOMPH_EXCEPTION_LOCATION);
88  }
89  // PARANOID check - Make sure the vector x has the right size
90  if (x.size() != 2)
91  {
92  // Create an output stream
93  std::ostringstream error_message_stream;
94 
95  // Create the error message
96  error_message_stream
97  << "Can only deal with two input vectors. You have " << x.size()
98  << " vectors." << std::endl;
99 
100  // Throw an error
101  throw OomphLibError(error_message_stream.str(),
102  OOMPH_CURRENT_FUNCTION,
103  OOMPH_EXCEPTION_LOCATION);
104  }
105  // PARANOID check - Make sure the vector soln has the right size
106  if (soln.size() != 2)
107  {
108  // Create an output stream
109  std::ostringstream error_message_stream;
110 
111  // Create the error message
112  error_message_stream
113  << "Can only deal with two output vectors. You have " << soln.size()
114  << " output vectors." << std::endl;
115 
116  // Throw an error
117  throw OomphLibError(error_message_stream.str(),
118  OOMPH_CURRENT_FUNCTION,
119  OOMPH_EXCEPTION_LOCATION);
120  }
121 #endif
122 
123  //-----------------------------------------------------------------------
124  // Suppose we have a complex matrix, A, and two complex vectors, x and
125  // soln. We wish to compute the product A*x=soln (note, * does not mean
126  // we are using complex conjugates here, it is simply used to indicate
127  // a multiplication). To do this we must make use of the fact that we
128  // possess the real and imaginary separately. As a result, it is computed
129  // using:
130  // soln = A*x,
131  // = (A_r + i*A_c)*(x_r + i*x_c),
132  // = [A_r*x_r - A_c*x_c] + i*[A_r*x_c + A_c*x_r],
133  // ==> real(soln) = A_r*x_r - A_c*x_c,
134  // & imag(soln) = A_r*x_c + A_c*x_r,
135  // where the subscripts _r and _c are used to identify the real and
136  // imaginary part, respectively.
137  //-----------------------------------------------------------------------
138 
139  // Store the value of A_r*x_r in the real part of soln
140  matrices_pt[0]->multiply(x[0], soln[0]);
141 
142  // Store the value of A_r*x_c in the imaginary part of soln
143  matrices_pt[0]->multiply(x[1], soln[1]);
144 
145  // Create a temporary vector
146  DoubleVector temp(this->distribution_pt(), 0.0);
147 
148  // Calculate the value of A_c*x_c
149  matrices_pt[1]->multiply(x[1], temp);
150 
151  // Subtract the value of temp from soln[0] to get the real part of soln
152  soln[0] -= temp;
153 
154  // Calculate the value of A_c*x_r
155  matrices_pt[1]->multiply(x[0], temp);
156 
157  // Add the value of temp to soln[1] to get the imaginary part of soln
158  soln[1] += temp;
159  } // End of complex_matrix_multiplication
160 
161  /// Self-test to check that all the dimensions of the inputs to
162  /// solve helper are consistent and everything that needs to be built, is.
163  template<typename MATRIX>
165  CRDoubleMatrix* const& real_matrix_pt,
166  CRDoubleMatrix* const& imag_matrix_pt,
167  const Vector<DoubleVector>& rhs,
168  Vector<DoubleVector>& solution,
169  const double& n_dof);
170 
171  protected:
172  /// When a derived class object is being used as a smoother in
173  /// the MG algorithm the residual norm
174  /// does not need to be calculated. This boolean is used as a flag to
175  /// indicate this in solve_helper(...)
177  };
178 
179  //==================================================================
180  /// Self-test to be called inside solve_helper to ensure
181  /// that all inputs are consistent and everything that needs to
182  /// be built, is.
183  //==================================================================
184  template<typename MATRIX>
186  CRDoubleMatrix* const& real_matrix_pt,
187  CRDoubleMatrix* const& imag_matrix_pt,
188  const Vector<DoubleVector>& rhs,
189  Vector<DoubleVector>& solution,
190  const double& n_dof)
191  {
192  // Number of dof types should be 2 (real & imaginary)
193  unsigned n_dof_types = 2;
194 
195  // Create a vector to hold the matrices
196  Vector<CRDoubleMatrix*> matrix_storage_pt(2, 0);
197 
198  // Assign the first entry in matrix_storage_pt
199  matrix_storage_pt[0] = real_matrix_pt;
200 
201  // Assign the second entry in matrix_storage_pt
202  matrix_storage_pt[1] = imag_matrix_pt;
203 
204  // Loop over the real and imaginary parts
205  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
206  {
207  // Check if the matrix is distributable. If it is then it should
208  // not be distributed
209  if (dynamic_cast<DistributableLinearAlgebraObject*>(
210  matrix_storage_pt[dof_type]) != 0)
211  {
212  if (dynamic_cast<DistributableLinearAlgebraObject*>(
213  matrix_storage_pt[dof_type])
214  ->distributed())
215  {
216  std::ostringstream error_message_stream;
217  error_message_stream << "The matrix must not be distributed.";
218  throw OomphLibError(error_message_stream.str(),
219  OOMPH_CURRENT_FUNCTION,
220  OOMPH_EXCEPTION_LOCATION);
221  }
222  }
223  // Check that this rhs distribution is setup
224  if (!(rhs[dof_type].built()))
225  {
226  std::ostringstream error_message_stream;
227  error_message_stream << "The rhs vector distribution must be setup.";
228  throw OomphLibError(error_message_stream.str(),
229  OOMPH_CURRENT_FUNCTION,
230  OOMPH_EXCEPTION_LOCATION);
231  }
232  // Check that the rhs has the right number of global rows
233  if (rhs[dof_type].nrow() != n_dof)
234  {
235  std::ostringstream error_message_stream;
236  error_message_stream << "RHS does not have the same dimension as the "
237  << "linear system";
238  throw OomphLibError(error_message_stream.str(),
239  OOMPH_CURRENT_FUNCTION,
240  OOMPH_EXCEPTION_LOCATION);
241  }
242  // Check that the rhs is not distributed
243  if (rhs[dof_type].distribution_pt()->distributed())
244  {
245  std::ostringstream error_message_stream;
246  error_message_stream << "The rhs vector must not be distributed.";
247  throw OomphLibError(error_message_stream.str(),
248  OOMPH_CURRENT_FUNCTION,
249  OOMPH_EXCEPTION_LOCATION);
250  }
251  // Check that if the result is setup it matches the distribution
252  // of the rhs
253  if (solution[dof_type].built())
254  {
255  if (!(*rhs[dof_type].distribution_pt() ==
256  *solution[dof_type].distribution_pt()))
257  {
258  std::ostringstream error_message_stream;
259  error_message_stream << "If the result distribution is setup then it "
260  << "must be the same as the rhs distribution";
261  throw OomphLibError(error_message_stream.str(),
262  OOMPH_CURRENT_FUNCTION,
263  OOMPH_EXCEPTION_LOCATION);
264  }
265  } // if ((solution[0].built())||(solution[1].built()))
266  } // for (unsigned dof_type=0;dof_type<n_dof_types;dof_type)
267  } // End of check_validity_of_solve_helper_inputs
268 
269 
270  /// ////////////////////////////////////////////////////////////////////
271  /// ////////////////////////////////////////////////////////////////////
272  /// ////////////////////////////////////////////////////////////////////
273 
274 
275  //=========================================================================
276  /// Damped Jacobi "solver" templated by matrix type. The "solver"
277  /// exists in many different incarnations: It's an IterativeLinearSolver,
278  /// and a Smoother, all of which use the same basic iteration.
279  //=========================================================================
280  template<typename MATRIX>
282  {
283  public:
284  /// Constructor (empty)
285  ComplexDampedJacobi(const double& omega = 0.5)
286  : Matrix_can_be_deleted(true),
287  Matrix_real_pt(0),
288  Matrix_imag_pt(0),
289  Omega(omega){};
290 
291  /// Empty destructor
293  {
294  // Run clean_up_memory
295  clean_up_memory();
296  } // End of ~ComplexDampedJacobi
297 
298  /// Cleanup data that's stored for resolve (if any has been stored)
300  {
301  // If the real matrix pointer isn't null AND we're allowed to delete
302  // the matrix which is only when we create the matrix ourselves
303  if ((Matrix_real_pt != 0) && (Matrix_can_be_deleted))
304  {
305  // Delete the matrix
306  delete Matrix_real_pt;
307 
308  // Assign the associated pointer the value NULL
309  Matrix_real_pt = 0;
310  }
311 
312  // If the real matrix pointer isn't null AND we're allowed to delete
313  // the matrix which is only when we create the matrix ourselves
314  if ((Matrix_imag_pt != 0) && (Matrix_can_be_deleted))
315  {
316  // Delete the matrix
317  delete Matrix_imag_pt;
318 
319  // Assign the associated pointer the value NULL
320  Matrix_imag_pt = 0;
321  }
322  } // End of clean_up_memory
323 
324  /// Broken copy constructor
326 
327  /// Broken assignment operator
328  void operator=(const ComplexDampedJacobi&) = delete;
329 
330  /// Function to calculate the value of Omega by passing in the
331  /// value of k and h [see Elman et al. "A multigrid method enhanced by
332  /// Krylov subspace iteration for discrete Helmholtz equations", p.1303]
333  void calculate_omega(const double& k, const double& h)
334  {
335  // Create storage for the parameter kh
336  double kh = k * h;
337 
338  // Store the value of pi
339  double pi = MathematicalConstants::Pi;
340 
341  // Calculate the value of Omega
342  double omega = (12.0 - 4.0 * pow(kh, 2.0)) / (18.0 - 3.0 * pow(kh, 2.0));
343 
344  // Set the tolerance for how close omega can be to 0
345  double tolerance = 1.0e-03;
346 
347  // Only store the value of omega if it lies in the range (0,1]. If it
348  // isn't it will produce a divergent scheme. Note, to avoid letting
349  // omega become too small we make sure it is greater than some tolerance
350  if ((omega > tolerance) && !(omega > 1))
351  {
352  // Since omega lies in the specified range, store it
353  Omega = omega;
354  }
355  // On the coarsest grids: if the wavenumber is greater than pi and
356  // kh>2cos(pi*h/2) then we can choose omega from the range (0,omega_2)
357  // where omega_2 is defined in [Elman et al."A multigrid method
358  // enhanced by Krylov subspace iteration for discrete Helmholtz
359  // equations", p.1295]
360  else if ((k > pi) && (kh > 2 * cos(pi * h / 2)))
361  {
362  // Calculate the numerator of omega_2
363  double omega_2 = (2.0 - pow(kh, 2.0));
364 
365  // Divide by the denominator of the fraction
366  omega_2 /= (2.0 * pow(sin(pi * h / 2), 2.0) - 0.5 * pow(kh, 2.0));
367 
368  // If 2/3 lies in the range (0,omega_2), use it
369  if (omega_2 > (2.0 / 3.0))
370  {
371  // Assign Omega the damping factor used for the Poisson equation
372  Omega = 2.0 / 3.0;
373  }
374  // If omega_2 is less than 2/3 use the midpoint of (tolerance,omega_2)
375  else
376  {
377  // Set the value of Omega
378  Omega = 0.5 * (tolerance + omega_2);
379  }
380  }
381  // Since the value of kh must be fairly large, make the value of
382  // omega small to compensate
383  else
384  {
385  // Since kh doesn't lie in the chosen range, set it to some small value
386  Omega = 0.2;
387  }
388  } // End of calculate_omega
389 
390  /// Get access to the value of Omega (lvalue)
391  double& omega()
392  {
393  // Return the value of Omega
394  return Omega;
395  } // End of omega
396 
397  /// Setup: Pass pointer to the matrix and store in cast form
399  {
400  // Assume the matrices have been passed in from the outside so we must
401  // not delete it. This is needed to avoid pre- and post-smoothers
402  // deleting the same matrix in the MG solver. If this was originally
403  // set to TRUE then this will be sorted out in the other functions
404  // from which this was called
405  Matrix_can_be_deleted = false;
406 
407  // Assign the real matrix pointers
408  Matrix_real_pt = helmholtz_matrix_pt[0];
409 
410  // Assign the imaginary matrix pointers
411  Matrix_imag_pt = helmholtz_matrix_pt[1];
412 
413 #ifdef PARANOID
414  // PARANOID check that if the matrix is distributable. If it is not then
415  // it should not be distributed
416  if (Matrix_real_pt->nrow() != Matrix_imag_pt->nrow())
417  {
418  std::ostringstream error_message_stream;
419  error_message_stream << "Incompatible real and complex matrix sizes.";
420  throw OomphLibError(error_message_stream.str(),
421  OOMPH_CURRENT_FUNCTION,
422  OOMPH_EXCEPTION_LOCATION);
423  }
424 #endif
425 
426  //--------------------------------------------------------------------
427  // We need the matrix diagonals to calculate inv(D) (where D is the
428  // diagonal of A) as it remains the same for each use of the iterative
429  // scheme. To avoid unnecessary computations we store it now so it can
430  // simply be called in each iteration.
431  //--------------------------------------------------------------------
432 
433  // Grab the diagonal entries of the real part of the system matrix
435 
436  // Grab the diagonal entries of the imaginary part of the system matrix
438 
439  // Find the number of entries in the vectors containing the diagonal
440  // entries of both the real and the imaginary matrix
441  unsigned n_dof = Matrix_diagonal_real.size();
442 
443  // Make a dummy vector to store the entries of Matrix_diagonal_inv_sq
444  Matrix_diagonal_inv_sq.resize(n_dof, 0.0);
445 
446  // Create a helper variable to store A_r(j,j), for some j
447  double dummy_r;
448 
449  // Create a helper variable to store A_c(j,j), for some j
450  double dummy_c;
451 
452  // Loop over the entries of Matrix_diagonal_inv_sq and set the i-th
453  // entry to 1/(A_r(i,i)^2 + A_c(i,i)^2)
454  for (unsigned i = 0; i < n_dof; i++)
455  {
456  // Store the value of A_r(i,i)
457  dummy_r = Matrix_diagonal_real[i];
458 
459  // Store the value of A_c(i,i)
460  dummy_c = Matrix_diagonal_imag[i];
461 
462  // Store the value of 1/(A_r(i,i)^2 + A_c(i,i)^2)
464  1.0 / (dummy_r * dummy_r + dummy_c * dummy_c);
465  }
466  } // End of complex_smoother_setup
467 
468  /// The smoother_solve function performs fixed number of iterations
469  /// on the system A*result=rhs. The number of (smoothing) iterations is
470  /// the same as the max. number of iterations in the underlying
471  /// IterativeLinearSolver class.
473  Vector<DoubleVector>& solution)
474  {
475  // If you use a smoother but you don't want to calculate the residual
476  Use_as_smoother = true;
477 
478  // Call the helper function
479  complex_solve_helper(rhs, solution);
480  } // End of complex_smoother_solve
481 
482  /// Use damped Jacobi iteration as an IterativeLinearSolver:
483  /// This obtains the Jacobian matrix J and the residual vector r
484  /// (needed for the Newton method) from the problem's get_jacobian
485  /// function and returns the result of Jx=r.
486  void solve(Problem* const& problem_pt, DoubleVector& result)
487  {
488  BrokenCopy::broken_assign("ComplexDampedJacobi");
489  }
490 
491  /// Number of iterations taken
492  unsigned iterations() const
493  {
494  return Iterations;
495  }
496 
497  private:
498  /// This is where the actual work is done
500  Vector<DoubleVector>& solution);
501 
502  /// Boolean flag to indicate if the matrices pointed to by
503  /// Matrix_real_pt and Matrix_imag_pt can be deleted.
505 
506  /// Pointer to the real part of the system matrix
508 
509  /// Pointer to the real part of the system matrix
511 
512  /// Vector containing the diagonal entries of A_r (real(A))
514 
515  /// Vector containing the diagonal entries of A_c (imag(A))
517 
518  /// Vector whose j'th entry is given by 1/(A_r(j,j)^2 + A_c(j,j)^2)
520 
521  /// Number of iterations taken
522  unsigned Iterations;
523 
524  /// Damping factor
525  double Omega;
526  };
527 
528  //======================================================================
529  /// This is where the actual work is done.
530  //======================================================================
531  template<typename MATRIX>
533  const Vector<DoubleVector>& rhs, Vector<DoubleVector>& solution)
534  {
535  // Get number of dofs
536  unsigned n_dof = Matrix_real_pt->nrow();
537 
538 #ifdef PARANOID
539  // Upcast the matrix to the appropriate type
540  CRDoubleMatrix* tmp_rmatrix_pt =
541  dynamic_cast<CRDoubleMatrix*>(Matrix_real_pt);
542 
543  // Upcast the matrix to the appropriate type
544  CRDoubleMatrix* tmp_imatrix_pt =
545  dynamic_cast<CRDoubleMatrix*>(Matrix_imag_pt);
546 
547  // PARANOID Run the self-tests to check the inputs are correct
548  this->check_validity_of_solve_helper_inputs<MATRIX>(
549  tmp_rmatrix_pt, tmp_imatrix_pt, rhs, solution, n_dof);
550 
551  // We don't need the real matrix pointer any more
552  tmp_rmatrix_pt = 0;
553 
554  // We don't need the imaginary matrix pointer any more
555  tmp_imatrix_pt = 0;
556 #endif
557 
558  // Setup the solution if it is not
559  if ((!solution[0].distribution_pt()->built()) ||
560  (!solution[1].distribution_pt()->built()))
561  {
562  solution[0].build(rhs[0].distribution_pt(), 0.0);
563  solution[1].build(rhs[1].distribution_pt(), 0.0);
564  }
565 
566  // Initialise timer
567  double t_start = TimingHelpers::timer();
568 
569  // Copy the real and imaginary part of the solution vector
570  DoubleVector x_real(solution[0]);
571  DoubleVector x_imag(solution[1]);
572 
573  // Copy the real and imaginary part of the RHS vector
574  DoubleVector rhs_real(rhs[0]);
575  DoubleVector rhs_imag(rhs[1]);
576 
577  // Create storage for the real and imaginary part of the constant term
578  DoubleVector constant_term_real(this->distribution_pt(), 0.0);
579  DoubleVector constant_term_imag(this->distribution_pt(), 0.0);
580 
581  // Create storage for the real and imaginary part of the residual vector.
582  // These aren't used/built if we're inside the multigrid solver
583  DoubleVector residual_real;
584  DoubleVector residual_imag;
585 
586  // Create storage for the norm of the real and imaginary parts of the
587  // residual vector. These aren't used if we're inside the multigrid
588  // solver
589  double res_real_norm = 0.0;
590  double res_imag_norm = 0.0;
591 
592  // Variable to hold the current residual norm. This isn't used if
593  // we're inside the multigrid solver
594  double norm_res = 0.0;
595 
596  // Variables to hold the initial residual norm. This isn't used if
597  // we're inside the multigrid solver
598  double norm_f = 0.0;
599 
600  // Initialise the value of Iterations
601  Iterations = 0;
602 
603  //------------------------------------------------------------------------
604  // Initial guess isn't necessarily zero (restricted solution from finer
605  // grids) therefore both x and the residual need to be assigned values
606  // from inputs. The constant term at the end of the stationary iteration
607  // is also calculated here since it does not change at all throughout the
608  // iterative process:
609  //------------------------------------------------------------------------
610 
611  // Loop over all the entries of each vector
612  for (unsigned i = 0; i < n_dof; i++)
613  {
614  // Calculate the numerator of the i'th entry in the real component of
615  // the constant term
616  constant_term_real[i] = (Matrix_diagonal_real[i] * rhs_real[i] +
617  Matrix_diagonal_imag[i] * rhs_imag[i]);
618 
619  // Divide by the denominator
620  constant_term_real[i] *= Matrix_diagonal_inv_sq[i];
621 
622  // Scale by the damping factor
623  constant_term_real[i] *= Omega;
624 
625  // Calculate the numerator of the i'th entry in the imaginary component of
626  // the constant term
627  constant_term_imag[i] = (Matrix_diagonal_real[i] * rhs_imag[i] -
628  Matrix_diagonal_imag[i] * rhs_real[i]);
629 
630  // Divide by the denominator
631  constant_term_imag[i] *= Matrix_diagonal_inv_sq[i];
632 
633  // Scale by the damping factor
634  constant_term_imag[i] *= Omega;
635  }
636 
637  // Create 4 temporary vectors to store the various matrix-vector products
638  // required. The appropriate combination has been appended to the name. For
639  // instance, the product A_r*x_c corresponds to temp_vec_rc (real*imag)
640  DoubleVector temp_vec_rr(this->distribution_pt(), 0.0);
641  DoubleVector temp_vec_cc(this->distribution_pt(), 0.0);
642  DoubleVector temp_vec_rc(this->distribution_pt(), 0.0);
643  DoubleVector temp_vec_cr(this->distribution_pt(), 0.0);
644 
645  // Calculate the different combinations of A*x (or A*e depending on the
646  // level in the heirarchy) in the complex case (i.e. A_r*x_r, A_c*x_c,
647  // A_r*x_c and A_c*x_r)
648  Matrix_real_pt->multiply(x_real, temp_vec_rr);
649  Matrix_imag_pt->multiply(x_imag, temp_vec_cc);
650  Matrix_real_pt->multiply(x_imag, temp_vec_rc);
651  Matrix_imag_pt->multiply(x_real, temp_vec_cr);
652 
653  //---------------------------------------------------------------------------
654  // Calculating the residual r=b-Ax in the complex case requires more care
655  // than the real case. The correct calculation can be derived rather easily.
656  // Consider the splitting of A, x and b into their complex components:
657  // r = b - A*x
658  // = (b_r + i*b_c) - (A_r + i*A_c)*(x_r + i*x_c)
659  // = [b_r - A_r*x_r + A_c*x_c] + i*[b_c - A_r*x_c - A_c*x_r]
660  // ==> real(r) = b_r - A_r*x_r + A_c*x_c
661  // & imag(r) = b_c - A_r*x_c - A_c*x_r.
662  //---------------------------------------------------------------------------
663 
664  // Calculate the residual only if we're not inside the multigrid solver
665  if (!Use_as_smoother)
666  {
667  // Create storage for the real and imaginary part of the residual vector
668  residual_real.build(this->distribution_pt(), 0.0);
669  residual_imag.build(this->distribution_pt(), 0.0);
670 
671  // Real part of the residual:
672  residual_real = rhs_real;
673  residual_real -= temp_vec_rr;
674  residual_real += temp_vec_cc;
675 
676  // Imaginary part of the residual:
677  residual_imag = rhs_imag;
678  residual_imag -= temp_vec_rc;
679  residual_imag -= temp_vec_cr;
680 
681  // Calculate the 2-norm (noting that the 2-norm of a complex vector a of
682  // length n is simply the square root of the sum of the squares of each
683  // real and imaginary component). That is:
684  /// \f[ \|a\|_2^2=\sum_{i=0}^{i=n-1}\Re(a[i])^2+\Im(a[i])^2. \f]
685  // can be written as:
686  /// \f[ \|a\|_2^2=\|\Re(a)\|_2^2+\|\Im(a)\|_2^2. \f]
687  double res_real_norm = residual_real.norm();
688  double res_imag_norm = residual_imag.norm();
689  double norm_res =
690  sqrt(res_real_norm * res_real_norm + res_imag_norm * res_imag_norm);
691 
692  // If required will document convergence history to screen
693  // or file (if stream is open)
694  if (Doc_convergence_history)
695  {
696  if (!Output_file_stream.is_open())
697  {
698  oomph_info << Iterations << " " << norm_res << std::endl;
699  }
700  else
701  {
702  Output_file_stream << Iterations << " " << norm_res << std::endl;
703  }
704  } // if (Doc_convergence_history)
705  } // if (!Use_as_smoother)
706 
707  // Two temporary vectors to store the value of A_r*x_r - A_c*x_c and
708  // A_c*x_r + A_r*x_c in each iteration
709  DoubleVector temp_vec_real(this->distribution_pt(), 0.0);
710  DoubleVector temp_vec_imag(this->distribution_pt(), 0.0);
711 
712  // Calculate A_r*x_r - A_c*x_c
713  temp_vec_real = temp_vec_rr;
714  temp_vec_real -= temp_vec_cc;
715 
716  // Calculate A_c*x_r + A_r*x_c
717  temp_vec_imag = temp_vec_cr;
718  temp_vec_imag += temp_vec_rc;
719 
720  // Outermost loop: Run up to Max_iter times
721  for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
722  {
723  // Loop over each degree of freedom and update
724  // the current approximation
725  for (unsigned i = 0; i < n_dof; i++)
726  {
727  // Make a couple of dummy variables to help with calculations
728  double dummy_r = 0.0;
729  double dummy_c = 0.0;
730 
731  // Calculate one part of the correction term (real)
732  dummy_r = (Matrix_diagonal_real[i] * temp_vec_real[i] +
733  Matrix_diagonal_imag[i] * temp_vec_imag[i]);
734 
735  // Calculate one part of the correction term (imaginary)
736  dummy_c = (Matrix_diagonal_real[i] * temp_vec_imag[i] -
737  Matrix_diagonal_imag[i] * temp_vec_real[i]);
738 
739  // Scale by Omega/([A(i,i)_r]^2+[A(i,i)_c]^2)
740  dummy_r *= Omega * Matrix_diagonal_inv_sq[i];
741  dummy_c *= Omega * Matrix_diagonal_inv_sq[i];
742 
743  // Update the value of x_real
744  x_real[i] -= dummy_r;
745  x_imag[i] -= dummy_c;
746  }
747 
748  // Update the value of x (or e depending on the level in the heirarchy)
749  x_real += constant_term_real;
750  x_imag += constant_term_imag;
751 
752  // Calculate the different combinations of A*x (or A*e depending on the
753  // level in the heirarchy) in the complex case (i.e. A_r*x_r, A_c*x_c,
754  // A_r*x_c and A_c*x_r)
755  Matrix_real_pt->multiply(x_real, temp_vec_rr);
756  Matrix_imag_pt->multiply(x_imag, temp_vec_cc);
757  Matrix_real_pt->multiply(x_imag, temp_vec_rc);
758  Matrix_imag_pt->multiply(x_real, temp_vec_cr);
759 
760  // Calculate A_r*x_r - A_c*x_c
761  temp_vec_real = temp_vec_rr;
762  temp_vec_real -= temp_vec_cc;
763 
764  // Calculate A_c*x_r + A_r*x_c
765  temp_vec_imag = temp_vec_cr;
766  temp_vec_imag += temp_vec_rc;
767 
768  // Calculate the residual only if we're not inside the multigrid solver
769  if (!Use_as_smoother)
770  {
771  // Calculate the residual
772  residual_real = rhs_real;
773  residual_real -= temp_vec_rr;
774  residual_real += temp_vec_cc;
775 
776  // Calculate the imaginary part of the residual vector
777  residual_imag = rhs_imag;
778  residual_imag -= temp_vec_rc;
779  residual_imag -= temp_vec_cr;
780 
781  // Calculate the 2-norm using the method mentioned previously
782  res_real_norm = residual_real.norm();
783  res_imag_norm = residual_imag.norm();
784  norm_res =
785  sqrt(res_real_norm * res_real_norm + res_imag_norm * res_imag_norm) /
786  norm_f;
787 
788  // If required, this will document convergence history to
789  // screen or file (if the stream is open)
790  if (Doc_convergence_history)
791  {
792  if (!Output_file_stream.is_open())
793  {
794  oomph_info << Iterations << " " << norm_res << std::endl;
795  }
796  else
797  {
798  Output_file_stream << Iterations << " " << norm_res << std::endl;
799  }
800  } // if (Doc_convergence_history)
801 
802  // Check the tolerance only if the residual norm is being computed
803  if (norm_res < Tolerance)
804  {
805  // Break out of the for-loop
806  break;
807  }
808  } // if (!Use_as_smoother)
809  } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
810 
811  // Calculate the residual only if we're not inside the multigrid solver
812  if (!Use_as_smoother)
813  {
814  // If time documentation is enabled
815  if (Doc_time)
816  {
817  oomph_info << "\n(Complex) damped Jacobi converged. Residual norm: "
818  << norm_res
819  << "\nNumber of iterations to convergence: " << Iterations
820  << "\n"
821  << std::endl;
822  }
823  } // if (!Use_as_smoother)
824 
825  // Copy the solution into the solution vector
826  solution[0] = x_real;
827  solution[1] = x_imag;
828 
829  // Doc time for solver
830  double t_end = TimingHelpers::timer();
831  Solution_time = t_end - t_start;
832  if (Doc_time)
833  {
834  oomph_info << "Time for solve with (complex) damped Jacobi [sec]: "
835  << Solution_time << std::endl;
836  }
837 
838  // If the solver failed to converge and the user asked for an error if
839  // this happened
840  if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
841  {
842  std::string error_message =
843  "Solver failed to converge and you requested ";
844  error_message += "an error on convergence failures.";
845  throw OomphLibError(
846  error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
847  }
848  } // End of complex_solve_helper function
849 
850  //======================================================================
851  /// The GMRES method rewritten for complex matrices
852  //======================================================================
853  template<typename MATRIX>
855  {
856  public:
857  /// Constructor
859  : Iterations(0),
861  Resolving(false),
863  {
864  } // End of ComplexGMRES constructor
865 
866  /// Empty destructor
868  {
869  // Run clean_up_memory
870  clean_up_memory();
871  } // End of ~ComplexGMRES
872 
873  /// Broken copy constructor
874  ComplexGMRES(const ComplexGMRES&) = delete;
875 
876  /// Broken assignment operator
877  void operator=(const ComplexGMRES&) = delete;
878 
879  /// Overload disable resolve so that it cleans up memory too
881  {
882  // Disable resolve using the LinearSolver function
884 
885  // Clean up anything kept in memory
886  clean_up_memory();
887  } // End of disable resolve
888 
889  /// Solver: Takes pointer to problem and returns the results vector
890  /// which contains the solution of the linear system defined by
891  /// the problem's fully assembled Jacobian and residual vector.
892  void solve(Problem* const& problem_pt, DoubleVector& result)
893  {
894  // Write the error message into a string
895  std::string error_message = "Solve function for class\n\n";
896  error_message += "ComplexGMRES\n\n";
897  error_message += "is deliberately broken to avoid the accidental \n";
898  error_message += "use of the inappropriate C++ default. If you \n";
899  error_message += "really need one for this class, write it yourself...\n";
900 
901  // Throw an error
902  throw OomphLibError(
903  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
904  } // End of solve
905 
906  /// Linear-algebra-type solver: Takes pointer to a matrix
907  /// and rhs vector and returns the solution of the linear system
908  /// Call the broken base-class version. If you want this, please
909  /// implement it
910  void solve(DoubleMatrixBase* const& matrix_pt,
911  const Vector<double>& rhs,
912  Vector<double>& result)
913  {
914  LinearSolver::solve(matrix_pt, rhs, result);
915  } // End of solve
916 
917  /// Number of iterations taken
918  unsigned iterations() const
919  {
920  // Return the number of iterations used
921  return Iterations;
922  } // End of iterations
923 
924  /// Setup: Pass pointer to the matrix and store in cast form
926  {
927  // Assume the matrices have been passed in from the outside so we must
928  // not delete it. This is needed to avoid pre- and post-smoothers
929  // deleting the same matrix in the MG solver. If this was originally
930  // set to TRUE then this will be sorted out in the other functions
931  // from which this was called
932  Matrix_can_be_deleted = false;
933 
934 #ifdef PARANOID
935  // PARANOID check - Make sure the input has the right number of matrices
936  if (helmholtz_matrix_pt.size() != 2)
937  {
938  std::ostringstream error_message_stream;
939  error_message_stream << "Can only deal with two matrices. You have "
940  << helmholtz_matrix_pt.size() << " matrices."
941  << std::endl;
942 
943  // Throw an error
944  throw OomphLibError(error_message_stream.str(),
945  OOMPH_CURRENT_FUNCTION,
946  OOMPH_EXCEPTION_LOCATION);
947  }
948 #endif
949 
950  // Resize the storage for the system matrices
951  Matrices_storage_pt.resize(2, 0);
952 
953  // Assign the real matrix pointer
954  Matrices_storage_pt[0] = helmholtz_matrix_pt[0];
955 
956  // Assign the imaginary matrix pointers
957  Matrices_storage_pt[1] = helmholtz_matrix_pt[1];
958 
959 #ifdef PARANOID
960  // PARANOID check - Make sure that the constituent matrices have the same
961  // number of rows
963  {
964  std::ostringstream error_message_stream;
965  error_message_stream << "Incompatible real and imag. matrix sizes.";
966  throw OomphLibError(error_message_stream.str(),
967  OOMPH_CURRENT_FUNCTION,
968  OOMPH_EXCEPTION_LOCATION);
969  }
970  // PARANOID check - Make sure that the constituent matrices have the same
971  // number of columns
972  if (Matrices_storage_pt[0]->ncol() != Matrices_storage_pt[1]->ncol())
973  {
974  std::ostringstream error_message_stream;
975  error_message_stream << "Incompatible real and imag. matrix sizes.";
976  throw OomphLibError(error_message_stream.str(),
977  OOMPH_CURRENT_FUNCTION,
978  OOMPH_EXCEPTION_LOCATION);
979  }
980 #endif
981  } // End of complex_smoother_setup
982 
983  /// The smoother_solve function performs fixed number of iterations
984  /// on the system A*result=rhs. The number of (smoothing) iterations is
985  /// the same as the max. number of iterations in the underlying
986  /// IterativeLinearSolver class.
988  Vector<DoubleVector>& solution)
989  {
990  // If you use a smoother but you don't want to calculate the residual
991  Use_as_smoother = true;
992 
993  // Use the helper function where the work is actually done
994  complex_solve_helper(rhs, solution);
995  } // End of complex_smoother_solve
996 
997  private:
998  /// Cleanup data that's stored for resolve (if any has been stored)
1000  {
1001  // If the real matrix pointer isn't null AND we're allowed to delete
1002  // the matrix which is only when we create the matrix ourselves
1003  if ((Matrices_storage_pt[0] != 0) && (Matrix_can_be_deleted))
1004  {
1005  // Delete the matrix
1006  delete Matrices_storage_pt[0];
1007 
1008  // Assign the associated pointer the value NULL
1009  Matrices_storage_pt[0] = 0;
1010  }
1011 
1012  // If the real matrix pointer isn't null AND we're allowed to delete
1013  // the matrix which is only when we create the matrix ourselves
1014  if ((Matrices_storage_pt[1] != 0) && (Matrix_can_be_deleted))
1015  {
1016  // Delete the matrix
1017  delete Matrices_storage_pt[1];
1018 
1019  // Assign the associated pointer the value NULL
1020  Matrices_storage_pt[1] = 0;
1021  }
1022  } // End of clean_up_memory
1023 
1024  /// This is where the actual work is done
1026  Vector<DoubleVector>& solution);
1027 
1028  /// Helper function to update the result vector
1029  void update(const unsigned& k,
1030  const Vector<Vector<std::complex<double>>>& hessenberg,
1031  const Vector<std::complex<double>>& s,
1032  const Vector<Vector<DoubleVector>>& v,
1034  {
1035  // Make a local copy of s
1037 
1038  //-----------------------------------------------------------------
1039  // The Hessenberg matrix should be an upper triangular matrix at
1040  // this point (although from its storage it would appear to be a
1041  // lower triangular matrix since the indexing has been reversed)
1042  // so finding the minimiser of J(y)=min||s-R_m*y|| where R_m is
1043  // the matrix R in the QR factorisation of the Hessenberg matrix.
1044  // Therefore, to obtain y we simply need to use a backwards
1045  // substitution. Note: The implementation here may appear to be
1046  // somewhat confusing as the indexing in the Hessenberg matrix is
1047  // reversed. This implementation of a backwards substitution does
1048  // not run along the columns of the triangular matrix but rather
1049  // up the rows.
1050  //-----------------------------------------------------------------
1051  // The outer loop is a loop over the columns of the Hessenberg matrix
1052  // since the indexing is reversed
1053  for (int i = int(k); i >= 0; i--)
1054  {
1055  // Divide the i-th entry of y by the i-th diagonal entry of H
1056  y[i] /= hessenberg[i][i];
1057 
1058  // The inner loop is a loop over the rows of the Hessenberg matrix
1059  for (int j = i - 1; j >= 0; j--)
1060  {
1061  // Update the j-th entry of y
1062  y[j] -= hessenberg[i][j] * y[i];
1063  }
1064  } // for (int i=int(k);i>=0;i--)
1065 
1066  // Calculate the number of entries in x (simply use the real part as
1067  // both the real and imaginary part should have the same length)
1068  unsigned n_row = x[0].nrow();
1069 
1070  // We assume here that the vector x (which is passed in) is actually x_0
1071  // so we simply need to update its entries to calculate the solution, x
1072  // which is given by x=x_0+Vy.
1073  for (unsigned j = 0; j <= k; j++)
1074  {
1075  // For fast access (real part)
1076  const double* vj_r_pt = v[j][0].values_pt();
1077 
1078  // For fast access (imaginary part)
1079  const double* vj_c_pt = v[j][1].values_pt();
1080 
1081  // Loop over the entries in x and update them
1082  for (unsigned i = 0; i < n_row; i++)
1083  {
1084  // Update the real part of the i-th entry in x
1085  x[0][i] += (vj_r_pt[i] * y[j].real()) - (vj_c_pt[i] * y[j].imag());
1086 
1087  // Update the imaginary part of the i-th entry in x
1088  x[1][i] += (vj_c_pt[i] * y[j].real()) + (vj_r_pt[i] * y[j].imag());
1089  }
1090  } // for (unsigned j=0;j<=k;j++)
1091  } // End of update
1092 
1093  /// Helper function: Generate a plane rotation. This is done by
1094  /// finding the value of \f$ \cos(\theta) \f$ (i.e. cs) and the value of
1095  /// \f$ \sin(\theta) \f$ (i.e. sn) such that:
1096  /// \f[ \begin{bmatrix} \overline{\cos\theta} & \overline{\sin\theta} \cr -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \\ dy \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}, \f]
1097  /// where \f$ r=\sqrt{pow(|dx|,2)+pow(|dy|,2)} \f$. The values of a and b
1098  /// are given by:
1099  /// The values of dx and dy are given by:
1100  /// \f[ \cos\theta=\dfrac{dx}{\sqrt{|dx|^2+|dy|^2}}, \f]
1101  /// and
1102  /// \f[ \sin\theta=\dfrac{dy}{\sqrt{|dx|^2+|dy|^2}}. \f]
1103  /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.193.
1104  /// We also check to see that sn is always a real (nonnegative) number. See
1105  /// pp.193-194 for an explanation.
1106  void generate_plane_rotation(std::complex<double>& dx,
1107  std::complex<double>& dy,
1108  std::complex<double>& cs,
1109  std::complex<double>& sn)
1110  {
1111  // If dy=0 then we do not need to apply a rotation
1112  if (dy == 0.0)
1113  {
1114  // Using theta=0 gives cos(theta)=1
1115  cs = 1.0;
1116 
1117  // Using theta=0 gives sin(theta)=0
1118  sn = 0.0;
1119  }
1120  // If dx or dy is large using the original form of calculting cs and sn is
1121  // naive since this may overflow or underflow so instead we calculate
1122  // r=sqrt(pow(|dx|,2)+pow(|dy|,2)) using r=|dy|sqrt(1+pow(|dx|/|dy|,2)) if
1123  // |dy|>|dx| [see <A
1124  // HREF=https://en.wikipedia.org/wiki/Hypot">Hypot</A>.].
1125  else if (std::abs(dy) > std::abs(dx))
1126  {
1127  // Since |dy|>|dx| calculate the ratio |dx|/|dy|
1128  std::complex<double> temp = dx / dy;
1129 
1130  // Calculate the value of sin(theta) using:
1131  // sin(theta)=dy/sqrt(pow(|dx|,2)+pow(|dy|,2))
1132  // =(dy/|dy|)/sqrt(1+pow(|dx|/|dy|,2)).
1133  sn = (dy / std::abs(dy)) / sqrt(1.0 + pow(std::abs(temp), 2.0));
1134 
1135  // Calculate the value of cos(theta) using:
1136  // cos(theta)=dx/sqrt(pow(|dy|,2)+pow(|dx|,2))
1137  // =(dx/|dy|)/sqrt(1+pow(|dx|/|dy|,2))
1138  // =(dx/dy)*sin(theta).
1139  cs = temp * sn;
1140  }
1141  // Otherwise, we have |dx|>=|dy| so to, again, avoid overflow or underflow
1142  // calculate the values of cs and sn using the method above
1143  else
1144  {
1145  // Since |dx|>=|dy| calculate the ratio dy/dx
1146  std::complex<double> temp = dy / dx;
1147 
1148  // Calculate the value of cos(theta) using:
1149  // cos(theta)=dx/sqrt(pow(|dx|,2)+pow(|dy|,2))
1150  // =(dx/|dx|)/sqrt(1+pow(|dy|/|dx|,2)).
1151  cs = (dx / std::abs(dx)) / sqrt(1.0 + pow(std::abs(temp), 2.0));
1152 
1153  // Calculate the value of sin(theta) using:
1154  // sin(theta)=dy/sqrt(pow(|dx|,2)+pow(|dy|,2))
1155  // =(dy/|dx|)/sqrt(1+pow(|dy|/|dx|,2))
1156  // =(dy/dx)*cos(theta).
1157  sn = temp * cs;
1158  }
1159 
1160  // Set the tolerance for sin(theta)
1161  double tolerance = 1.0e-15;
1162 
1163  // Make sure sn is real and nonnegative (it should be!)
1164  if ((std::fabs(sn.imag()) > tolerance) || (sn.real() < 0))
1165  {
1166  // Create an output stream
1167  std::ostringstream error_message_stream;
1168 
1169  // Create an error message
1170  error_message_stream << "The value of sin(theta) is not real "
1171  << "and/or nonnegative. Value is: " << sn
1172  << std::endl;
1173 
1174  // Throw an error
1175  throw OomphLibError(error_message_stream.str(),
1176  OOMPH_CURRENT_FUNCTION,
1177  OOMPH_EXCEPTION_LOCATION);
1178  }
1179  } // End of generate_plane_rotation
1180 
1181  /// Helper function: Apply plane rotation. This is done using the
1182  /// update:
1183  /// \f[ \begin{bmatrix} dx \\ dy \end{bmatrix} \leftarrow \begin{bmatrix} \overline{\cos\theta} & \overline{\sin\theta} \\ -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \\ dy \end{bmatrix}. \f]
1184  /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.193.
1185  void apply_plane_rotation(std::complex<double>& dx,
1186  std::complex<double>& dy,
1187  std::complex<double>& cs,
1188  std::complex<double>& sn)
1189  {
1190  // Calculate the value of dx but don't update it yet
1191  std::complex<double> temp = std::conj(cs) * dx + std::conj(sn) * dy;
1192 
1193  // Set the value of dy
1194  dy = -sn * dx + cs * dy;
1195 
1196  // Set the value of dx using the correct values of dx and dy
1197  dx = temp;
1198  } // End of apply_plane_rotation
1199 
1200  /// Number of iterations taken
1201  unsigned Iterations;
1202 
1203  /// Vector of pointers to the real and imaginary part of the system matrix
1205 
1206  /// Boolean flag to indicate if the solve is done in re-solve mode,
1207  /// bypassing setup of matrix and preconditioner
1209 
1210  /// Boolean flag to indicate if the real and imaginary system
1211  /// matrices can be deleted
1213  };
1214 
1215  //======================================================================
1216  /// This is where the actual work is done
1217  //======================================================================
1218  template<typename MATRIX>
1220  const Vector<DoubleVector>& rhs, Vector<DoubleVector>& solution)
1221  {
1222  // Set the number of dof types (real and imaginary for this solver)
1223  unsigned n_dof_types = 2;
1224 
1225  // Get the number of dofs (note, the total number of dofs in the problem
1226  // is 2*n_row but if the constituent vectors and matrices were stored in
1227  // complex objects there would only be (n_row) rows so we use that)
1228  unsigned n_row = Matrices_storage_pt[0]->nrow();
1229 
1230  // Make sure Max_iter isn't greater than n_dof. The user cannot use this
1231  // many iterations when using Krylov subspace methods
1232  if (Max_iter > n_row)
1233  {
1234  // Create an output string stream
1235  std::ostringstream error_message_stream;
1236 
1237  // Create the error message
1238  error_message_stream << "The maximum number of iterations cannot exceed "
1239  << "the number of rows in the problem."
1240  << "\nMaximum number of iterations: " << Max_iter
1241  << "\nNumber of rows: " << n_row << std::endl;
1242 
1243  // Throw the error message
1244  throw OomphLibError(error_message_stream.str(),
1245  OOMPH_CURRENT_FUNCTION,
1246  OOMPH_EXCEPTION_LOCATION);
1247  }
1248 
1249  // Loop over the real and imaginary parts
1250  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
1251  {
1252 #ifdef PARANOID
1253  // PARANOID check that if the matrix is distributable then it should not
1254  // be then it should not be distributed
1255  if (dynamic_cast<DistributableLinearAlgebraObject*>(
1256  Matrices_storage_pt[dof_type]) != 0)
1257  {
1258  if (dynamic_cast<DistributableLinearAlgebraObject*>(
1259  Matrices_storage_pt[dof_type])
1260  ->distributed())
1261  {
1262  std::ostringstream error_message_stream;
1263  error_message_stream << "The matrix must not be distributed.";
1264  throw OomphLibError(error_message_stream.str(),
1265  OOMPH_CURRENT_FUNCTION,
1266  OOMPH_EXCEPTION_LOCATION);
1267  }
1268  }
1269  // PARANOID check that this rhs distribution is setup
1270  if (!rhs[dof_type].built())
1271  {
1272  std::ostringstream error_message_stream;
1273  error_message_stream << "The rhs vector distribution must be setup.";
1274  throw OomphLibError(error_message_stream.str(),
1275  OOMPH_CURRENT_FUNCTION,
1276  OOMPH_EXCEPTION_LOCATION);
1277  }
1278  // PARANOID check that the rhs has the right number of global rows
1279  if (rhs[dof_type].nrow() != n_row)
1280  {
1281  std::ostringstream error_message_stream;
1282  error_message_stream << "RHS does not have the same dimension as the"
1283  << " linear system";
1284  throw OomphLibError(error_message_stream.str(),
1285  OOMPH_CURRENT_FUNCTION,
1286  OOMPH_EXCEPTION_LOCATION);
1287  }
1288  // PARANOID check that the rhs is not distributed
1289  if (rhs[dof_type].distribution_pt()->distributed())
1290  {
1291  std::ostringstream error_message_stream;
1292  error_message_stream << "The rhs vector must not be distributed.";
1293  throw OomphLibError(error_message_stream.str(),
1294  OOMPH_CURRENT_FUNCTION,
1295  OOMPH_EXCEPTION_LOCATION);
1296  }
1297  // PARANOID check that if the result is setup it matches the distribution
1298  // of the rhs
1299  if (solution[dof_type].built())
1300  {
1301  if (!(*rhs[dof_type].distribution_pt() ==
1302  *solution[dof_type].distribution_pt()))
1303  {
1304  std::ostringstream error_message_stream;
1305  error_message_stream << "If the result distribution is setup then it "
1306  << "must be the same as the rhs distribution";
1307  throw OomphLibError(error_message_stream.str(),
1308  OOMPH_CURRENT_FUNCTION,
1309  OOMPH_EXCEPTION_LOCATION);
1310  }
1311  } // if (solution[dof_type].built())
1312 #endif
1313  // Set up the solution distribution if it's not already distributed
1314  if (!solution[dof_type].built())
1315  {
1316  // Build the distribution
1317  solution[dof_type].build(this->distribution_pt(), 0.0);
1318  }
1319  // Otherwise initialise all entries to zero
1320  else
1321  {
1322  // Initialise the entries in the k-th vector in solution to zero
1323  solution[dof_type].initialise(0.0);
1324  }
1325  } // for (unsigned dof_type=0;dof_type<n_dof_types;dof_type++)
1326 
1327  // Start the solver timer
1328  double t_start = TimingHelpers::timer();
1329 
1330  // Storage for the relative residual
1331  double resid;
1332 
1333  // Initialise vectors (since these are not distributed vectors we template
1334  // one vector by the type std::complex<double> instead of using two vectors,
1335  // each templated by the type double
1336 
1337  // Vector, s, used in the minimisation problem: J(y)=min||s-R_m*y||
1338  // [see Saad Y."Iterative methods for sparse linear systems", p.176.]
1339  Vector<std::complex<double>> s(n_row + 1, std::complex<double>(0.0, 0.0));
1340 
1341  // Vector to store the value of cos(theta) when using the Givens rotation
1342  Vector<std::complex<double>> cs(n_row + 1, std::complex<double>(0.0, 0.0));
1343 
1344  // Vector to store the value of sin(theta) when using the Givens rotation
1345  Vector<std::complex<double>> sn(n_row + 1, std::complex<double>(0.0, 0.0));
1346 
1347  // Create a vector of DoubleVectors (this is a distributed vector so we have
1348  // to create two separate DoubleVector objects to cope with the arithmetic)
1349  Vector<DoubleVector> w(n_dof_types);
1350 
1351  // Build the distribution of both vectors
1352  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
1353  {
1354  // Build the distribution of the k-th constituent vector
1355  w[dof_type].build(this->distribution_pt(), 0.0);
1356  }
1357 
1358  // Create a vector of DoubleVectors to store the RHS of b-Jx=Mr. We assume
1359  // both x=0 and that a preconditioner is not applied by which we deduce b=r
1360  Vector<DoubleVector> r(n_dof_types);
1361 
1362  // Build the distribution of both vectors
1363  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
1364  {
1365  // Build the distribution of the k-th constituent vector
1366  r[dof_type].build(this->distribution_pt(), 0.0);
1367  }
1368 
1369  // Store the value of b (the RHS vector) in r
1370  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
1371  {
1372  // Build the distribution of the k-th constituent vector
1373  r[dof_type] = rhs[dof_type];
1374  }
1375 
1376  // Calculate the norm of the real part of r
1377  double norm_r = r[0].norm();
1378 
1379  // Calculate the norm of the imaginary part of r
1380  double norm_c = r[1].norm();
1381 
1382  // Compute norm(r)
1383  double normb = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
1384 
1385  // Set the value of beta (the initial residual)
1386  double beta = normb;
1387 
1388  // Compute the initial relative residual. If the entries of the RHS vector
1389  // are all zero then set normb equal to one. This is because we divide the
1390  // value of the norm at later stages by normb and dividing by zero is not
1391  // definied
1392  if (normb == 0.0)
1393  {
1394  // Set the value of normb
1395  normb = 1.0;
1396  }
1397 
1398  // Calculate the ratio between the initial norm and the current norm.
1399  // Since we haven't completed an iteration yet this will simply be one
1400  // unless normb was zero, in which case resid will have value zero
1401  resid = beta / normb;
1402 
1403  // If required, will document convergence history to screen or file (if
1404  // stream open)
1405  if (Doc_convergence_history)
1406  {
1407  // If an output file which is open isn't provided then output to screen
1408  if (!Output_file_stream.is_open())
1409  {
1410  // Output the residual value to the screen
1411  oomph_info << 0 << " " << resid << std::endl;
1412  }
1413  // Otherwise, output to file
1414  else
1415  {
1416  // Output the residual value to file
1417  Output_file_stream << 0 << " " << resid << std::endl;
1418  }
1419  } // if (Doc_convergence_history)
1420 
1421  // If the GMRES algorithm converges immediately
1422  if (resid <= Tolerance)
1423  {
1424  // If time documentation is enabled
1425  if (Doc_time)
1426  {
1427  // Notify the user
1428  oomph_info << "GMRES converged immediately. Normalised residual norm: "
1429  << resid << std::endl;
1430  }
1431 
1432  // Finish running the solver
1433  return;
1434  } // if (resid<=Tolerance)
1435 
1436  // Initialise a vector of orthogonal basis vectors
1438 
1439  // Resize the number of vectors needed
1440  v.resize(n_row + 1);
1441 
1442  // Resize each Vector of DoubleVectors to store the real and imaginary
1443  // part of a given vector
1444  for (unsigned dof_type = 0; dof_type < n_row + 1; dof_type++)
1445  {
1446  // Create two DoubleVector objects in each Vector
1447  v[dof_type].resize(n_dof_types);
1448  }
1449 
1450  // Initialise the upper hessenberg matrix. Since we are not using
1451  // distributed vectors here, the algebra is best done using entries
1452  // of the type std::complex<double>. NOTE: For implementation purposes
1453  // the upper Hessenberg matrix indices are swapped so the matrix is
1454  // effectively transposed
1455  Vector<Vector<std::complex<double>>> hessenberg(n_row + 1);
1456 
1457  // Build the zeroth basis vector
1458  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
1459  {
1460  // Build the k-th part of the zeroth vector. Here k=0 and k=1 correspond
1461  // to the real and imaginary part of the zeroth vector, respectively
1462  v[0][dof_type].build(this->distribution_pt(), 0.0);
1463  }
1464 
1465  // Loop over the real and imaginary parts of v
1466  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
1467  {
1468  // For fast access
1469  double* v0_pt = v[0][dof_type].values_pt();
1470 
1471  // For fast access
1472  const double* r_pt = r[dof_type].values_pt();
1473 
1474  // Set the zeroth basis vector v[0] to r/beta
1475  for (unsigned i = 0; i < n_row; i++)
1476  {
1477  // Assign the i-th entry of the zeroth basis vector
1478  v0_pt[i] = r_pt[i] / beta;
1479  }
1480  } // for (unsigned k=0;k<n_dof_types;k++)
1481 
1482  // Set the first entry in the minimisation problem RHS vector (is meant
1483  // to the vector beta*e_1 initially, where e_1 is the unit vector with
1484  // one in its first entry)
1485  s[0] = beta;
1486 
1487  // Compute the next step of the iterative scheme
1488  for (unsigned j = 0; j < Max_iter; j++)
1489  {
1490  // Resize the next column of the upper hessenberg matrix
1491  hessenberg[j].resize(j + 2, std::complex<double>(0.0, 0.0));
1492 
1493  // Calculate w=J*v_j. Note, we cannot use inbuilt complex matrix algebra
1494  // here as we're using distributed vectors
1495  complex_matrix_multiplication(Matrices_storage_pt, v[j], w);
1496 
1497  // For fast access
1498  double* w_r_pt = w[0].values_pt();
1499 
1500  // For fast access
1501  double* w_c_pt = w[1].values_pt();
1502 
1503  // Loop over all of the entries on and above the principal subdiagonal of
1504  // the Hessenberg matrix in the j-th column (remembering that
1505  // the indices of the upper Hessenberg matrix are swapped for the purpose
1506  // of implementation)
1507  for (unsigned i = 0; i < j + 1; i++)
1508  {
1509  // For fast access
1510  const double* vi_r_pt = v[i][0].values_pt();
1511 
1512  // For fast access
1513  const double* vi_c_pt = v[i][1].values_pt();
1514 
1515  // Loop over the entries of v and w
1516  for (unsigned k = 0; k < n_row; k++)
1517  {
1518  // Store the appropriate entry in v as a complex value
1519  std::complex<double> complex_v(vi_r_pt[k], vi_c_pt[k]);
1520 
1521  // Store the appropriate entry in w as a complex value
1522  std::complex<double> complex_w(w_r_pt[k], w_c_pt[k]);
1523 
1524  // Update the value of H(i,j) noting we're computing a complex
1525  // inner product here
1526  hessenberg[j][i] += std::conj(complex_v) * complex_w;
1527  }
1528 
1529  // Orthonormalise w against all previous orthogonal vectors, v_i by
1530  // looping over its entries and updating them
1531  for (unsigned k = 0; k < n_row; k++)
1532  {
1533  // Update the real part of the k-th entry of w
1534  w_r_pt[k] -= (hessenberg[j][i].real() * vi_r_pt[k] -
1535  hessenberg[j][i].imag() * vi_c_pt[k]);
1536 
1537  // Update the imaginary part of the k-th entry of w
1538  w_c_pt[k] -= (hessenberg[j][i].real() * vi_c_pt[k] +
1539  hessenberg[j][i].imag() * vi_r_pt[k]);
1540  }
1541  } // for (unsigned i=0;i<j+1;i++)
1542 
1543  // Calculate the norm of the real part of w
1544  norm_r = w[0].norm();
1545 
1546  // Calculate the norm of the imaginary part of w
1547  norm_c = w[1].norm();
1548 
1549  // Calculate the norm of the vector w using norm_r and norm_c and assign
1550  // its value to the appropriate entry in the Hessenberg matrix
1551  hessenberg[j][j + 1] = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
1552 
1553  // Build the real part of the next orthogonal vector
1554  v[j + 1][0].build(this->distribution_pt(), 0.0);
1555 
1556  // Build the imaginary part of the next orthogonal vector
1557  v[j + 1][1].build(this->distribution_pt(), 0.0);
1558 
1559  // Check if the value of hessenberg[j][j+1] is zero. If it
1560  // isn't then we update the next entries in v
1561  if (hessenberg[j][j + 1] != 0.0)
1562  {
1563  // For fast access
1564  double* v_r_pt = v[j + 1][0].values_pt();
1565 
1566  // For fast access
1567  double* v_c_pt = v[j + 1][1].values_pt();
1568 
1569  // For fast access
1570  const double* w_r_pt = w[0].values_pt();
1571 
1572  // For fast access
1573  const double* w_c_pt = w[1].values_pt();
1574 
1575  // Notice, the value of H(j,j+1), as calculated above, is clearly a real
1576  // number. As such, calculating the division
1577  // v_{j+1}=w_{j}/h_{j+1,j},
1578  // here is simple, i.e. we don't need to worry about cross terms in the
1579  // algebra. To avoid computing h_{j+1,j} several times we precompute it
1580  double h_subdiag_val = hessenberg[j][j + 1].real();
1581 
1582  // Loop over the entries of the new orthogonal vector and set its values
1583  for (unsigned k = 0; k < n_row; k++)
1584  {
1585  // The i-th entry of the real component is given by
1586  v_r_pt[k] = w_r_pt[k] / h_subdiag_val;
1587 
1588  // Similarly, the i-th entry of the imaginary component is given by
1589  v_c_pt[k] = w_c_pt[k] / h_subdiag_val;
1590  }
1591  }
1592  // Otherwise, we have to jump to the next part of the algorithm; if
1593  // the value of hessenberg[j][j+1] is zero then the norm of the latest
1594  // orthogonal vector is zero. This is only possible if the entries
1595  // in w are all zero. As a result, the Krylov space of A and r_0 has
1596  // been spanned by the previously calculated orthogonal vectors
1597  else
1598  {
1599  // Book says "Set m=j and jump to step 11" (p.172)...
1600  // Do something here!
1601  oomph_info << "Subdiagonal Hessenberg entry is zero."
1602  << "Do something here..." << std::endl;
1603  } // if (hessenberg[j][j+1]!=0.0)
1604 
1605  // Loop over the entries in the Hessenberg matrix and calculate the
1606  // entries of the Givens rotation matrices
1607  for (unsigned k = 0; k < j; k++)
1608  {
1609  // Apply the plane rotation to all of the previous entries in the
1610  // (j)-th column (remembering the indexing is reversed)
1611  apply_plane_rotation(
1612  hessenberg[j][k], hessenberg[j][k + 1], cs[k], sn[k]);
1613  }
1614 
1615  // Now calculate the entries of the latest Givens rotation matrix
1616  generate_plane_rotation(
1617  hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
1618 
1619  // Apply the plane rotation using the newly calculated entries
1620  apply_plane_rotation(
1621  hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
1622 
1623  // Apply a plane rotation to the corresponding entry in the vector
1624  // s used in the minimisation problem, J(y)=min||s-R_m*y||
1625  apply_plane_rotation(s[j], s[j + 1], cs[j], sn[j]);
1626 
1627  // Compute current residual using equation (6.42) in Saad Y, "Iterative
1628  // methods for sparse linear systems", p.177.]. Note, since s has complex
1629  // entries we have to use std::abs instead of std::fabs
1630  beta = std::abs(s[j + 1]);
1631 
1632  // Compute the relative residual
1633  resid = beta / normb;
1634 
1635  // If required will document convergence history to screen or file (if
1636  // stream open)
1637  if (Doc_convergence_history)
1638  {
1639  // If an output file which is open isn't provided then output to screen
1640  if (!Output_file_stream.is_open())
1641  {
1642  // Output the residual value to the screen
1643  oomph_info << j + 1 << " " << resid << std::endl;
1644  }
1645  // Otherwise, output to file
1646  else
1647  {
1648  // Output the residual value to file
1649  Output_file_stream << j + 1 << " " << resid << std::endl;
1650  }
1651  } // if (Doc_convergence_history)
1652 
1653  // If the required tolerance has been met
1654  if (resid < Tolerance)
1655  {
1656  // Store the number of iterations taken
1657  Iterations = j + 1;
1658 
1659  // Update the result vector using the result, x=x_0+V_m*y (where V_m
1660  // is given by v here)
1661  update(j, hessenberg, s, v, solution);
1662 
1663  // If time documentation was enabled
1664  if (Doc_time)
1665  {
1666  // Output the current normalised residual norm
1667  oomph_info << "\nGMRES converged (1). Normalised residual norm: "
1668  << resid << std::endl;
1669 
1670  // Output the number of iterations it took for convergence
1671  oomph_info << "Number of iterations to convergence: " << j + 1 << "\n"
1672  << std::endl;
1673  }
1674 
1675  // Stop the timer
1676  double t_end = TimingHelpers::timer();
1677 
1678  // Calculate the time taken to calculate the solution
1679  Solution_time = t_end - t_start;
1680 
1681  // If time documentation was enabled
1682  if (Doc_time)
1683  {
1684  // Output the time taken to solve the problem using GMRES
1685  oomph_info << "Time for solve with GMRES [sec]: " << Solution_time
1686  << std::endl;
1687  }
1688 
1689  // As we've met the tolerance for the solver and everything that should
1690  // be documented, has been, finish using the solver
1691  return;
1692  } // if (resid<Tolerance)
1693  } // for (unsigned j=0;j<Max_iter;j++)
1694 
1695  // Store the number of iterations taken
1696  Iterations = Max_iter;
1697 
1698  // Only update if we actually did something
1699  if (Max_iter > 0)
1700  {
1701  // Update the result vector using the result, x=x_0+V_m*y (where V_m
1702  // is given by v here)
1703  update(Max_iter - 1, hessenberg, s, v, solution);
1704  }
1705 
1706  // Stop the timer
1707  double t_end = TimingHelpers::timer();
1708 
1709  // Calculate the time taken to calculate the solution
1710  Solution_time = t_end - t_start;
1711 
1712  // If time documentation was enabled
1713  if (Doc_time)
1714  {
1715  // Output the time taken to solve the problem using GMRES
1716  oomph_info << "Time for solve with GMRES [sec]: " << Solution_time
1717  << std::endl;
1718  }
1719 
1720  // Finish using the solver
1721  return;
1722  } // End of complex_solve_helper
1723 
1724 
1725  /// ////////////////////////////////////////////////////////////////////
1726  /// ////////////////////////////////////////////////////////////////////
1727  /// ////////////////////////////////////////////////////////////////////
1728 
1729 
1730  //======================================================================
1731  /// The GMRES method for the Helmholtz solver.
1732  //======================================================================
1733  template<typename MATRIX>
1735  public BlockPreconditioner<MATRIX>
1736  {
1737  public:
1738  /// Constructor
1741  Iterations(0),
1742  Resolving(false),
1743  Matrix_can_be_deleted(true)
1744  {
1745  Preconditioner_LHS = true;
1746  }
1747 
1748  /// Destructor (cleanup storage)
1750  {
1751  clean_up_memory();
1752  }
1753 
1754  /// Broken copy constructor
1756 
1757  /// Broken assignment operator
1758  void operator=(const HelmholtzGMRESMG&) = delete;
1759 
1760  /// Overload disable resolve so that it cleans up memory too
1762  {
1764  clean_up_memory();
1765  }
1766 
1767  /// Implementation of the pure virtual base class function. The
1768  /// function has been broken because this is meant to be used as a linear
1769  /// solver
1771  {
1772  // Open an output stream
1773  std::ostringstream error_message_stream;
1774 
1775  // Create an error message
1776  error_message_stream << "Preconditioner_solve(...) is broken. "
1777  << "HelmholtzGMRESMG is only meant to be used as "
1778  << "a linear solver.\n";
1779 
1780  // Throw the error message
1781  throw OomphLibError(error_message_stream.str(),
1782  OOMPH_CURRENT_FUNCTION,
1783  OOMPH_EXCEPTION_LOCATION);
1784  }
1785 
1786  /// Implementation of the pure virtual base class function. This
1787  /// accompanies the preconditioner_solve function and so is also broken
1788  void setup()
1789  {
1790  // Open an output stream
1791  std::ostringstream error_message_stream;
1792 
1793  // Create an error message
1794  error_message_stream << "This function is broken. HelmholtzGMRESMG is "
1795  << "only meant to be used as a linear solver.\n";
1796 
1797  // Throw the error message
1798  throw OomphLibError(error_message_stream.str(),
1799  OOMPH_CURRENT_FUNCTION,
1800  OOMPH_EXCEPTION_LOCATION);
1801  }
1802 
1803  /// Solver: Takes pointer to problem and returns the results vector
1804  /// which contains the solution of the linear system defined by
1805  /// the problem's fully assembled Jacobian and residual vector.
1806  void solve(Problem* const& problem_pt, DoubleVector& result)
1807  {
1808 #ifdef OOMPH_HAS_MPI
1809  // Make sure that this is running in serial. Can't guarantee it'll
1810  // work when the problem is distributed over several processors
1811  if (MPI_Helpers::communicator_pt()->nproc() > 1)
1812  {
1813  // Throw a warning
1814  OomphLibWarning("Can't guarantee the MG solver will work in parallel!",
1815  OOMPH_CURRENT_FUNCTION,
1816  OOMPH_EXCEPTION_LOCATION);
1817  }
1818 #endif
1819 
1820  // Find # of degrees of freedom (variables)
1821  unsigned n_dof = problem_pt->ndof();
1822 
1823  // Initialise timer
1824  double t_start = TimingHelpers::timer();
1825 
1826  // We're not re-solving
1827  Resolving = false;
1828 
1829  // Get rid of any previously stored data
1830  clean_up_memory();
1831 
1832  // Grab the communicator from the MGProblem object and assign it
1833  this->set_comm_pt(problem_pt->communicator_pt());
1834 
1835  // Setup the distribution
1837  problem_pt->communicator_pt(), n_dof, false);
1838 
1839  // Build the internal distribution in this way because both the
1840  // IterativeLinearSolver and BlockPreconditioner class have base-
1841  // class subobjects of type oomph::DistributableLinearAlgebraObject
1843 
1844  // Get Jacobian matrix in format specified by template parameter
1845  // and nonlinear residual vector
1846  MATRIX* matrix_pt = new MATRIX;
1847  DoubleVector f;
1848  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
1849  {
1850  if (dynamic_cast<CRDoubleMatrix*>(matrix_pt) != 0)
1851  {
1852  dynamic_cast<CRDoubleMatrix*>(matrix_pt)->build(
1855  }
1856  }
1857 
1858  // Get the Jacobian and residuals vector
1859  problem_pt->get_jacobian(f, *matrix_pt);
1860 
1861  // We've made the matrix, we can delete it...
1862  Matrix_can_be_deleted = true;
1863 
1864  // Replace the current matrix used in Preconditioner by the new matrix
1865  this->set_matrix_pt(matrix_pt);
1866 
1867  // The preconditioner works with one mesh; set it! Since we only use
1868  // the block preconditioner on the finest level, we use the mesh from
1869  // that level
1870  this->set_nmesh(1);
1871 
1872  // Elements in actual pml layer are trivially wrapped versions of
1873  // their bulk counterparts. Technically they are different elements
1874  // so we have to allow different element types
1875  bool allow_different_element_types_in_mesh = true;
1876  this->set_mesh(
1877  0, problem_pt->mesh_pt(), allow_different_element_types_in_mesh);
1878 
1879  // Set up the generic block look up scheme
1880  this->block_setup();
1881 
1882  // Extract the number of blocks.
1883  unsigned nblock_types = this->nblock_types();
1884 
1885 #ifdef PARANOID
1886  // PARANOID check - there must only be two block types
1887  if (nblock_types != 2)
1888  {
1889  // Create the error message
1890  std::stringstream tmp;
1891  tmp << "There are supposed to be two block types.\nYours has "
1892  << nblock_types << std::endl;
1893 
1894  // Throw an error
1895  throw OomphLibError(
1896  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1897  }
1898 #endif
1899 
1900  // Resize the storage for the system matrices
1901  Matrices_storage_pt.resize(2, 0);
1902 
1903  // Loop over the rows of the block matrix
1904  for (unsigned i = 0; i < nblock_types; i++)
1905  {
1906  // Fix the column index
1907  unsigned j = 0;
1908 
1909  // Create new CRDoubleMatrix objects
1911 
1912  // Extract the required blocks, i.e. the first column
1913  this->get_block(i, j, *Matrices_storage_pt[i]);
1914  }
1915 
1916  // Doc time for setup
1917  double t_end = TimingHelpers::timer();
1918  Jacobian_setup_time = t_end - t_start;
1919 
1920  if (Doc_time)
1921  {
1922  oomph_info << "\nTime for setup of block Jacobian [sec]: "
1923  << Jacobian_setup_time << std::endl;
1924  }
1925 
1926  // Call linear algebra-style solver
1927  // If the result distribution is wrong, then redistribute
1928  // before the solve and return to original distribution
1929  // afterwards
1930  if ((!(*result.distribution_pt() ==
1932  (result.built()))
1933  {
1934  // Make a distribution object
1935  LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
1936 
1937  // Build the result vector distribution
1939 
1940  // Solve the problem
1941  solve_helper(matrix_pt, f, result);
1942 
1943  // Redistribute the vector
1944  result.redistribute(&temp_global_dist);
1945  }
1946  // Otherwise just solve
1947  else
1948  {
1949  // Solve
1950  solve_helper(matrix_pt, f, result);
1951  }
1952 
1953  // Kill matrix unless it's still required for resolve
1954  if (!Enable_resolve)
1955  {
1956  // Clean up anything left in memory
1957  clean_up_memory();
1958  }
1959  } // End of solve
1960 
1961  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
1962  /// vector and returns the solution of the linear system.
1964  const DoubleVector& rhs,
1965  DoubleVector& solution)
1966  {
1967  // Open an output stream
1968  std::ostringstream error_message_stream;
1969 
1970  // Create an error message
1971  error_message_stream
1972  << "This function is broken. The block preconditioner "
1973  << "needs access to the underlying mesh.\n";
1974 
1975  // Throw the error message
1976  throw OomphLibError(error_message_stream.str(),
1977  OOMPH_CURRENT_FUNCTION,
1978  OOMPH_EXCEPTION_LOCATION);
1979  }
1980 
1981 
1982  /// Linear-algebra-type solver: Takes pointer to a matrix
1983  /// and rhs vector and returns the solution of the linear system
1984  /// Call the broken base-class version. If you want this, please
1985  /// implement it
1987  const Vector<double>& rhs,
1988  Vector<double>& result)
1989  {
1990  LinearSolver::solve(matrix_pt, rhs, result);
1991  }
1992 
1993  /// Re-solve the system defined by the last assembled Jacobian
1994  /// and the rhs vector specified here. Solution is returned in the
1995  /// vector result.
1996  void resolve(const DoubleVector& rhs, DoubleVector& result)
1997  {
1998  // We are re-solving
1999  Resolving = true;
2000 
2001 #ifdef PARANOID
2002  if ((Matrices_storage_pt[0] == 0) || (Matrices_storage_pt[1] == 0))
2003  {
2004  throw OomphLibError("No matrix was stored -- cannot re-solve",
2005  OOMPH_CURRENT_FUNCTION,
2006  OOMPH_EXCEPTION_LOCATION);
2007  }
2008 #endif
2009 
2010  // Set up a dummy matrix. As we're resolving this won't be used in
2011  // solve_helper but we need to pass a matrix in to fill the input.
2012  // The matrices used in the calculations have already been stored
2014 
2015  // Call the helper function
2016  solve_helper(matrix_pt, rhs, result);
2017 
2018  // Delete the matrix
2019  delete matrix_pt;
2020 
2021  // Make it a null pointer
2022  matrix_pt = 0;
2023 
2024  // Reset re-solving flag
2025  Resolving = false;
2026  }
2027 
2028  /// Number of iterations taken
2029  unsigned iterations() const
2030  {
2031  return Iterations;
2032  }
2033 
2034  /// Set left preconditioning (the default)
2036  {
2037  Preconditioner_LHS = true;
2038  }
2039 
2040  /// Enable right preconditioning
2042  {
2043  Preconditioner_LHS = false;
2044  }
2045 
2046  protected:
2047  /// General interface to solve function
2049  const DoubleVector& rhs,
2050  DoubleVector& solution);
2051 
2052  /// Cleanup data that's stored for resolve (if any has been stored)
2054  {
2055  // If the matrix storage has been resized
2056  if (Matrices_storage_pt.size() > 0)
2057  {
2058  // If the real matrix pointer isn't null AND we're allowed to delete
2059  // the matrix which is only when we create the matrix ourselves
2060  if ((Matrices_storage_pt[0] != 0) && (Matrix_can_be_deleted))
2061  {
2062  // Delete the matrix
2063  delete Matrices_storage_pt[0];
2064 
2065  // Assign the associated pointer the value NULL
2066  Matrices_storage_pt[0] = 0;
2067  }
2068 
2069  // If the real matrix pointer isn't null AND we're allowed to delete
2070  // the matrix which is only when we create the matrix ourselves
2071  if ((Matrices_storage_pt[1] != 0) && (Matrix_can_be_deleted))
2072  {
2073  // Delete the matrix
2074  delete Matrices_storage_pt[1];
2075 
2076  // Assign the associated pointer the value NULL
2077  Matrices_storage_pt[1] = 0;
2078  }
2079  }
2080  } // End of clean_up_memory
2081 
2082  /// Helper function to calculate a complex matrix-vector product.
2083  /// Assumes the matrix has been provided as a Vector of length two; the
2084  /// first entry containing the real part of the system matrix and the
2085  /// second entry containing the imaginary part
2087  Vector<CRDoubleMatrix*> const matrices_pt,
2088  const Vector<DoubleVector>& x,
2089  Vector<DoubleVector>& soln)
2090  {
2091 #ifdef PARANOID
2092  // PARANOID check - Make sure the input matrix has the right size
2093  if (matrices_pt.size() != 2)
2094  {
2095  // Create an output stream
2096  std::ostringstream error_message_stream;
2097 
2098  // Create the error message
2099  error_message_stream << "Can only deal with two matrices. You have "
2100  << matrices_pt.size() << " matrices." << std::endl;
2101 
2102  // Throw an error
2103  throw OomphLibError(error_message_stream.str(),
2104  OOMPH_CURRENT_FUNCTION,
2105  OOMPH_EXCEPTION_LOCATION);
2106  }
2107  // PARANOID check - Make sure the vector x has the right size
2108  if (x.size() != 2)
2109  {
2110  // Create an output stream
2111  std::ostringstream error_message_stream;
2112 
2113  // Create the error message
2114  error_message_stream
2115  << "Can only deal with two input vectors. You have " << x.size()
2116  << " vectors." << std::endl;
2117 
2118  // Throw an error
2119  throw OomphLibError(error_message_stream.str(),
2120  OOMPH_CURRENT_FUNCTION,
2121  OOMPH_EXCEPTION_LOCATION);
2122  }
2123  // PARANOID check - Make sure the vector soln has the right size
2124  if (soln.size() != 2)
2125  {
2126  // Create an output stream
2127  std::ostringstream error_message_stream;
2128 
2129  // Create the error message
2130  error_message_stream
2131  << "Can only deal with two output vectors. You have " << soln.size()
2132  << " output vectors." << std::endl;
2133 
2134  // Throw an error
2135  throw OomphLibError(error_message_stream.str(),
2136  OOMPH_CURRENT_FUNCTION,
2137  OOMPH_EXCEPTION_LOCATION);
2138  }
2139 #endif
2140 
2141  // NOTE: We assume all vectors have been distributed at this point but
2142  // code can be written at a later time to build the vectors if they're
2143  // not already built.
2144 
2145  //-----------------------------------------------------------------------
2146  // Suppose we have a complex matrix, A, and two complex vectors, x and
2147  // soln. We wish to compute the product A*x=soln (note, * does not mean
2148  // we are using complex conjugates here, it is simply used to indicate
2149  // a multiplication). To do this we must make use of the fact that we
2150  // possess the real and imaginary separately. As a result, it is computed
2151  // using:
2152  // soln = A*x,
2153  // = (A_r + i*A_c)*(x_r + i*x_c),
2154  // = [A_r*x_r - A_c*x_c] + i*[A_r*x_c + A_c*x_r],
2155  // ==> real(soln) = A_r*x_r - A_c*x_c,
2156  // & imag(soln) = A_r*x_c + A_c*x_r,
2157  // where the subscripts _r and _c are used to identify the real and
2158  // imaginary part, respectively.
2159  //-----------------------------------------------------------------------
2160 
2161  // Store the value of A_r*x_r in the real part of soln
2162  matrices_pt[0]->multiply(x[0], soln[0]);
2163 
2164  // Store the value of A_r*x_c in the imaginary part of soln
2165  matrices_pt[0]->multiply(x[1], soln[1]);
2166 
2167  // Create a temporary vector
2169 
2170  // Calculate the value of A_c*x_c
2171  matrices_pt[1]->multiply(x[1], temp);
2172 
2173  // Subtract the value of temp from soln[0] to get the real part of soln
2174  soln[0] -= temp;
2175 
2176  // Calculate the value of A_c*x_r
2177  matrices_pt[1]->multiply(x[0], temp);
2178 
2179  // Add the value of temp to soln[1] to get the imaginary part of soln
2180  soln[1] += temp;
2181  } // End of complex_matrix_multiplication
2182 
2183  /// Helper function to update the result vector
2184  void update(const unsigned& k,
2185  const Vector<Vector<std::complex<double>>>& hessenberg,
2186  const Vector<std::complex<double>>& s,
2187  const Vector<Vector<DoubleVector>>& v,
2189  {
2190  // Make a local copy of s
2192 
2193  //-----------------------------------------------------------------
2194  // The Hessenberg matrix should be an upper triangular matrix at
2195  // this point (although from its storage it would appear to be a
2196  // lower triangular matrix since the indexing has been reversed)
2197  // so finding the minimiser of J(y)=min||s-R_m*y|| where R_m is
2198  // the matrix R in the QR factorisation of the Hessenberg matrix.
2199  // Therefore, to obtain y we simply need to use a backwards
2200  // substitution. Note: The implementation here may appear to be
2201  // somewhat confusing as the indexing in the Hessenberg matrix is
2202  // reversed. This implementation of a backwards substitution does
2203  // not run along the columns of the triangular matrix but rather
2204  // up the rows.
2205  //-----------------------------------------------------------------
2206 
2207  // The outer loop is a loop over the columns of the Hessenberg matrix
2208  // since the indexing is reversed
2209  for (int i = int(k); i >= 0; i--)
2210  {
2211  // Divide the i-th entry of y by the i-th diagonal entry of H
2212  y[i] /= hessenberg[i][i];
2213 
2214  // The inner loop is a loop over the rows of the Hessenberg matrix
2215  for (int j = i - 1; j >= 0; j--)
2216  {
2217  // Update the j-th entry of y
2218  y[j] -= hessenberg[i][j] * y[i];
2219  }
2220  } // for (int i=int(k);i>=0;i--)
2221 
2222  // Calculate the number of entries in x (simply use the real part as
2223  // both the real and imaginary part should have the same length)
2224  unsigned n_row = x[0].nrow();
2225 
2226  // Build a temporary vector with entries initialised to 0.0
2227  Vector<DoubleVector> block_z(2);
2228 
2229  // Build a temporary vector with entries initialised to 0.0
2230  Vector<DoubleVector> block_temp(2);
2231 
2232  // Build the distributions
2233  for (unsigned dof_type = 0; dof_type < 2; dof_type++)
2234  {
2235  // Build the (dof_type)-th vector
2236  block_z[dof_type].build(x[0].distribution_pt(), 0.0);
2237 
2238  // Build the (dof_type)-th vector
2239  block_temp[dof_type].build(x[0].distribution_pt(), 0.0);
2240  }
2241 
2242  // Get access to the underlying values
2243  double* block_temp_r_pt = block_temp[0].values_pt();
2244 
2245  // Get access to the underlying values
2246  double* block_temp_c_pt = block_temp[1].values_pt();
2247 
2248  // Calculate x=Vy
2249  for (unsigned j = 0; j <= k; j++)
2250  {
2251  // Get access to j-th column of Z_m
2252  const double* vj_r_pt = v[j][0].values_pt();
2253 
2254  // Get access to j-th column of Z_m
2255  const double* vj_c_pt = v[j][1].values_pt();
2256 
2257  // Loop over the entries in x and update them
2258  for (unsigned i = 0; i < n_row; i++)
2259  {
2260  // Update the real part of the i-th entry in x
2261  block_temp_r_pt[i] +=
2262  (vj_r_pt[i] * y[j].real()) - (vj_c_pt[i] * y[j].imag());
2263 
2264  // Update the imaginary part of the i-th entry in x
2265  block_temp_c_pt[i] +=
2266  (vj_c_pt[i] * y[j].real()) + (vj_r_pt[i] * y[j].imag());
2267  }
2268  } // for (unsigned j=0;j<=k;j++)
2269 
2270  // If we're using LHS preconditioning
2271  if (Preconditioner_LHS)
2272  {
2273  // Since we're using LHS preconditioning the preconditioner is applied
2274  // to the matrix and RHS vector so we simply update the value of x
2275  for (unsigned dof_type = 0; dof_type < 2; dof_type++)
2276  {
2277  // Update
2278  x[dof_type] += block_temp[dof_type];
2279  }
2280  }
2281  // If we're using RHS preconditioning
2282  else
2283  {
2284  // Create a temporary vector
2286 
2287  // Copy block vectors block_temp back to temp
2288  this->return_block_vectors(block_temp, temp);
2289 
2290  // Create a temporary vector
2292 
2293  // Copy block vectors block_z back to z
2294  this->return_block_vectors(block_z, z);
2295 
2296  // Since we're using RHS preconditioning the preconditioner is applied
2297  // to the solution vector
2299 
2300  // Split up the solution vector into DoubleVectors, whose entries are
2301  // arranged to match the matrix blocks and assign it
2302  this->get_block_vectors(z, block_z);
2303 
2304  // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
2305  // sparse linear systems", p.284]
2306  for (unsigned dof_type = 0; dof_type < 2; dof_type++)
2307  {
2308  // Update
2309  x[dof_type] += block_z[dof_type];
2310  }
2311  } // if(Preconditioner_LHS) else
2312  } // End of update
2313 
2314  /// Helper function: Generate a plane rotation. This is done by
2315  /// finding the value of \f$ \cos(\theta) \f$ (i.e. cs) and the value of
2316  /// \f$ \sin(\theta) \f$ (i.e. sn) such that:
2317  /// \f[ \begin{bmatrix} \overline{\cos\theta} & \overline{\sin\theta} \\ -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \\ dy \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}, \f]
2318  /// where \f$ r=\sqrt{pow(|dx|,2)+pow(|dy|,2)} \f$. The values of a and b
2319  /// are given by:
2320  /// The values of dx and dy are given by:
2321  /// \f[ \cos\theta=\dfrac{dx}{\sqrt{|dx|^2+|dy|^2}}, \f]
2322  /// and
2323  /// \f[ \sin\theta=\dfrac{dy}{\sqrt{|dx|^2+|dy|^2}}. \f]
2324  /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.193.
2325  /// We also check to see that sn is always a real (nonnegative) number. See
2326  /// pp.193-194 for an explanation.
2327  void generate_plane_rotation(std::complex<double>& dx,
2328  std::complex<double>& dy,
2329  std::complex<double>& cs,
2330  std::complex<double>& sn)
2331  {
2332  // If dy=0 then we do not need to apply a rotation
2333  if (dy == 0.0)
2334  {
2335  // Using theta=0 gives cos(theta)=1
2336  cs = 1.0;
2337 
2338  // Using theta=0 gives sin(theta)=0
2339  sn = 0.0;
2340  }
2341  // If dx or dy is large using the original form of calculting cs and sn is
2342  // naive since this may overflow or underflow so instead we calculate
2343  // r=sqrt(pow(|dx|,2)+pow(|dy|,2)) using r=|dy|sqrt(1+pow(|dx|/|dy|,2)) if
2344  // |dy|>|dx| [see <A
2345  // HREF=https://en.wikipedia.org/wiki/Hypot">Hypot</A>.].
2346  else if (std::abs(dy) > std::abs(dx))
2347  {
2348  // Since |dy|>|dx| calculate the ratio |dx|/|dy|
2349  std::complex<double> temp = dx / dy;
2350 
2351  // Calculate the value of sin(theta) using:
2352  // sin(theta)=dy/sqrt(pow(|dx|,2)+pow(|dy|,2))
2353  // =(dy/|dy|)/sqrt(1+pow(|dx|/|dy|,2)).
2354  sn = (dy / std::abs(dy)) / sqrt(1.0 + pow(std::abs(temp), 2.0));
2355 
2356  // Calculate the value of cos(theta) using:
2357  // cos(theta)=dx/sqrt(pow(|dy|,2)+pow(|dx|,2))
2358  // =(dx/|dy|)/sqrt(1+pow(|dx|/|dy|,2))
2359  // =(dx/dy)*sin(theta).
2360  cs = temp * sn;
2361  }
2362  // Otherwise, we have |dx|>=|dy| so to, again, avoid overflow or underflow
2363  // calculate the values of cs and sn using the method above
2364  else
2365  {
2366  // Since |dx|>=|dy| calculate the ratio dy/dx
2367  std::complex<double> temp = dy / dx;
2368 
2369  // Calculate the value of cos(theta) using:
2370  // cos(theta)=dx/sqrt(pow(|dx|,2)+pow(|dy|,2))
2371  // =(dx/|dx|)/sqrt(1+pow(|dy|/|dx|,2)).
2372  cs = (dx / std::abs(dx)) / sqrt(1.0 + pow(std::abs(temp), 2.0));
2373 
2374  // Calculate the value of sin(theta) using:
2375  // sin(theta)=dy/sqrt(pow(|dx|,2)+pow(|dy|,2))
2376  // =(dy/|dx|)/sqrt(1+pow(|dy|/|dx|,2))
2377  // =(dy/dx)*cos(theta).
2378  sn = temp * cs;
2379  }
2380 
2381  // Set the tolerance for sin(theta)
2382  double tolerance = 1.0e-15;
2383 
2384  // Make sure sn is real and nonnegative (it should be!)
2385  if ((std::fabs(sn.imag()) > tolerance) || (sn.real() < 0))
2386  {
2387  // Create an output stream
2388  std::ostringstream error_message_stream;
2389 
2390  // Create an error message
2391  error_message_stream << "The value of sin(theta) is not real "
2392  << "and/or nonnegative. Value is: " << sn
2393  << std::endl;
2394 
2395  // Throw an error
2396  throw OomphLibError(error_message_stream.str(),
2397  OOMPH_CURRENT_FUNCTION,
2398  OOMPH_EXCEPTION_LOCATION);
2399  }
2400  } // End of generate_plane_rotation
2401 
2402  /// Helper function: Apply plane rotation. This is done using the
2403  /// update:
2404  /// \f[ \begin{bmatrix} dx \\ dy \end{bmatrix} \leftarrow \begin{bmatrix} \overline{\cos\theta} & \overline{\sin\theta} \\ -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \\ dy \end{bmatrix}. \f]
2405  /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.193.
2406  void apply_plane_rotation(std::complex<double>& dx,
2407  std::complex<double>& dy,
2408  std::complex<double>& cs,
2409  std::complex<double>& sn)
2410  {
2411  // Calculate the value of dx but don't update it yet
2412  std::complex<double> temp = std::conj(cs) * dx + std::conj(sn) * dy;
2413 
2414  // Set the value of dy
2415  dy = -sn * dx + cs * dy;
2416 
2417  // Set the value of dx using the correct values of dx and dy
2418  dx = temp;
2419  } // End of apply_plane_rotation
2420 
2421  /// Number of iterations taken
2422  unsigned Iterations;
2423 
2424  /// Vector of pointers to the real and imaginary part of the system matrix
2426 
2427  /// Boolean flag to indicate if the solve is done in re-solve mode,
2428  /// bypassing setup of matrix and preconditioner
2430 
2431  /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
2432  /// can be deleted.
2434 
2435  /// boolean indicating use of left hand preconditioning (if true)
2436  /// or right hand preconditioning (if false)
2438  };
2439 
2440  /// ////////////////////////////////////////////////////////////////////
2441  /// ////////////////////////////////////////////////////////////////////
2442  /// ////////////////////////////////////////////////////////////////////
2443 
2444  //=============================================================================
2445  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
2446  /// and returns the solution of the linear system.
2447  /// based on the algorithm presented in Templates for the
2448  /// Solution of Linear Systems: Building Blocks for Iterative Methods,
2449  /// Barrett, Berry et al, SIAM, 2006 and the implementation in the IML++
2450  /// library : http://math.nist.gov/iml++/
2451  //=============================================================================
2452  template<typename MATRIX>
2454  DoubleMatrixBase* const& matrix_pt,
2455  const DoubleVector& rhs,
2456  DoubleVector& solution)
2457  {
2458  // Set the number of dof types (real and imaginary for this solver)
2459  unsigned n_dof_types = this->ndof_types();
2460 
2461 #ifdef PARANOID
2462  // This only works for 2 dof types
2463  if (n_dof_types != 2)
2464  {
2465  // Create an output stream
2466  std::stringstream error_message_stream;
2467 
2468  // Create the error message
2469  error_message_stream << "This preconditioner only works for problems "
2470  << "with 2 dof types\nYours has " << n_dof_types;
2471 
2472  // Throw the error message
2473  throw OomphLibError(error_message_stream.str(),
2474  OOMPH_CURRENT_FUNCTION,
2475  OOMPH_EXCEPTION_LOCATION);
2476  }
2477 #endif
2478 
2479  // Get the number of dofs (note, the total number of dofs in the problem
2480  // is 2*n_row but if the constituent vectors and matrices were stored in
2481  // complex objects there would only be (n_row) rows so we use that)
2482  unsigned n_row = Matrices_storage_pt[0]->nrow();
2483 
2484  // Make sure Max_iter isn't greater than n_dof. The user cannot use this
2485  // many iterations when using Krylov subspace methods
2486  if (Max_iter > n_row)
2487  {
2488  // Create an output string stream
2489  std::ostringstream error_message_stream;
2490 
2491  // Create the error message
2492  error_message_stream << "The maximum number of iterations cannot exceed "
2493  << "the number of rows in the problem."
2494  << "\nMaximum number of iterations: " << Max_iter
2495  << "\nNumber of rows: " << n_row << std::endl;
2496 
2497  // Throw the error message
2498  throw OomphLibError(error_message_stream.str(),
2499  OOMPH_CURRENT_FUNCTION,
2500  OOMPH_EXCEPTION_LOCATION);
2501  }
2502 
2503 #ifdef PARANOID
2504  // Loop over the real and imaginary parts
2505  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2506  {
2507  // PARANOID check that if the matrix is distributable then it should not
2508  // be then it should not be distributed
2509  if (dynamic_cast<DistributableLinearAlgebraObject*>(
2510  Matrices_storage_pt[dof_type]) != 0)
2511  {
2512  if (dynamic_cast<DistributableLinearAlgebraObject*>(
2513  Matrices_storage_pt[dof_type])
2514  ->distributed())
2515  {
2516  std::ostringstream error_message_stream;
2517  error_message_stream << "The matrix must not be distributed.";
2518  throw OomphLibError(error_message_stream.str(),
2519  OOMPH_CURRENT_FUNCTION,
2520  OOMPH_EXCEPTION_LOCATION);
2521  }
2522  }
2523  }
2524  // PARANOID check that this rhs distribution is setup
2525  if (!rhs.built())
2526  {
2527  std::ostringstream error_message_stream;
2528  error_message_stream << "The rhs vector distribution must be setup.";
2529  throw OomphLibError(error_message_stream.str(),
2530  OOMPH_CURRENT_FUNCTION,
2531  OOMPH_EXCEPTION_LOCATION);
2532  }
2533  // PARANOID check that the rhs has the right number of global rows
2534  if (rhs.nrow() != 2 * n_row)
2535  {
2536  std::ostringstream error_message_stream;
2537  error_message_stream << "RHS does not have the same dimension as the"
2538  << " linear system";
2539  throw OomphLibError(error_message_stream.str(),
2540  OOMPH_CURRENT_FUNCTION,
2541  OOMPH_EXCEPTION_LOCATION);
2542  }
2543  // PARANOID check that the rhs is not distributed
2544  if (rhs.distribution_pt()->distributed())
2545  {
2546  std::ostringstream error_message_stream;
2547  error_message_stream << "The rhs vector must not be distributed.";
2548  throw OomphLibError(error_message_stream.str(),
2549  OOMPH_CURRENT_FUNCTION,
2550  OOMPH_EXCEPTION_LOCATION);
2551  }
2552  // PARANOID check that if the result is setup it matches the distribution
2553  // of the rhs
2554  if (solution.built())
2555  {
2556  if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
2557  {
2558  std::ostringstream error_message_stream;
2559  error_message_stream << "If the result distribution is setup then it "
2560  << "must be the same as the rhs distribution";
2561  throw OomphLibError(error_message_stream.str(),
2562  OOMPH_CURRENT_FUNCTION,
2563  OOMPH_EXCEPTION_LOCATION);
2564  }
2565  } // if (solution[dof_type].built())
2566 #endif
2567 
2568  // Set up the solution distribution if it's not already distributed
2569  if (!solution.built())
2570  {
2571  // Build the distribution
2573  }
2574  // Otherwise initialise all entries to zero
2575  else
2576  {
2577  // Initialise the entries in the k-th vector in solution to zero
2578  solution.initialise(0.0);
2579  }
2580 
2581  // Create a vector of DoubleVectors (this is a distributed vector so we have
2582  // to create two separate DoubleVector objects to cope with the arithmetic)
2583  Vector<DoubleVector> block_solution(n_dof_types);
2584 
2585  // Create a vector of DoubleVectors (this is a distributed vector so we have
2586  // to create two separate DoubleVector objects to cope with the arithmetic)
2587  Vector<DoubleVector> block_rhs(n_dof_types);
2588 
2589  // Build the distribution of both vectors
2590  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2591  {
2592  // Build the distribution of the k-th constituent vector
2593  block_solution[dof_type].build(this->block_distribution_pt(dof_type),
2594  0.0);
2595 
2596  // Build the distribution of the k-th constituent vector
2597  block_rhs[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
2598  }
2599 
2600  // Grab the solution vector in block form
2601  this->get_block_vectors(solution, block_solution);
2602 
2603  // Grab the RHS vector in block form
2604  this->get_block_vectors(rhs, block_rhs);
2605 
2606  // Start the solver timer
2607  double t_start = TimingHelpers::timer();
2608 
2609  // Storage for the relative residual
2610  double resid;
2611 
2612  // Initialise vectors (since these are not distributed vectors we template
2613  // one vector by the type std::complex<double> instead of using two vectors,
2614  // each templated by the type double
2615 
2616  // Vector, s, used in the minimisation problem: J(y)=min||s-R_m*y||
2617  // [see Saad Y."Iterative methods for sparse linear systems", p.176.]
2618  Vector<std::complex<double>> s(n_row + 1, std::complex<double>(0.0, 0.0));
2619 
2620  // Vector to store the value of cos(theta) when using the Givens rotation
2621  Vector<std::complex<double>> cs(n_row + 1, std::complex<double>(0.0, 0.0));
2622 
2623  // Vector to store the value of sin(theta) when using the Givens rotation
2624  Vector<std::complex<double>> sn(n_row + 1, std::complex<double>(0.0, 0.0));
2625 
2626  // Create a vector of DoubleVectors (this is a distributed vector so we have
2627  // to create two separate DoubleVector objects to cope with the arithmetic)
2628  Vector<DoubleVector> block_w(n_dof_types);
2629 
2630  // Build the distribution of both vectors
2631  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2632  {
2633  // Build the distribution of the k-th constituent vector
2634  block_w[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
2635  }
2636 
2637  // Set up the preconditioner only if we're not re-solving
2638  if (!Resolving)
2639  {
2640  // Only set up the preconditioner before solve if required
2641  if (Setup_preconditioner_before_solve)
2642  {
2643  // Set up the preconditioner from the Jacobian matrix
2644  double t_start_prec = TimingHelpers::timer();
2645 
2646  // Use the setup function in the Preconditioner class
2647  preconditioner_pt()->setup(dynamic_cast<MATRIX*>(matrix_pt));
2648 
2649  // Doc time for setup of preconditioner
2650  double t_end_prec = TimingHelpers::timer();
2651  Preconditioner_setup_time = t_end_prec - t_start_prec;
2652 
2653  // If time documentation is enabled
2654  if (Doc_time)
2655  {
2656  // Output the time taken
2657  oomph_info << "Time for setup of preconditioner [sec]: "
2658  << Preconditioner_setup_time << std::endl;
2659  }
2660  }
2661  }
2662  else
2663  {
2664  // If time documentation is enabled
2665  if (Doc_time)
2666  {
2667  // Notify the user
2668  oomph_info << "Setup of preconditioner is bypassed in resolve mode"
2669  << std::endl;
2670  }
2671  } // if (!Resolving) else
2672 
2673  // Create a vector of DoubleVectors to store the RHS of b-Jx=Mr. We assume
2674  // both x=0 and that a preconditioner is not applied by which we deduce b=r
2675  Vector<DoubleVector> block_r(n_dof_types);
2676 
2677  // Build the distribution of both vectors
2678  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2679  {
2680  // Build the distribution of the k-th constituent vector
2681  block_r[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
2682  }
2683 
2684  // If we're using LHS preconditioning solve b-Jx=Mr for r (assumes x=0)
2685  // so calculate r=M^{-1}b otherwise set r=b (RHS prec.)
2686  if (Preconditioner_LHS)
2687  {
2688  // Create a vector of the same size as rhs
2690 
2691  // Copy the vectors in r to full_r
2692  this->return_block_vectors(block_r, r);
2693 
2694  // Use the preconditioner
2695  preconditioner_pt()->preconditioner_solve(rhs, r);
2696 
2697  // Copy the vector full_r into the vectors in r
2698  this->get_block_vectors(r, block_r);
2699  }
2700  else
2701  {
2702  // Store the value of b (the RHS vector) in r
2703  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2704  {
2705  // Copy the entries of rhs into r
2706  block_r[dof_type] = block_rhs[dof_type];
2707  }
2708  } // if(Preconditioner_LHS)
2709 
2710  // Calculate the norm of the real part of r
2711  double norm_r = block_r[0].norm();
2712 
2713  // Calculate the norm of the imaginary part of r
2714  double norm_c = block_r[1].norm();
2715 
2716  // Compute norm(r)
2717  double normb = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
2718 
2719  // Set the value of beta (the initial residual)
2720  double beta = normb;
2721 
2722  // Compute the initial relative residual. If the entries of the RHS vector
2723  // are all zero then set normb equal to one. This is because we divide the
2724  // value of the norm at later stages by normb and dividing by zero is not
2725  // definied
2726  if (normb == 0.0)
2727  {
2728  // Set the value of normb
2729  normb = 1.0;
2730  }
2731 
2732  // Calculate the ratio between the initial norm and the current norm.
2733  // Since we haven't completed an iteration yet this will simply be one
2734  // unless normb was zero, in which case resid will have value zero
2735  resid = beta / normb;
2736 
2737  // If required, will document convergence history to screen or file (if
2738  // stream open)
2739  if (Doc_convergence_history)
2740  {
2741  // If an output file which is open isn't provided then output to screen
2742  if (!Output_file_stream.is_open())
2743  {
2744  // Output the residual value to the screen
2745  oomph_info << 0 << " " << resid << std::endl;
2746  }
2747  // Otherwise, output to file
2748  else
2749  {
2750  // Output the residual value to file
2751  Output_file_stream << 0 << " " << resid << std::endl;
2752  }
2753  } // if (Doc_convergence_history)
2754 
2755  // If the GMRES algorithm converges immediately
2756  if (resid <= Tolerance)
2757  {
2758  // If time documentation is enabled
2759  if (Doc_time)
2760  {
2761  // Notify the user
2762  oomph_info << "GMRES converged immediately. Normalised residual norm: "
2763  << resid << std::endl;
2764  }
2765 
2766  // Finish running the solver
2767  return;
2768  } // if (resid<=Tolerance)
2769 
2770  // Initialise a vector of orthogonal basis vectors
2771  Vector<Vector<DoubleVector>> block_v;
2772 
2773  // Resize the number of vectors needed
2774  block_v.resize(n_row + 1);
2775 
2776  // Resize each Vector of DoubleVectors to store the real and imaginary
2777  // part of a given vector
2778  for (unsigned dof_type = 0; dof_type < n_row + 1; dof_type++)
2779  {
2780  // Create two DoubleVector objects in each Vector
2781  block_v[dof_type].resize(n_dof_types);
2782  }
2783 
2784  // Initialise the upper hessenberg matrix. Since we are not using
2785  // distributed vectors here, the algebra is best done using entries
2786  // of the type std::complex<double>. NOTE: For implementation purposes
2787  // the upper Hessenberg matrix indices are swapped so the matrix is
2788  // effectively transposed
2789  Vector<Vector<std::complex<double>>> hessenberg(n_row + 1);
2790 
2791  // Build the zeroth basis vector
2792  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2793  {
2794  // Build the k-th part of the zeroth vector. Here k=0 and k=1 correspond
2795  // to the real and imaginary part of the zeroth vector, respectively
2796  block_v[0][dof_type].build(this->block_distribution_pt(dof_type), 0.0);
2797  }
2798 
2799  // Loop over the real and imaginary parts of v
2800  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2801  {
2802  // For fast access
2803  double* v0_pt = block_v[0][dof_type].values_pt();
2804 
2805  // For fast access
2806  const double* block_r_pt = block_r[dof_type].values_pt();
2807 
2808  // Set the zeroth basis vector v[0] to r/beta
2809  for (unsigned i = 0; i < n_row; i++)
2810  {
2811  // Assign the i-th entry of the zeroth basis vector
2812  v0_pt[i] = block_r_pt[i] / beta;
2813  }
2814  } // for (unsigned k=0;k<n_dof_types;k++)
2815 
2816  // Set the first entry in the minimisation problem RHS vector (is meant
2817  // to the vector beta*e_1 initially, where e_1 is the unit vector with
2818  // one in its first entry)
2819  s[0] = beta;
2820 
2821  // Compute the next step of the iterative scheme
2822  for (unsigned j = 0; j < Max_iter; j++)
2823  {
2824  // Resize the next column of the upper hessenberg matrix
2825  hessenberg[j].resize(j + 2, std::complex<double>(0.0, 0.0));
2826 
2827  // Calculate w=M^{-1}(Jv[j]) (LHS prec.) or w=JM^{-1}v (RHS prec.)
2828  {
2829  // Create a temporary vector
2831 
2832  // Create a temporary vector
2834 
2835  // Create a temporary vector
2837 
2838  // Create a temporary vector of DoubleVectors
2839  Vector<DoubleVector> block_temp(2);
2840 
2841  // Create two DoubleVectors
2842  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2843  {
2844  block_temp[dof_type].build(this->block_distribution_pt(dof_type),
2845  0.0);
2846  }
2847 
2848  // If we're using LHS preconditioning
2849  if (Preconditioner_LHS)
2850  {
2851  // Solve Jv[j]=Mw for w. Note, we cannot use inbuilt complex matrix
2852  // algebra here as we're using distributed vectors
2853  complex_matrix_multiplication(
2854  Matrices_storage_pt, block_v[j], block_temp);
2855 
2856  // Copy block_temp into temp
2857  this->return_block_vectors(block_temp, temp);
2858 
2859  // Copy block_w into w
2860  this->return_block_vectors(block_w, w);
2861 
2862  // Apply the preconditioner
2863  this->preconditioner_pt()->preconditioner_solve(temp, w);
2864 
2865  // Copy w into block_w
2866  this->get_block_vectors(w, block_w);
2867  }
2868  // If we're using RHS preconditioning
2869  else
2870  {
2871  // Copy the real and imaginary part of v[j] into one vector, vj
2872  this->return_block_vectors(block_v[j], vj);
2873 
2874  // Use w=JM^{-1}v by saad p270
2875  this->preconditioner_pt()->preconditioner_solve(vj, temp);
2876 
2877  // Copy w into block_w
2878  this->get_block_vectors(temp, block_temp);
2879 
2880  // Solve Jv[j] = Mw for w. Note, we cannot use inbuilt complex matrix
2881  // algebra here as we're using distributed vectors
2882  complex_matrix_multiplication(
2883  Matrices_storage_pt, block_temp, block_w);
2884  }
2885  } // Calculate w=M^{-1}(Jv[j]) (LHS prec.) or w=JM^{-1}v (RHS prec.)
2886 
2887  // For fast access
2888  double* block_w_r_pt = block_w[0].values_pt();
2889 
2890  // For fast access
2891  double* block_w_c_pt = block_w[1].values_pt();
2892 
2893  // Loop over all of the entries on and above the principal subdiagonal of
2894  // the Hessenberg matrix in the j-th column (remembering that
2895  // the indices of the upper Hessenberg matrix are swapped for the purpose
2896  // of implementation)
2897  for (unsigned i = 0; i < j + 1; i++)
2898  {
2899  // For fast access
2900  const double* vi_r_pt = block_v[i][0].values_pt();
2901 
2902  // For fast access
2903  const double* vi_c_pt = block_v[i][1].values_pt();
2904 
2905  // Loop over the entries of v and w
2906  for (unsigned k = 0; k < n_row; k++)
2907  {
2908  // Store the appropriate entry in v as a complex value
2909  std::complex<double> complex_v(vi_r_pt[k], vi_c_pt[k]);
2910 
2911  // Store the appropriate entry in w as a complex value
2912  std::complex<double> complex_w(block_w_r_pt[k], block_w_c_pt[k]);
2913 
2914  // Update the value of H(i,j) noting we're computing a complex
2915  // inner product here (the ordering is very important here!)
2916  hessenberg[j][i] += std::conj(complex_v) * complex_w;
2917  }
2918 
2919  // Orthonormalise w against all previous orthogonal vectors, v_i by
2920  // looping over its entries and updating them
2921  for (unsigned k = 0; k < n_row; k++)
2922  {
2923  // Update the real part of the k-th entry of w
2924  block_w_r_pt[k] -= (hessenberg[j][i].real() * vi_r_pt[k] -
2925  hessenberg[j][i].imag() * vi_c_pt[k]);
2926 
2927  // Update the imaginary part of the k-th entry of w
2928  block_w_c_pt[k] -= (hessenberg[j][i].real() * vi_c_pt[k] +
2929  hessenberg[j][i].imag() * vi_r_pt[k]);
2930  }
2931  } // for (unsigned i=0;i<j+1;i++)
2932 
2933  // Calculate the norm of the real part of w
2934  norm_r = block_w[0].norm();
2935 
2936  // Calculate the norm of the imaginary part of w
2937  norm_c = block_w[1].norm();
2938 
2939  // Calculate the norm of the vector w using norm_r and norm_c and assign
2940  // its value to the appropriate entry in the Hessenberg matrix
2941  hessenberg[j][j + 1] = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
2942 
2943  // Build the (j+1)-th basis vector
2944  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
2945  {
2946  // Build the k-th part of the zeroth vector. Here k=0 and k=1 correspond
2947  // to the real and imaginary part of the zeroth vector, respectively
2948  block_v[j + 1][dof_type].build(this->block_distribution_pt(dof_type),
2949  0.0);
2950  }
2951 
2952  // Check if the value of hessenberg[j][j+1] is zero. If it
2953  // isn't then we update the next entries in v
2954  if (hessenberg[j][j + 1] != 0.0)
2955  {
2956  // For fast access
2957  double* v_r_pt = block_v[j + 1][0].values_pt();
2958 
2959  // For fast access
2960  double* v_c_pt = block_v[j + 1][1].values_pt();
2961 
2962  // For fast access
2963  const double* block_w_r_pt = block_w[0].values_pt();
2964 
2965  // For fast access
2966  const double* block_w_c_pt = block_w[1].values_pt();
2967 
2968  // Notice, the value of H(j,j+1), as calculated above, is clearly a real
2969  // number. As such, calculating the division
2970  // v_{j+1}=w_{j}/h_{j+1,j},
2971  // here is simple, i.e. we don't need to worry about cross terms in the
2972  // algebra. To avoid computing h_{j+1,j} several times we precompute it
2973  double h_subdiag_val = hessenberg[j][j + 1].real();
2974 
2975  // Loop over the entries of the new orthogonal vector and set its values
2976  for (unsigned k = 0; k < n_row; k++)
2977  {
2978  // The i-th entry of the real component is given by
2979  v_r_pt[k] = block_w_r_pt[k] / h_subdiag_val;
2980 
2981  // Similarly, the i-th entry of the imaginary component is given by
2982  v_c_pt[k] = block_w_c_pt[k] / h_subdiag_val;
2983  }
2984  }
2985  // Otherwise, we have to jump to the next part of the algorithm; if
2986  // the value of hessenberg[j][j+1] is zero then the norm of the latest
2987  // orthogonal vector is zero. This is only possible if the entries
2988  // in w are all zero. As a result, the Krylov space of A and r_0 has
2989  // been spanned by the previously calculated orthogonal vectors
2990  else
2991  {
2992  // Book says "Set m=j and jump to step 11" (p.172)...
2993  // Do something here!
2994  oomph_info << "Subdiagonal Hessenberg entry is zero. "
2995  << "Do something here..." << std::endl;
2996  } // if (hessenberg[j][j+1]!=0.0)
2997 
2998  // Loop over the entries in the Hessenberg matrix and calculate the
2999  // entries of the Givens rotation matrices
3000  for (unsigned k = 0; k < j; k++)
3001  {
3002  // Apply the plane rotation to all of the previous entries in the
3003  // (j)-th column (remembering the indexing is reversed)
3004  apply_plane_rotation(
3005  hessenberg[j][k], hessenberg[j][k + 1], cs[k], sn[k]);
3006  }
3007 
3008  // Now calculate the entries of the latest Givens rotation matrix
3009  generate_plane_rotation(
3010  hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
3011 
3012  // Apply the plane rotation using the newly calculated entries
3013  apply_plane_rotation(
3014  hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
3015 
3016  // Apply a plane rotation to the corresponding entry in the vector
3017  // s used in the minimisation problem, J(y)=min||s-R_m*y||
3018  apply_plane_rotation(s[j], s[j + 1], cs[j], sn[j]);
3019 
3020  // Compute current residual using equation (6.42) in Saad Y, "Iterative
3021  // methods for sparse linear systems", p.177.]. Note, since s has complex
3022  // entries we have to use std::abs instead of std::fabs
3023  beta = std::abs(s[j + 1]);
3024 
3025  // Compute the relative residual
3026  resid = beta / normb;
3027 
3028  // If required will document convergence history to screen or file (if
3029  // stream open)
3030  if (Doc_convergence_history)
3031  {
3032  // If an output file which is open isn't provided then output to screen
3033  if (!Output_file_stream.is_open())
3034  {
3035  // Output the residual value to the screen
3036  oomph_info << j + 1 << " " << resid << std::endl;
3037  }
3038  // Otherwise, output to file
3039  else
3040  {
3041  // Output the residual value to file
3042  Output_file_stream << j + 1 << " " << resid << std::endl;
3043  }
3044  } // if (Doc_convergence_history)
3045 
3046  // If the required tolerance has been met
3047  if (resid < Tolerance)
3048  {
3049  // Store the number of iterations taken
3050  Iterations = j + 1;
3051 
3052  // Update the result vector using the result, x=x_0+V_m*y (where V_m
3053  // is given by v here)
3054  update(j, hessenberg, s, block_v, block_solution);
3055 
3056  // Copy the vectors in block_solution to solution
3057  this->return_block_vectors(block_solution, solution);
3058 
3059  // If time documentation was enabled
3060  if (Doc_time)
3061  {
3062  // Output the current normalised residual norm
3063  oomph_info << "\nGMRES converged (1). Normalised residual norm: "
3064  << resid << std::endl;
3065 
3066  // Output the number of iterations it took for convergence
3067  oomph_info << "Number of iterations to convergence: " << j + 1 << "\n"
3068  << std::endl;
3069  }
3070 
3071  // Stop the timer
3072  double t_end = TimingHelpers::timer();
3073 
3074  // Calculate the time taken to calculate the solution
3075  Solution_time = t_end - t_start;
3076 
3077  // If time documentation was enabled
3078  if (Doc_time)
3079  {
3080  // Output the time taken to solve the problem using GMRES
3081  oomph_info << "Time for solve with GMRES [sec]: " << Solution_time
3082  << std::endl;
3083  }
3084 
3085  // As we've met the tolerance for the solver and everything that should
3086  // be documented, has been, finish using the solver
3087  return;
3088  } // if (resid<Tolerance)
3089  } // for (unsigned j=0;j<Max_iter;j++)
3090 
3091  // Store the number of iterations taken
3092  Iterations = Max_iter;
3093 
3094  // Only update if we actually did something
3095  if (Max_iter > 0)
3096  {
3097  // Update the result vector using the result, x=x_0+V_m*y (where V_m
3098  // is given by v here)
3099  update(Max_iter - 1, hessenberg, s, block_v, block_solution);
3100 
3101  // Copy the vectors in block_solution to solution
3102  this->return_block_vectors(block_solution, solution);
3103  }
3104 
3105  // Solve Mr=(b-Jx) for r
3106  {
3107  // Create a temporary vector of DoubleVectors
3108  Vector<DoubleVector> block_temp(2);
3109 
3110  // Create two DoubleVectors
3111  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3112  {
3113  // Build the distribution of the (dof_type)-th vector
3114  block_temp[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
3115  }
3116 
3117  // Calculate the value of Jx
3118  complex_matrix_multiplication(
3119  Matrices_storage_pt, block_solution, block_temp);
3120 
3121  // Get the values pointer of the vector (real)
3122  double* block_temp_r_pt = block_temp[0].values_pt();
3123 
3124  // Get the values pointer of the vector (imaginary)
3125  double* block_temp_c_pt = block_temp[1].values_pt();
3126 
3127  // Get the values pointer of the RHS vector (real)
3128  const double* block_rhs_r_pt = block_rhs[0].values_pt();
3129 
3130  // Get the values pointer of the RHS vector (imaginary)
3131  const double* block_rhs_c_pt = block_rhs[1].values_pt();
3132 
3133  // Loop over the dofs
3134  for (unsigned i = 0; i < n_row; i++)
3135  {
3136  // Calculate b-Jx (real)
3137  block_temp_r_pt[i] = block_rhs_r_pt[i] - block_temp_r_pt[i];
3138 
3139  // Calculate b-Jx (imaginary)
3140  block_temp_c_pt[i] = block_rhs_c_pt[i] - block_temp_c_pt[i];
3141  }
3142 
3143  // If we're using LHS preconditioning
3144  if (Preconditioner_LHS)
3145  {
3146  // Create a temporary DoubleVectors
3148 
3149  // Create a vector of the same size as rhs
3151 
3152  // Copy the vectors in r to full_r
3153  this->return_block_vectors(block_temp, temp);
3154 
3155  // Copy the vectors in r to full_r
3156  this->return_block_vectors(block_r, r);
3157 
3158  // Apply the preconditioner
3159  preconditioner_pt()->preconditioner_solve(temp, r);
3160  }
3161  }
3162 
3163  // Compute the current residual
3164  beta = 0.0;
3165 
3166  // Get access to the values pointer (real)
3167  norm_r = block_r[0].norm();
3168 
3169  // Get access to the values pointer (imaginary)
3170  norm_c = block_r[1].norm();
3171 
3172  // Calculate the full norm
3173  beta = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
3174 
3175  // Calculate the relative residual
3176  resid = beta / normb;
3177 
3178  // If the relative residual lies within tolerance
3179  if (resid < Tolerance)
3180  {
3181  // If time documentation is enabled
3182  if (Doc_time)
3183  {
3184  // Notify the user
3185  oomph_info << "\nGMRES converged (2). Normalised residual norm: "
3186  << resid
3187  << "\nNumber of iterations to convergence: " << Iterations
3188  << "\n"
3189  << std::endl;
3190  }
3191 
3192  // End the timer
3193  double t_end = TimingHelpers::timer();
3194 
3195  // Calculate the time taken for the solver
3196  Solution_time = t_end - t_start;
3197 
3198  // If time documentation is enabled
3199  if (Doc_time)
3200  {
3201  oomph_info << "Time for solve with GMRES [sec]: " << Solution_time
3202  << std::endl;
3203  }
3204  return;
3205  }
3206 
3207  // Otherwise GMRES failed convergence
3208  oomph_info << "\nGMRES did not converge to required tolerance! "
3209  << "\nReturning with normalised residual norm: " << resid
3210  << "\nafter " << Max_iter << " iterations.\n"
3211  << std::endl;
3212 
3213  // Throw an error if requested
3214  if (Throw_error_after_max_iter)
3215  {
3216  std::string err = "Solver failed to converge and you requested an error";
3217  err += " on convergence failures.";
3218  throw OomphLibError(
3219  err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
3220  }
3221 
3222  // Finish using the solver
3223  return;
3224  } // End of solve_helper
3225 
3226 
3227  /// ////////////////////////////////////////////////////////////////////
3228  /// ////////////////////////////////////////////////////////////////////
3229  /// ////////////////////////////////////////////////////////////////////
3230 
3231 
3232  //======================================================================
3233  /// The FGMRES method, i.e. the flexible variant of the GMRES
3234  /// method which allows for nonconstant preconditioners [see Saad Y,
3235  /// "Iterative methods for sparse linear systems", p.287]. Note, FGMRES
3236  /// can only cater to right preconditioning; if the user tries to switch
3237  /// to left preconditioning they will be notified of this
3238  //======================================================================
3239  template<typename MATRIX>
3240  class HelmholtzFGMRESMG : public virtual HelmholtzGMRESMG<MATRIX>
3241  {
3242  public:
3243  /// Constructor (empty)
3245  {
3246  // Can only use RHS preconditioning
3247  this->Preconditioner_LHS = false;
3248  };
3249 
3250  /// Destructor (cleanup storage)
3252  {
3253  // Call the clean up function in the base class GMRES
3254  this->clean_up_memory();
3255  }
3256 
3257  /// Broken copy constructor
3259 
3260  /// Broken assignment operator
3261  void operator=(const HelmholtzFGMRESMG&) = delete;
3262 
3263  /// Overloaded function to let the user know that left
3264  /// preconditioning is not possible with FGMRES, only right preconditioning
3266  {
3267  // Create an output stream
3268  std::ostringstream error_message_stream;
3269 
3270  // Create an error message
3271  error_message_stream << "FGMRES cannot use left preconditioning. It is "
3272  << "only capable of using right preconditioning."
3273  << std::endl;
3274 
3275  // Throw the error message
3276  throw OomphLibError(error_message_stream.str(),
3277  OOMPH_CURRENT_FUNCTION,
3278  OOMPH_EXCEPTION_LOCATION);
3279  } // End of set_preconditioner_LHS
3280 
3281  /// Solver: Takes pointer to problem and returns the results vector
3282  /// which contains the solution of the linear system defined by
3283  /// the problem's fully assembled Jacobian and residual vector.
3284  void solve(Problem* const& problem_pt, DoubleVector& result)
3285  {
3286 #ifdef OOMPH_HAS_MPI
3287  // Make sure that this is running in serial. Can't guarantee it'll
3288  // work when the problem is distributed over several processors
3289  if (MPI_Helpers::communicator_pt()->nproc() > 1)
3290  {
3291  // Throw a warning
3292  OomphLibWarning("Can't guarantee the MG solver will work in parallel!",
3293  OOMPH_CURRENT_FUNCTION,
3294  OOMPH_EXCEPTION_LOCATION);
3295  }
3296 #endif
3297 
3298  // Find # of degrees of freedom (variables)
3299  unsigned n_dof = problem_pt->ndof();
3300 
3301  // Initialise timer
3302  double t_start = TimingHelpers::timer();
3303 
3304  // We're not re-solving
3305  this->Resolving = false;
3306 
3307  // Get rid of any previously stored data
3308  this->clean_up_memory();
3309 
3310  // Grab the communicator from the MGProblem object and assign it
3311  this->set_comm_pt(problem_pt->communicator_pt());
3312 
3313  // Setup the distribution
3315  problem_pt->communicator_pt(), n_dof, false);
3316 
3317  // Build the internal distribution in this way because both the
3318  // IterativeLinearSolver and BlockPreconditioner class have base-
3319  // class subobjects of type oomph::DistributableLinearAlgebraObject
3321 
3322  // Get Jacobian matrix in format specified by template parameter
3323  // and nonlinear residual vector
3324  MATRIX* matrix_pt = new MATRIX;
3325  DoubleVector f;
3326  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
3327  {
3328  if (dynamic_cast<CRDoubleMatrix*>(matrix_pt) != 0)
3329  {
3330  dynamic_cast<CRDoubleMatrix*>(matrix_pt)->build(
3333  }
3334  }
3335 
3336  // Get the Jacobian and residuals vector
3337  problem_pt->get_jacobian(f, *matrix_pt);
3338 
3339  // We've made the matrix, we can delete it...
3340  this->Matrix_can_be_deleted = true;
3341 
3342  // Replace the current matrix used in Preconditioner by the new matrix
3343  this->set_matrix_pt(matrix_pt);
3344 
3345  // The preconditioner works with one mesh; set it! Since we only use
3346  // the block preconditioner on the finest level, we use the mesh from
3347  // that level
3348  this->set_nmesh(1);
3349 
3350  // Elements in actual pml layer are trivially wrapped versions of
3351  // their bulk counterparts. Technically they are different elements
3352  // so we have to allow different element types
3353  bool allow_different_element_types_in_mesh = true;
3354  this->set_mesh(
3355  0, problem_pt->mesh_pt(), allow_different_element_types_in_mesh);
3356 
3357  // Set up the generic block look up scheme
3358  this->block_setup();
3359 
3360  // Extract the number of blocks.
3361  unsigned nblock_types = this->nblock_types();
3362 
3363 #ifdef PARANOID
3364  // PARANOID check - there must only be two block types
3365  if (nblock_types != 2)
3366  {
3367  // Create the error message
3368  std::stringstream tmp;
3369  tmp << "There are supposed to be two block types.\nYours has "
3370  << nblock_types << std::endl;
3371 
3372  // Throw an error
3373  throw OomphLibError(
3374  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3375  }
3376 #endif
3377 
3378  // Resize the storage for the system matrices
3379  this->Matrices_storage_pt.resize(2, 0);
3380 
3381  // Loop over the rows of the block matrix
3382  for (unsigned i = 0; i < nblock_types; i++)
3383  {
3384  // Fix the column index
3385  unsigned j = 0;
3386 
3387  // Create new CRDoubleMatrix objects
3388  this->Matrices_storage_pt[i] = new CRDoubleMatrix;
3389 
3390  // Extract the required blocks, i.e. the first column
3391  this->get_block(i, j, *(this->Matrices_storage_pt[i]));
3392  }
3393 
3394  // Doc time for setup
3395  double t_end = TimingHelpers::timer();
3396  this->Jacobian_setup_time = t_end - t_start;
3397 
3398  if (this->Doc_time)
3399  {
3400  oomph_info << "\nTime for setup of block Jacobian [sec]: "
3401  << this->Jacobian_setup_time << std::endl;
3402  }
3403 
3404  // Call linear algebra-style solver
3405  // If the result distribution is wrong, then redistribute
3406  // before the solve and return to original distribution
3407  // afterwards
3408  if ((!(*result.distribution_pt() ==
3410  (result.built()))
3411  {
3412  // Make a distribution object
3413  LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
3414 
3415  // Build the result vector distribution
3417 
3418  // Solve the problem
3419  solve_helper(matrix_pt, f, result);
3420 
3421  // Redistribute the vector
3422  result.redistribute(&temp_global_dist);
3423  }
3424  // Otherwise just solve
3425  else
3426  {
3427  // Solve
3428  solve_helper(matrix_pt, f, result);
3429  }
3430 
3431  // Kill matrix unless it's still required for resolve
3432  if (!(this->Enable_resolve))
3433  {
3434  // Clean up anything left in memory
3435  this->clean_up_memory();
3436  }
3437  } // End of solve
3438 
3439  private:
3440  /// General interface to solve function
3442  const DoubleVector& rhs,
3443  DoubleVector& solution);
3444 
3445  /// Helper function to update the result vector
3446  void update(const unsigned& k,
3447  const Vector<Vector<std::complex<double>>>& hessenberg,
3448  const Vector<std::complex<double>>& s,
3449  const Vector<Vector<DoubleVector>>& z_m,
3451  {
3452  // Make a local copy of s
3454 
3455  //-----------------------------------------------------------------
3456  // The Hessenberg matrix should be an upper triangular matrix at
3457  // this point (although from its storage it would appear to be a
3458  // lower triangular matrix since the indexing has been reversed)
3459  // so finding the minimiser of J(y)=min||s-R_m*y|| where R_m is
3460  // the matrix R in the QR factorisation of the Hessenberg matrix.
3461  // Therefore, to obtain y we simply need to use a backwards
3462  // substitution. Note: The implementation here may appear to be
3463  // somewhat confusing as the indexing in the Hessenberg matrix is
3464  // reversed. This implementation of a backwards substitution does
3465  // not run along the columns of the triangular matrix but rather
3466  // up the rows.
3467  //-----------------------------------------------------------------
3468 
3469  // The outer loop is a loop over the columns of the Hessenberg matrix
3470  // since the indexing is reversed
3471  for (int i = int(k); i >= 0; i--)
3472  {
3473  // Divide the i-th entry of y by the i-th diagonal entry of H
3474  y[i] /= hessenberg[i][i];
3475 
3476  // The inner loop is a loop over the rows of the Hessenberg matrix
3477  for (int j = i - 1; j >= 0; j--)
3478  {
3479  // Update the j-th entry of y
3480  y[j] -= hessenberg[i][j] * y[i];
3481  }
3482  } // for (int i=int(k);i>=0;i--)
3483 
3484  // Calculate the number of entries in x (simply use the real part as
3485  // both the real and imaginary part should have the same length)
3486  unsigned n_row = x[0].nrow();
3487 
3488  // Build a temporary vector with entries initialised to 0.0
3489  Vector<DoubleVector> block_update(2);
3490 
3491  // Build the distributions
3492  for (unsigned dof_type = 0; dof_type < 2; dof_type++)
3493  {
3494  // Build the (dof_type)-th vector of block_update
3495  block_update[dof_type].build(x[0].distribution_pt(), 0.0);
3496  }
3497 
3498  // Get access to the underlying values
3499  double* block_update_r_pt = block_update[0].values_pt();
3500 
3501  // Get access to the underlying values
3502  double* block_update_c_pt = block_update[1].values_pt();
3503 
3504  // Calculate x=Vy
3505  for (unsigned j = 0; j <= k; j++)
3506  {
3507  // Get access to j-th column of Z_m
3508  const double* z_mj_r_pt = z_m[j][0].values_pt();
3509 
3510  // Get access to j-th column of Z_m
3511  const double* z_mj_c_pt = z_m[j][1].values_pt();
3512 
3513  // Loop over the entries in x and update them
3514  for (unsigned i = 0; i < n_row; i++)
3515  {
3516  // Update the real part of the i-th entry in x
3517  block_update_r_pt[i] +=
3518  (z_mj_r_pt[i] * y[j].real()) - (z_mj_c_pt[i] * y[j].imag());
3519 
3520  // Update the imaginary part of the i-th entry in x
3521  block_update_c_pt[i] +=
3522  (z_mj_c_pt[i] * y[j].real()) + (z_mj_r_pt[i] * y[j].imag());
3523  }
3524  } // for (unsigned j=0;j<=k;j++)
3525 
3526  // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
3527  // sparse linear systems", p.284]
3528  for (unsigned dof_type = 0; dof_type < 2; dof_type++)
3529  {
3530  // Update the solution
3531  x[dof_type] += block_update[dof_type];
3532  }
3533  } // End of update
3534  };
3535 
3536 
3537  /// ///////////////////////////////////////////////////////////////////
3538  /// ///////////////////////////////////////////////////////////////////
3539  /// ///////////////////////////////////////////////////////////////////
3540 
3541 
3542  //=============================================================================
3543  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
3544  /// and returns the solution of the linear system.
3545  /// based on the algorithm presented in Templates for the
3546  /// Solution of Linear Systems: Building Blocks for Iterative Methods,
3547  /// Barrett, Berry et al, SIAM, 2006 and the implementation in the IML++
3548  /// library : http://math.nist.gov/iml++/
3549  //=============================================================================
3550  template<typename MATRIX>
3552  DoubleMatrixBase* const& matrix_pt,
3553  const DoubleVector& rhs,
3554  DoubleVector& solution)
3555  {
3556  // Set the number of dof types (real and imaginary for this solver)
3557  unsigned n_dof_types = this->ndof_types();
3558 
3559 #ifdef PARANOID
3560  // This only works for 2 dof types
3561  if (n_dof_types != 2)
3562  {
3563  // Create an output stream
3564  std::stringstream error_message_stream;
3565 
3566  // Create the error message
3567  error_message_stream << "This preconditioner only works for problems "
3568  << "with 2 dof types\nYours has " << n_dof_types;
3569 
3570  // Throw the error message
3571  throw OomphLibError(error_message_stream.str(),
3572  OOMPH_CURRENT_FUNCTION,
3573  OOMPH_EXCEPTION_LOCATION);
3574  }
3575 #endif
3576 
3577  // Get the number of dofs (note, the total number of dofs in the problem
3578  // is 2*n_row but if the constituent vectors and matrices were stored in
3579  // complex objects there would only be (n_row) rows so we use that)
3580  unsigned n_row = this->Matrices_storage_pt[0]->nrow();
3581 
3582  // Make sure Max_iter isn't greater than n_dof. The user cannot use this
3583  // many iterations when using Krylov subspace methods
3584  if (this->Max_iter > n_row)
3585  {
3586  // Create an output string stream
3587  std::ostringstream error_message_stream;
3588 
3589  // Create the error message
3590  error_message_stream << "The maximum number of iterations cannot exceed "
3591  << "the number of rows in the problem."
3592  << "\nMaximum number of iterations: "
3593  << this->Max_iter << "\nNumber of rows: " << n_row
3594  << std::endl;
3595 
3596  // Throw the error message
3597  throw OomphLibError(error_message_stream.str(),
3598  OOMPH_CURRENT_FUNCTION,
3599  OOMPH_EXCEPTION_LOCATION);
3600  }
3601 
3602 #ifdef PARANOID
3603  // Loop over the real and imaginary parts
3604  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3605  {
3606  // PARANOID check that if the matrix is distributable then it should not
3607  // be then it should not be distributed
3608  if (dynamic_cast<DistributableLinearAlgebraObject*>(
3609  this->Matrices_storage_pt[dof_type]) != 0)
3610  {
3611  if (dynamic_cast<DistributableLinearAlgebraObject*>(
3612  this->Matrices_storage_pt[dof_type])
3613  ->distributed())
3614  {
3615  std::ostringstream error_message_stream;
3616  error_message_stream << "The matrix must not be distributed.";
3617  throw OomphLibError(error_message_stream.str(),
3618  OOMPH_CURRENT_FUNCTION,
3619  OOMPH_EXCEPTION_LOCATION);
3620  }
3621  }
3622  }
3623  // PARANOID check that this rhs distribution is setup
3624  if (!rhs.built())
3625  {
3626  std::ostringstream error_message_stream;
3627  error_message_stream << "The rhs vector distribution must be setup.";
3628  throw OomphLibError(error_message_stream.str(),
3629  OOMPH_CURRENT_FUNCTION,
3630  OOMPH_EXCEPTION_LOCATION);
3631  }
3632  // PARANOID check that the rhs has the right number of global rows
3633  if (rhs.nrow() != 2 * n_row)
3634  {
3635  std::ostringstream error_message_stream;
3636  error_message_stream << "RHS does not have the same dimension as the"
3637  << " linear system";
3638  throw OomphLibError(error_message_stream.str(),
3639  OOMPH_CURRENT_FUNCTION,
3640  OOMPH_EXCEPTION_LOCATION);
3641  }
3642  // PARANOID check that the rhs is not distributed
3643  if (rhs.distribution_pt()->distributed())
3644  {
3645  std::ostringstream error_message_stream;
3646  error_message_stream << "The rhs vector must not be distributed.";
3647  throw OomphLibError(error_message_stream.str(),
3648  OOMPH_CURRENT_FUNCTION,
3649  OOMPH_EXCEPTION_LOCATION);
3650  }
3651  // PARANOID check that if the result is setup it matches the distribution
3652  // of the rhs
3653  if (solution.built())
3654  {
3655  if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
3656  {
3657  std::ostringstream error_message_stream;
3658  error_message_stream << "If the result distribution is setup then it "
3659  << "must be the same as the rhs distribution";
3660  throw OomphLibError(error_message_stream.str(),
3661  OOMPH_CURRENT_FUNCTION,
3662  OOMPH_EXCEPTION_LOCATION);
3663  }
3664  } // if (solution[dof_type].built())
3665 #endif
3666 
3667  // Set up the solution distribution if it's not already distributed
3668  if (!solution.built())
3669  {
3670  // Build the distribution
3672  }
3673  // Otherwise initialise all entries to zero
3674  else
3675  {
3676  // Initialise the entries in the k-th vector in solution to zero
3677  solution.initialise(0.0);
3678  }
3679 
3680  // Create a vector of DoubleVectors (this is a distributed vector so we have
3681  // to create two separate DoubleVector objects to cope with the arithmetic)
3682  Vector<DoubleVector> block_solution(n_dof_types);
3683 
3684  // Create a vector of DoubleVectors (this is a distributed vector so we have
3685  // to create two separate DoubleVector objects to cope with the arithmetic)
3686  Vector<DoubleVector> block_rhs(n_dof_types);
3687 
3688  // Build the distribution of both vectors
3689  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3690  {
3691  // Build the distribution of the k-th constituent vector
3692  block_solution[dof_type].build(this->block_distribution_pt(dof_type),
3693  0.0);
3694 
3695  // Build the distribution of the k-th constituent vector
3696  block_rhs[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
3697  }
3698 
3699  // Grab the solution vector in block form
3700  this->get_block_vectors(solution, block_solution);
3701 
3702  // Grab the RHS vector in block form
3703  this->get_block_vectors(rhs, block_rhs);
3704 
3705  // Start the solver timer
3706  double t_start = TimingHelpers::timer();
3707 
3708  // Storage for the relative residual
3709  double resid;
3710 
3711  // Initialise vectors (since these are not distributed vectors we template
3712  // one vector by the type std::complex<double> instead of using two vectors,
3713  // each templated by the type double
3714 
3715  // Vector, s, used in the minimisation problem: J(y)=min||s-R_m*y||
3716  // [see Saad Y."Iterative methods for sparse linear systems", p.176.]
3717  Vector<std::complex<double>> s(n_row + 1, std::complex<double>(0.0, 0.0));
3718 
3719  // Vector to store the value of cos(theta) when using the Givens rotation
3720  Vector<std::complex<double>> cs(n_row + 1, std::complex<double>(0.0, 0.0));
3721 
3722  // Vector to store the value of sin(theta) when using the Givens rotation
3723  Vector<std::complex<double>> sn(n_row + 1, std::complex<double>(0.0, 0.0));
3724 
3725  // Create a vector of DoubleVectors (this is a distributed vector so we have
3726  // to create two separate DoubleVector objects to cope with the arithmetic)
3727  Vector<DoubleVector> block_w(n_dof_types);
3728 
3729  // Build the distribution of both vectors
3730  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3731  {
3732  // Build the distribution of the k-th constituent vector
3733  block_w[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
3734  }
3735 
3736  // Set up the preconditioner only if we're not re-solving
3737  if (!(this->Resolving))
3738  {
3739  // Only set up the preconditioner before solve if required
3740  if (this->Setup_preconditioner_before_solve)
3741  {
3742  // Set up the preconditioner from the Jacobian matrix
3743  double t_start_prec = TimingHelpers::timer();
3744 
3745  // Use the setup function in the Preconditioner class
3746  this->preconditioner_pt()->setup(dynamic_cast<MATRIX*>(matrix_pt));
3747 
3748  // Doc time for setup of preconditioner
3749  double t_end_prec = TimingHelpers::timer();
3750  this->Preconditioner_setup_time = t_end_prec - t_start_prec;
3751 
3752  // If time documentation is enabled
3753  if (this->Doc_time)
3754  {
3755  // Output the time taken
3756  oomph_info << "Time for setup of preconditioner [sec]: "
3757  << this->Preconditioner_setup_time << std::endl;
3758  }
3759  }
3760  }
3761  else
3762  {
3763  // If time documentation is enabled
3764  if (this->Doc_time)
3765  {
3766  // Notify the user
3767  oomph_info << "Setup of preconditioner is bypassed in resolve mode"
3768  << std::endl;
3769  }
3770  } // if (!Resolving) else
3771 
3772  // Create a vector of DoubleVectors to store the RHS of b-Jx=Mr. We assume
3773  // both x=0 and that a preconditioner is not applied by which we deduce b=r
3774  Vector<DoubleVector> block_r(n_dof_types);
3775 
3776  // Build the distribution of both vectors
3777  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3778  {
3779  // Build the distribution of the k-th constituent vector
3780  block_r[dof_type].build(this->block_distribution_pt(dof_type), 0.0);
3781  }
3782 
3783  // Store the value of b (the RHS vector) in r
3784  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3785  {
3786  // Copy the entries of rhs into r
3787  block_r[dof_type] = block_rhs[dof_type];
3788  }
3789 
3790  // Calculate the norm of the real part of r
3791  double norm_r = block_r[0].norm();
3792 
3793  // Calculate the norm of the imaginary part of r
3794  double norm_c = block_r[1].norm();
3795 
3796  // Compute norm(r)
3797  double normb = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
3798 
3799  // Set the value of beta (the initial residual)
3800  double beta = normb;
3801 
3802  // Compute the initial relative residual. If the entries of the RHS vector
3803  // are all zero then set normb equal to one. This is because we divide the
3804  // value of the norm at later stages by normb and dividing by zero is not
3805  // definied
3806  if (normb == 0.0)
3807  {
3808  // Set the value of normb
3809  normb = 1.0;
3810  }
3811 
3812  // Calculate the ratio between the initial norm and the current norm.
3813  // Since we haven't completed an iteration yet this will simply be one
3814  // unless normb was zero, in which case resid will have value zero
3815  resid = beta / normb;
3816 
3817  // If required, will document convergence history to screen or file (if
3818  // stream open)
3819  if (this->Doc_convergence_history)
3820  {
3821  // If an output file which is open isn't provided then output to screen
3822  if (!(this->Output_file_stream.is_open()))
3823  {
3824  // Output the residual value to the screen
3825  oomph_info << 0 << " " << resid << std::endl;
3826  }
3827  // Otherwise, output to file
3828  else
3829  {
3830  // Output the residual value to file
3831  this->Output_file_stream << 0 << " " << resid << std::endl;
3832  }
3833  } // if (Doc_convergence_history)
3834 
3835  // If the GMRES algorithm converges immediately
3836  if (resid <= this->Tolerance)
3837  {
3838  // If time documentation is enabled
3839  if (this->Doc_time)
3840  {
3841  // Notify the user
3842  oomph_info << "FGMRES converged immediately. Normalised residual norm: "
3843  << resid << std::endl;
3844  }
3845 
3846  // Finish running the solver
3847  return;
3848  } // if (resid<=Tolerance)
3849 
3850  // Initialise a vector of orthogonal basis vectors
3851  Vector<Vector<DoubleVector>> block_v;
3852 
3853  // Resize the number of vectors needed
3854  block_v.resize(n_row + 1);
3855 
3856  // Create a vector of DoubleVectors (stores the preconditioned vectors)
3857  Vector<Vector<DoubleVector>> block_z;
3858 
3859  // Resize the number of vectors needed
3860  block_z.resize(n_row + 1);
3861 
3862  // Resize each Vector of DoubleVectors to store the real and imaginary
3863  // part of a given vector
3864  for (unsigned i = 0; i < n_row + 1; i++)
3865  {
3866  // Create space for two DoubleVector objects in each Vector
3867  block_v[i].resize(n_dof_types);
3868 
3869  // Create space for two DoubleVector objects in each Vector
3870  block_z[i].resize(n_dof_types);
3871  }
3872 
3873  // Initialise the upper hessenberg matrix. Since we are not using
3874  // distributed vectors here, the algebra is best done using entries
3875  // of the type std::complex<double>. NOTE: For implementation purposes
3876  // the upper Hessenberg matrix indices are swapped so the matrix is
3877  // effectively transposed
3878  Vector<Vector<std::complex<double>>> hessenberg(n_row + 1);
3879 
3880  // Build the zeroth basis vector
3881  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3882  {
3883  // Build the k-th part of the zeroth vector. Here k=0 and k=1 correspond
3884  // to the real and imaginary part of the zeroth vector, respectively
3885  block_v[0][dof_type].build(this->block_distribution_pt(dof_type), 0.0);
3886  }
3887 
3888  // Loop over the real and imaginary parts of v
3889  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3890  {
3891  // For fast access
3892  double* block_v0_pt = block_v[0][dof_type].values_pt();
3893 
3894  // For fast access
3895  const double* block_r_pt = block_r[dof_type].values_pt();
3896 
3897  // Set the zeroth basis vector v[0] to r/beta
3898  for (unsigned i = 0; i < n_row; i++)
3899  {
3900  // Assign the i-th entry of the zeroth basis vector
3901  block_v0_pt[i] = block_r_pt[i] / beta;
3902  }
3903  } // for (unsigned k=0;k<n_dof_types;k++)
3904 
3905  // Set the first entry in the minimisation problem RHS vector (is meant
3906  // to the vector beta*e_1 initially, where e_1 is the unit vector with
3907  // one in its first entry)
3908  s[0] = beta;
3909 
3910  // Compute the next step of the iterative scheme
3911  for (unsigned j = 0; j < this->Max_iter; j++)
3912  {
3913  // Resize the next column of the upper hessenberg matrix
3914  hessenberg[j].resize(j + 2, std::complex<double>(0.0, 0.0));
3915 
3916  // Calculate w_j=Jz_j where z_j=M^{-1}v_j (RHS prec.)
3917  {
3918  // Create a temporary vector
3920 
3921  // Create a temporary vector
3923 
3924  // Create a temporary vector
3926 
3927  // Create two DoubleVectors
3928  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
3929  {
3930  // Build the k-th part of the j-th preconditioning result vector
3931  block_z[j][dof_type].build(this->block_distribution_pt(dof_type),
3932  0.0);
3933  }
3934 
3935  // Copy the real and imaginary part of v[j] into one vector, vj
3936  this->return_block_vectors(block_v[j], vj);
3937 
3938  // Calculate z_j=M^{-1}v_j by saad p270
3939  this->preconditioner_pt()->preconditioner_solve(vj, zj);
3940 
3941  // Copy zj into z[j][0] and z[j][1]
3942  this->get_block_vectors(zj, block_z[j]);
3943 
3944  // Calculate w_j=J*(M^{-1}v_j). Note, we cannot use inbuilt complex
3945  // matrix algebra here as we're using distributed vectors
3946  this->complex_matrix_multiplication(
3947  this->Matrices_storage_pt, block_z[j], block_w);
3948  } // Calculate w=JM^{-1}v (RHS prec.)
3949 
3950  // For fast access
3951  double* block_w_r_pt = block_w[0].values_pt();
3952 
3953  // For fast access
3954  double* block_w_c_pt = block_w[1].values_pt();
3955 
3956  // Loop over all of the entries on and above the principal subdiagonal of
3957  // the Hessenberg matrix in the j-th column (remembering that
3958  // the indices of the upper Hessenberg matrix are swapped for the purpose
3959  // of implementation)
3960  for (unsigned i = 0; i < j + 1; i++)
3961  {
3962  // For fast access
3963  const double* block_vi_r_pt = block_v[i][0].values_pt();
3964 
3965  // For fast access
3966  const double* block_vi_c_pt = block_v[i][1].values_pt();
3967 
3968  // Loop over the entries of v and w
3969  for (unsigned k = 0; k < n_row; k++)
3970  {
3971  // Store the appropriate entry in v as a complex value
3972  std::complex<double> complex_v(block_vi_r_pt[k], block_vi_c_pt[k]);
3973 
3974  // Store the appropriate entry in w as a complex value
3975  std::complex<double> complex_w(block_w_r_pt[k], block_w_c_pt[k]);
3976 
3977  // Update the value of H(i,j) noting we're computing a complex
3978  // inner product here (the ordering is very important here!)
3979  hessenberg[j][i] += std::conj(complex_v) * complex_w;
3980  }
3981 
3982  // Orthonormalise w against all previous orthogonal vectors, v_i by
3983  // looping over its entries and updating them
3984  for (unsigned k = 0; k < n_row; k++)
3985  {
3986  // Update the real part of the k-th entry of w
3987  block_w_r_pt[k] -= (hessenberg[j][i].real() * block_vi_r_pt[k] -
3988  hessenberg[j][i].imag() * block_vi_c_pt[k]);
3989 
3990  // Update the imaginary part of the k-th entry of w
3991  block_w_c_pt[k] -= (hessenberg[j][i].real() * block_vi_c_pt[k] +
3992  hessenberg[j][i].imag() * block_vi_r_pt[k]);
3993  }
3994  } // for (unsigned i=0;i<j+1;i++)
3995 
3996  // Calculate the norm of the real part of w
3997  norm_r = block_w[0].norm();
3998 
3999  // Calculate the norm of the imaginary part of w
4000  norm_c = block_w[1].norm();
4001 
4002  // Calculate the norm of the vector w using norm_r and norm_c and assign
4003  // its value to the appropriate entry in the Hessenberg matrix
4004  hessenberg[j][j + 1] = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
4005 
4006  // Build the (j+1)-th basis vector
4007  for (unsigned dof_type = 0; dof_type < n_dof_types; dof_type++)
4008  {
4009  // Build the k-th part of the zeroth vector. Here k=0 and k=1 correspond
4010  // to the real and imaginary part of the zeroth vector, respectively
4011  block_v[j + 1][dof_type].build(this->block_distribution_pt(dof_type),
4012  0.0);
4013  }
4014 
4015  // Check if the value of hessenberg[j][j+1] is zero. If it
4016  // isn't then we update the next entries in v
4017  if (hessenberg[j][j + 1] != 0.0)
4018  {
4019  // For fast access
4020  double* block_v_r_pt = block_v[j + 1][0].values_pt();
4021 
4022  // For fast access
4023  double* block_v_c_pt = block_v[j + 1][1].values_pt();
4024 
4025  // For fast access
4026  const double* block_w_r_pt = block_w[0].values_pt();
4027 
4028  // For fast access
4029  const double* block_w_c_pt = block_w[1].values_pt();
4030 
4031  // Notice, the value of H(j,j+1), as calculated above, is clearly a real
4032  // number. As such, calculating the division
4033  // v_{j+1}=w_{j}/h_{j+1,j},
4034  // here is simple, i.e. we don't need to worry about cross terms in the
4035  // algebra. To avoid computing h_{j+1,j} several times we precompute it
4036  double h_subdiag_val = hessenberg[j][j + 1].real();
4037 
4038  // Loop over the entries of the new orthogonal vector and set its values
4039  for (unsigned k = 0; k < n_row; k++)
4040  {
4041  // The i-th entry of the real component is given by
4042  block_v_r_pt[k] = block_w_r_pt[k] / h_subdiag_val;
4043 
4044  // Similarly, the i-th entry of the imaginary component is given by
4045  block_v_c_pt[k] = block_w_c_pt[k] / h_subdiag_val;
4046  }
4047  }
4048  // Otherwise, we have to jump to the next part of the algorithm; if
4049  // the value of hessenberg[j][j+1] is zero then the norm of the latest
4050  // orthogonal vector is zero. This is only possible if the entries
4051  // in w are all zero. As a result, the Krylov space of A and r_0 has
4052  // been spanned by the previously calculated orthogonal vectors
4053  else
4054  {
4055  // Book says "Set m=j and jump to step 11" (p.172)...
4056  // Do something here!
4057  oomph_info << "Subdiagonal Hessenberg entry is zero. "
4058  << "Do something here..." << std::endl;
4059  } // if (hessenberg[j][j+1]!=0.0)
4060 
4061  // Loop over the entries in the Hessenberg matrix and calculate the
4062  // entries of the Givens rotation matrices
4063  for (unsigned k = 0; k < j; k++)
4064  {
4065  // Apply the plane rotation to all of the previous entries in the
4066  // (j)-th column (remembering the indexing is reversed)
4067  this->apply_plane_rotation(
4068  hessenberg[j][k], hessenberg[j][k + 1], cs[k], sn[k]);
4069  }
4070 
4071  // Now calculate the entries of the latest Givens rotation matrix
4072  this->generate_plane_rotation(
4073  hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
4074 
4075  // Apply the plane rotation using the newly calculated entries
4076  this->apply_plane_rotation(
4077  hessenberg[j][j], hessenberg[j][j + 1], cs[j], sn[j]);
4078 
4079  // Apply a plane rotation to the corresponding entry in the vector
4080  // s used in the minimisation problem, J(y)=min||s-R_m*y||
4081  this->apply_plane_rotation(s[j], s[j + 1], cs[j], sn[j]);
4082 
4083  // Compute current residual using equation (6.42) in Saad Y, "Iterative
4084  // methods for sparse linear systems", p.177.]. Note, since s has complex
4085  // entries we have to use std::abs instead of std::fabs
4086  beta = std::abs(s[j + 1]);
4087 
4088  // Compute the relative residual
4089  resid = beta / normb;
4090 
4091  // If required will document convergence history to screen or file (if
4092  // stream open)
4093  if (this->Doc_convergence_history)
4094  {
4095  // If an output file which is open isn't provided then output to screen
4096  if (!(this->Output_file_stream.is_open()))
4097  {
4098  // Output the residual value to the screen
4099  oomph_info << j + 1 << " " << resid << std::endl;
4100  }
4101  // Otherwise, output to file
4102  else
4103  {
4104  // Output the residual value to file
4105  this->Output_file_stream << j + 1 << " " << resid << std::endl;
4106  }
4107  } // if (Doc_convergence_history)
4108 
4109  // If the required tolerance has been met
4110  if (resid < this->Tolerance)
4111  {
4112  // Store the number of iterations taken
4113  this->Iterations = j + 1;
4114 
4115  // Update the result vector using the result, x=x_0+V_m*y (where V_m
4116  // is given by v here)
4117  update(j, hessenberg, s, block_z, block_solution);
4118 
4119  // Copy the vectors in block_solution to solution
4120  this->return_block_vectors(block_solution, solution);
4121 
4122  // If time documentation was enabled
4123  if (this->Doc_time)
4124  {
4125  // Output the current normalised residual norm
4126  oomph_info << "\nFGMRES converged (1). Normalised residual norm: "
4127  << resid << std::endl;
4128 
4129  // Output the number of iterations it took for convergence
4130  oomph_info << "Number of iterations to convergence: " << j + 1 << "\n"
4131  << std::endl;
4132  }
4133 
4134  // Stop the timer
4135  double t_end = TimingHelpers::timer();
4136 
4137  // Calculate the time taken to calculate the solution
4138  this->Solution_time = t_end - t_start;
4139 
4140  // If time documentation was enabled
4141  if (this->Doc_time)
4142  {
4143  // Output the time taken to solve the problem using GMRES
4144  oomph_info << "Time for solve with FGMRES [sec]: "
4145  << this->Solution_time << std::endl;
4146  }
4147 
4148  // As we've met the tolerance for the solver and everything that should
4149  // be documented, has been, finish using the solver
4150  return;
4151  } // if (resid<Tolerance)
4152  } // for (unsigned j=0;j<Max_iter;j++)
4153 
4154  // Store the number of iterations taken
4155  this->Iterations = this->Max_iter;
4156 
4157  // Only update if we actually did something
4158  if (this->Max_iter > 0)
4159  {
4160  // Update the result vector using the result, x=x_0+V_m*y (where V_m
4161  // is given by v here)
4162  update(this->Max_iter - 1, hessenberg, s, block_z, block_solution);
4163 
4164  // Copy the vectors in block_solution to solution
4165  this->return_block_vectors(block_solution, solution);
4166  }
4167 
4168  // Compute the current residual
4169  beta = 0.0;
4170 
4171  // Get access to the values pointer (real)
4172  norm_r = block_r[0].norm();
4173 
4174  // Get access to the values pointer (imaginary)
4175  norm_c = block_r[1].norm();
4176 
4177  // Calculate the full norm
4178  beta = sqrt(pow(norm_r, 2.0) + pow(norm_c, 2.0));
4179 
4180  // Calculate the relative residual
4181  resid = beta / normb;
4182 
4183  // If the relative residual lies within tolerance
4184  if (resid < this->Tolerance)
4185  {
4186  // If time documentation is enabled
4187  if (this->Doc_time)
4188  {
4189  // Notify the user
4190  oomph_info << "\nFGMRES converged (2). Normalised residual norm: "
4191  << resid << "\nNumber of iterations to convergence: "
4192  << this->Iterations << "\n"
4193  << std::endl;
4194  }
4195 
4196  // End the timer
4197  double t_end = TimingHelpers::timer();
4198 
4199  // Calculate the time taken for the solver
4200  this->Solution_time = t_end - t_start;
4201 
4202  // If time documentation is enabled
4203  if (this->Doc_time)
4204  {
4205  oomph_info << "Time for solve with FGMRES [sec]: "
4206  << this->Solution_time << std::endl;
4207  }
4208  return;
4209  }
4210 
4211  // Otherwise GMRES failed convergence
4212  oomph_info << "\nFGMRES did not converge to required tolerance! "
4213  << "\nReturning with normalised residual norm: " << resid
4214  << "\nafter " << this->Max_iter << " iterations.\n"
4215  << std::endl;
4216 
4217  // Throw an error if requested
4218  if (this->Throw_error_after_max_iter)
4219  {
4220  std::string err = "Solver failed to converge and you requested an error";
4221  err += " on convergence failures.";
4222  throw OomphLibError(
4223  err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
4224  }
4225 
4226  // Finish using the solver
4227  return;
4228  } // End of solve_helper
4229 } // End of namespace oomph
4230 
4231 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
Block Preconditioner base class. The block structure of the overall problem is determined from the Me...
void return_block_vectors(const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const
Takes the vector of block vectors, s, and copies its entries into the naturally ordered vector,...
void get_block(const unsigned &i, const unsigned &j, MATRIX &output_matrix, const bool &ignore_replacement_block=false) const
Put block (i,j) into output_matrix. This block accounts for any coarsening of dof types and any repla...
unsigned nblock_types() const
Return the number of block types.
MATRIX * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*,...
void get_block_vectors(const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
Takes the naturally ordered vector and rearranges it into a vector of sub vectors corresponding to th...
void set_nmesh(const unsigned &n)
Specify the number of meshes required by this block preconditioner. Note: elements in different meshe...
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
void set_mesh(const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Set the i-th mesh for this block preconditioner. Note: The method set_nmesh(...) must be called befor...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1002
Vector< double > diagonal_entries() const
returns a Vector of diagonal entries of this matrix. This only works with square matrices....
Definition: matrices.cc:3465
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void complex_solve_helper(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
This is where the actual work is done.
~ComplexDampedJacobi()
Empty destructor.
void operator=(const ComplexDampedJacobi &)=delete
Broken assignment operator.
unsigned iterations() const
Number of iterations taken.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrices pointed to by Matrix_real_pt and Matrix_imag_pt can be delet...
void complex_smoother_solve(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
void complex_smoother_setup(Vector< CRDoubleMatrix * > helmholtz_matrix_pt)
Setup: Pass pointer to the matrix and store in cast form.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
Vector< double > Matrix_diagonal_inv_sq
Vector whose j'th entry is given by 1/(A_r(j,j)^2 + A_c(j,j)^2)
ComplexDampedJacobi(const double &omega=0.5)
Constructor (empty)
Vector< double > Matrix_diagonal_real
Vector containing the diagonal entries of A_r (real(A))
Vector< double > Matrix_diagonal_imag
Vector containing the diagonal entries of A_c (imag(A))
CRDoubleMatrix * Matrix_imag_pt
Pointer to the real part of the system matrix.
double & omega()
Get access to the value of Omega (lvalue)
double Omega
Damping factor.
CRDoubleMatrix * Matrix_real_pt
Pointer to the real part of the system matrix.
void solve(Problem *const &problem_pt, DoubleVector &result)
Use damped Jacobi iteration as an IterativeLinearSolver: This obtains the Jacobian matrix J and the r...
ComplexDampedJacobi(const ComplexDampedJacobi &)=delete
Broken copy constructor.
unsigned Iterations
Number of iterations taken.
void calculate_omega(const double &k, const double &h)
Function to calculate the value of Omega by passing in the value of k and h [see Elman et al....
The GMRES method rewritten for complex matrices.
void generate_plane_rotation(std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)
Helper function: Generate a plane rotation. This is done by finding the value of (i....
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
unsigned iterations() const
Number of iterations taken.
void complex_smoother_setup(Vector< CRDoubleMatrix * > helmholtz_matrix_pt)
Setup: Pass pointer to the matrix and store in cast form.
unsigned Iterations
Number of iterations taken.
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
ComplexGMRES(const ComplexGMRES &)=delete
Broken copy constructor.
void update(const unsigned &k, const Vector< Vector< std::complex< double >>> &hessenberg, const Vector< std::complex< double >> &s, const Vector< Vector< DoubleVector >> &v, Vector< DoubleVector > &x)
Helper function to update the result vector.
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void complex_solve_helper(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
This is where the actual work is done.
void complex_smoother_solve(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
void apply_plane_rotation(std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)
Helper function: Apply plane rotation. This is done using the update:
ComplexGMRES()
Constructor.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
bool Matrix_can_be_deleted
Boolean flag to indicate if the real and imaginary system matrices can be deleted.
void operator=(const ComplexGMRES &)=delete
Broken assignment operator.
~ComplexGMRES()
Empty destructor.
Vector< CRDoubleMatrix * > Matrices_storage_pt
Vector of pointers to the real and imaginary part of the system matrix.
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
bool distributed() const
distribution is serial or distributed
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned nrow() const
access function to the number of global rows.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
Abstract base class for matrices of doubles – adds abstract interfaces for solving,...
Definition: matrices.h:261
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void initialise(const double &v)
initialise the whole vector with value v
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
double norm() const
compute the 2 norm of this vector
bool built() const
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void update(const unsigned &k, const Vector< Vector< std::complex< double >>> &hessenberg, const Vector< std::complex< double >> &s, const Vector< Vector< DoubleVector >> &z_m, Vector< DoubleVector > &x)
Helper function to update the result vector.
HelmholtzFGMRESMG(const HelmholtzFGMRESMG &)=delete
Broken copy constructor.
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
virtual ~HelmholtzFGMRESMG()
Destructor (cleanup storage)
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
HelmholtzFGMRESMG()
Constructor (empty)
void operator=(const HelmholtzFGMRESMG &)=delete
Broken assignment operator.
void set_preconditioner_LHS()
Overloaded function to let the user know that left preconditioning is not possible with FGMRES,...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
Vector< CRDoubleMatrix * > Matrices_storage_pt
Vector of pointers to the real and imaginary part of the system matrix.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Implementation of the pure virtual base class function. The function has been broken because this is ...
void apply_plane_rotation(std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)
Helper function: Apply plane rotation. This is done using the update:
HelmholtzGMRESMG(const HelmholtzGMRESMG &)=delete
Broken copy constructor.
void update(const unsigned &k, const Vector< Vector< std::complex< double >>> &hessenberg, const Vector< std::complex< double >> &s, const Vector< Vector< DoubleVector >> &v, Vector< DoubleVector > &x)
Helper function to update the result vector.
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void set_preconditioner_LHS()
Set left preconditioning (the default)
virtual ~HelmholtzGMRESMG()
Destructor (cleanup storage)
void complex_matrix_multiplication(Vector< CRDoubleMatrix * > const matrices_pt, const Vector< DoubleVector > &x, Vector< DoubleVector > &soln)
Helper function to calculate a complex matrix-vector product. Assumes the matrix has been provided as...
void set_preconditioner_RHS()
Enable right preconditioning.
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
unsigned iterations() const
Number of iterations taken.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
unsigned Iterations
Number of iterations taken.
void setup()
Implementation of the pure virtual base class function. This accompanies the preconditioner_solve fun...
void operator=(const HelmholtzGMRESMG &)=delete
Broken assignment operator.
void generate_plane_rotation(std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)
Helper function: Generate a plane rotation. This is done by finding the value of (i....
Helmholtz smoother class: The smoother class is designed for the Helmholtz equation to be used in con...
void check_validity_of_solve_helper_inputs(CRDoubleMatrix *const &real_matrix_pt, CRDoubleMatrix *const &imag_matrix_pt, const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution, const double &n_dof)
Self-test to check that all the dimensions of the inputs to solve helper are consistent and everythin...
bool Use_as_smoother
When a derived class object is being used as a smoother in the MG algorithm the residual norm does no...
virtual void complex_smoother_setup(Vector< CRDoubleMatrix * > matrix_pt)=0
Setup the smoother for the matrix specified by the pointer.
HelmholtzSmoother()
Empty constructor.
virtual void complex_smoother_solve(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &result)=0
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
virtual ~HelmholtzSmoother()
Virtual empty destructor.
void complex_matrix_multiplication(Vector< CRDoubleMatrix * > matrices_pt, const Vector< DoubleVector > &x, Vector< DoubleVector > &soln)
Helper function to calculate a complex matrix-vector product. Assumes the matrix has been provided as...
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
double Jacobian_setup_time
Jacobian setup time.
double & tolerance()
Access to convergence tolerance.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:77
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:73
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
static OomphCommunicator * communicator_pt()
access to global communicator. This is the oomph-lib equivalent of MPI_COMM_WORLD
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
virtual void set_comm_pt(const OomphCommunicator *const comm_pt)
Set the communicator pointer.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
virtual void set_matrix_pt(DoubleMatrixBase *matrix_pt)
Set the matrix pointer.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1678
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
Definition: problem.cc:3921
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1246
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:167
unsigned Max_iter
Max. # of Newton iterations.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const double Pi
50 digits from maple
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...