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