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