iterative_linear_solver.cc
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-2022 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// The actual solve functions for Iterative solvers.
27
28
29// Config header generated by autoconfig
30#ifdef HAVE_CONFIG_H
31#include <oomph-lib-config.h>
32#endif
33
34// Oomph-lib includes
36
37// Required to force_ get templated builds of iterative solvers for
38// sumofmatrices class.
39#include "sum_of_matrices.h"
40
41
42namespace oomph
43{
44 //==================================================================
45 /// Default preconditioner for iterative solvers: The base
46 /// class for preconditioners is a fully functional (if trivial!)
47 /// preconditioner.
48 //==================================================================
50
51
52 /// ////////////////////////////////////////////////////////////////
53 /// ////////////////////////////////////////////////////////////////
54 /// ////////////////////////////////////////////////////////////////
55
56
57 //==================================================================
58 /// Re-solve the system defined by the last assembled Jacobian
59 /// and the rhs vector specified here. Solution is returned in
60 /// the vector result.
61 //==================================================================
62 template<typename MATRIX>
64 {
65 // We are re-solving
66 Resolving = true;
67
68#ifdef PARANOID
69 if (Matrix_pt == 0)
70 {
71 throw OomphLibError("No matrix was stored -- cannot re-solve",
72 OOMPH_CURRENT_FUNCTION,
73 OOMPH_EXCEPTION_LOCATION);
74 }
75#endif
76
77 // Call linear algebra-style solver
78 this->solve(Matrix_pt, rhs, result);
79
80 // Reset re-solving flag
81 Resolving = false;
82 }
83
84 //==================================================================
85 /// Solver: Takes pointer to problem and returns the results vector
86 /// which contains the solution of the linear system defined by
87 /// the problem's fully assembled Jacobian and residual vector.
88 //==================================================================
89 template<typename MATRIX>
90 void BiCGStab<MATRIX>::solve(Problem* const& problem_pt, DoubleVector& result)
91 {
92 // Initialise timer
93#ifdef OOMPH_HAS_MPI
94 double t_start = TimingHelpers::timer();
95#else
96 clock_t t_start = clock();
97#endif
98
99 // We're not re-solving
100 Resolving = false;
101
102 // Get rid of any previously stored data
104
105 // Get Jacobian matrix in format specified by template parameter
106 // and nonlinear residual vector
107 Matrix_pt = new MATRIX;
108 DoubleVector f;
109 problem_pt->get_jacobian(f, *Matrix_pt);
110
111
112 // We've made the matrix, we can delete it...
113 Matrix_can_be_deleted = true;
114
115 // Doc time for setup
116#ifdef OOMPH_HAS_MPI
117 double t_end = TimingHelpers::timer();
118 Jacobian_setup_time = t_end - t_start;
119#else
120 clock_t t_end = clock();
121 Jacobian_setup_time = double(t_end - t_start) / CLOCKS_PER_SEC;
122#endif
123
124 if (Doc_time)
125 {
126 oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
127 << std::endl;
128 }
129
130 // set the distribution
131 if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
132 {
133 // the solver has the same distribution as the matrix if possible
134 this->build_distribution(
135 dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt)
136 ->distribution_pt());
137 }
138 else
139 {
140 // the solver has the same distribution as the RHS
141 this->build_distribution(f.distribution_pt());
142 }
143
144 // Call linear algebra-style solver
145 if ((result.built()) &&
146 (!(*result.distribution_pt() == *this->distribution_pt())))
147 {
148 LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
149 result.build(this->distribution_pt(), 0.0);
150 this->solve_helper(Matrix_pt, f, result);
151 result.redistribute(&temp_global_dist);
152 }
153 else
154 {
155 this->solve_helper(Matrix_pt, f, result);
156 }
157
158 // Kill matrix unless it's still required for resolve
159 if (!Enable_resolve) clean_up_memory();
160 };
161
162
163 //==================================================================
164 /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
165 /// and returns the solution of the linear system.
166 /// Algorithm and variable names based on "Numerical Linear Algebra
167 /// for High-Performance Computers" by Dongarra, Duff, Sorensen & van
168 /// der Vorst. SIAM (1998), page 185.
169 //==================================================================
170 template<typename MATRIX>
172 const DoubleVector& rhs,
173 DoubleVector& solution)
174 {
175#ifdef PARANOID
176 // check that the rhs vector is setup
177 if (!rhs.built())
178 {
179 std::ostringstream error_message_stream;
180 error_message_stream << "The vectors rhs must be setup";
181 throw OomphLibError(error_message_stream.str(),
182 OOMPH_CURRENT_FUNCTION,
183 OOMPH_EXCEPTION_LOCATION);
184 }
185
186 // check that the matrix is square
187 if (matrix_pt->nrow() != matrix_pt->ncol())
188 {
189 std::ostringstream error_message_stream;
190 error_message_stream << "The matrix at matrix_pt must be square.";
191 throw OomphLibError(error_message_stream.str(),
192 OOMPH_CURRENT_FUNCTION,
193 OOMPH_EXCEPTION_LOCATION);
194 }
195
196 // check that the matrix and the rhs vector have the same nrow()
197 if (matrix_pt->nrow() != rhs.nrow())
198 {
199 std::ostringstream error_message_stream;
200 error_message_stream
201 << "The matrix and the rhs vector must have the same number of rows.";
202 throw OomphLibError(error_message_stream.str(),
203 OOMPH_CURRENT_FUNCTION,
204 OOMPH_EXCEPTION_LOCATION);
205 }
206
207 // if the matrix is distributable then it too should have the same
208 // communicator as the rhs vector
209 DistributableLinearAlgebraObject* dist_matrix_pt =
210 dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
211 if (dist_matrix_pt != 0)
212 {
213 if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
214 {
215 std::ostringstream error_message_stream;
216 error_message_stream
217 << "The matrix matrix_pt must have the same communicator as the "
218 "vectors"
219 << " rhs and result must have the same communicator";
220 throw OomphLibError(error_message_stream.str(),
221 OOMPH_CURRENT_FUNCTION,
222 OOMPH_EXCEPTION_LOCATION);
223 }
224 }
225 // if the matrix is not distributable then it the rhs vector should not be
226 // distributed
227 else
228 {
229 if (rhs.distribution_pt()->distributed())
230 {
231 std::ostringstream error_message_stream;
232 error_message_stream
233 << "The matrix (matrix_pt) is not distributable and therefore the rhs"
234 << " vector must not be distributed";
235 throw OomphLibError(error_message_stream.str(),
236 OOMPH_CURRENT_FUNCTION,
237 OOMPH_EXCEPTION_LOCATION);
238 }
239 }
240 // if the result vector is setup then check it has the same distribution
241 // as the rhs
242 if (solution.built())
243 {
244 if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
245 {
246 std::ostringstream error_message_stream;
247 error_message_stream << "The solution vector distribution has been "
248 "setup; it must have the "
249 << "same distribution as the rhs vector.";
250 throw OomphLibError(error_message_stream.str(),
251 OOMPH_CURRENT_FUNCTION,
252 OOMPH_EXCEPTION_LOCATION);
253 }
254 }
255#endif
256
257 // setup the solution if it is not
258 if (!solution.distribution_pt()->built())
259 {
260 solution.build(this->distribution_pt(), 0.0);
261 }
262 // zero
263 else
264 {
265 solution.initialise(0.0);
266 }
267
268 // Get number of dofs
269 // unsigned n_dof=rhs.size();
270 unsigned nrow_local = this->nrow_local();
271
272 // Time solver
273#ifdef OOMPH_HAS_MPI
274 double t_start = TimingHelpers::timer();
275#else
276 clock_t t_start = clock();
277#endif
278
279 // Initialise: Zero initial guess so the initial residual is
280 // equal to the RHS, i.e. the nonlinear residual
281 DoubleVector residual(rhs);
282 double residual_norm = residual.norm();
283 double rhs_norm = residual_norm;
284 if (rhs_norm == 0.0) rhs_norm = 1.0;
285 DoubleVector x(rhs.distribution_pt(), 0.0);
286
287 // Hat residual by copy operation
288 DoubleVector r_hat(residual);
289
290 // Normalised residual
291 double normalised_residual_norm = residual_norm / rhs_norm;
292
293 // if required will document convergence history to screen or file (if
294 // stream open)
295 if (Doc_convergence_history)
296 {
297 if (!Output_file_stream.is_open())
298 {
299 oomph_info << 0 << " " << normalised_residual_norm << std::endl;
300 }
301 else
302 {
303 Output_file_stream << 0 << " " << normalised_residual_norm << std::endl;
304 }
305 }
306
307 // Check immediate convergence
308 if (normalised_residual_norm < Tolerance)
309 {
310 if (Doc_time)
311 {
312 oomph_info << "BiCGStab converged immediately" << std::endl;
313 }
314 solution = x;
315
316 // Doc time for solver
317 double t_end = TimingHelpers::timer();
318 Solution_time = t_end - t_start;
319
320 if (Doc_time)
321 {
322 oomph_info << "Time for solve with BiCGStab [sec]: " << Solution_time
323 << std::endl;
324 }
325 return;
326 }
327
328 // Setup preconditioner only if we're not re-solving
329 if (!Resolving)
330 {
331 // only setup the preconditioner if required
332 if (Setup_preconditioner_before_solve)
333 {
334 // Setup preconditioner from the Jacobian matrix
335 double t_start_prec = TimingHelpers::timer();
336
337 preconditioner_pt()->setup(matrix_pt);
338
339 // Doc time for setup of preconditioner
340 double t_end_prec = TimingHelpers::timer();
341 Preconditioner_setup_time = t_end_prec - t_start_prec;
342
343 if (Doc_time)
344 {
345 oomph_info << "Time for setup of preconditioner [sec]: "
346 << Preconditioner_setup_time << std::endl;
347 }
348 }
349 }
350 else
351 {
352 if (Doc_time)
353 {
354 oomph_info << "Setup of preconditioner is bypassed in resolve mode"
355 << std::endl;
356 }
357 }
358
359 // Some auxiliary variables
360 double rho_prev = 0.0;
361 double alpha = 0.0;
362 double omega = 0.0;
363 double rho, beta, dot_prod, dot_prod_tt, dot_prod_ts;
364 double s_norm, r_norm;
365
366 // Some vectors
367 DoubleVector p(this->distribution_pt(), 0.0),
368 p_hat(this->distribution_pt(), 0.0), v(this->distribution_pt(), 0.0),
369 z(this->distribution_pt(), 0.0), t(this->distribution_pt(), 0.0),
370 s(this->distribution_pt(), 0.0);
371
372 // Loop over max. number of iterations
373 for (unsigned iter = 1; iter <= Max_iter; iter++)
374 {
375 // Dot product for rho
376 rho = r_hat.dot(residual);
377
378 // Breakdown?
379 if (rho == 0.0)
380 {
381 oomph_info << "BiCGStab has broken down after " << iter << " iterations"
382 << std::endl;
383 oomph_info << "Returning with current normalised residual of "
384 << normalised_residual_norm << std::endl;
385 }
386
387 // First step is different
388 if (iter == 1)
389 {
390 p = residual;
391 }
392 else
393 {
394 beta = (rho / rho_prev) * (alpha / omega);
395 for (unsigned i = 0; i < nrow_local; i++)
396 {
397 p[i] = residual[i] + beta * (p[i] - omega * v[i]);
398 }
399 }
400
401 // Apply precondtitioner: p_hat=P^-1*p
402 preconditioner_pt()->preconditioner_solve(p, p_hat);
403
404 // Matrix vector product: v=A*p_hat
405 matrix_pt->multiply(p_hat, v);
406 dot_prod = r_hat.dot(v);
407 alpha = rho / dot_prod;
408 for (unsigned i = 0; i < nrow_local; i++)
409 {
410 s[i] = residual[i] - alpha * v[i];
411 }
412 s_norm = s.norm();
413
414 // Normalised residual
415 normalised_residual_norm = s_norm / rhs_norm;
416
417 // if required will document convergence history to screen or file (if
418 // stream open)
419 if (Doc_convergence_history)
420 {
421 if (!Output_file_stream.is_open())
422 {
423 oomph_info << double(iter - 0.5) << " " << normalised_residual_norm
424 << std::endl;
425 }
426 else
427 {
428 Output_file_stream << double(iter - 0.5) << " "
429 << normalised_residual_norm << std::endl;
430 }
431 }
432
433 // Converged?
434 if (normalised_residual_norm < Tolerance)
435 {
436 for (unsigned i = 0; i < nrow_local; i++)
437 {
438 solution[i] = x[i] + alpha * p_hat[i];
439 }
440
441 if (Doc_time)
442 {
443 oomph_info << std::endl;
444 oomph_info << "BiCGStab converged. Normalised residual norm: "
445 << normalised_residual_norm << std::endl;
446 oomph_info << "Number of iterations to convergence: " << iter
447 << std::endl;
448 oomph_info << std::endl;
449 }
450
451 // Store number of iterations taken
452 Iterations = iter;
453
454 // Doc time for solver
455 double t_end = TimingHelpers::timer();
456 Solution_time = t_end - t_start;
457
458 if (Doc_time)
459 {
460 oomph_info << "Time for solve with BiCGStab [sec]: " << Solution_time
461 << std::endl;
462 }
463
464 return;
465 }
466
467 // Apply precondtitioner: z=P^-1*s
468 preconditioner_pt()->preconditioner_solve(s, z);
469
470 // Matrix vector product: t=A*z
471 matrix_pt->multiply(z, t);
472 dot_prod_ts = t.dot(s);
473 dot_prod_tt = t.dot(t);
474 omega = dot_prod_ts / dot_prod_tt;
475 for (unsigned i = 0; i < nrow_local; i++)
476 {
477 x[i] += alpha * p_hat[i] + omega * z[i];
478 residual[i] = s[i] - omega * t[i];
479 }
480 r_norm = residual.norm();
481 rho_prev = rho;
482
483 // Check convergence again
484 normalised_residual_norm = r_norm / rhs_norm;
485
486 // if required will document convergence history to screen or file (if
487 // stream open)
488 if (Doc_convergence_history)
489 {
490 if (!Output_file_stream.is_open())
491 {
492 oomph_info << iter << " " << normalised_residual_norm << std::endl;
493 }
494 else
495 {
496 Output_file_stream << iter << " " << normalised_residual_norm
497 << std::endl;
498 }
499 }
500
501
502 if (normalised_residual_norm < Tolerance)
503 {
504 if (Doc_time)
505 {
506 oomph_info << std::endl;
507 oomph_info << "BiCGStab converged. Normalised residual norm: "
508 << normalised_residual_norm << std::endl;
509 oomph_info << "Number of iterations to convergence: " << iter
510 << std::endl;
511 oomph_info << std::endl;
512 }
513 solution = x;
514
515 // Store the number of itertions taken.
516 Iterations = iter;
517
518 // Doc time for solver
519 double t_end = TimingHelpers::timer();
520 Solution_time = t_end - t_start;
521
522 if (Doc_time)
523 {
524 oomph_info << "Time for solve with BiCGStab [sec]: " << Solution_time
525 << std::endl;
526 }
527 return;
528 }
529
530
531 // Breakdown: Omega has to be >0 for to be able to continue
532 if (omega == 0.0)
533 {
534 oomph_info << std::endl;
535 oomph_info << "BiCGStab breakdown with omega=0.0. "
536 << "Normalised residual norm: " << normalised_residual_norm
537 << std::endl;
538 oomph_info << "Number of iterations so far: " << iter << std::endl;
539 oomph_info << std::endl;
540 solution = x;
541
542 // Store the number of itertions taken.
543 Iterations = iter;
544
545 // Doc time for solver
546 double t_end = TimingHelpers::timer();
547 Solution_time = t_end - t_start;
548
549 if (Doc_time)
550 {
551 oomph_info << "Time for solve with BiCGStab [sec]: " << Solution_time
552 << std::endl;
553 }
554 return;
555 }
556
557
558 } // end of iteration loop
559
560
561 // No convergence
562 oomph_info << std::endl;
563 oomph_info << "BiCGStab did not converge to required tolerance! "
564 << std::endl;
565 oomph_info << "Returning with normalised residual norm: "
566 << normalised_residual_norm << std::endl;
567 oomph_info << "after " << Max_iter << " iterations." << std::endl;
568 oomph_info << std::endl;
569
570 solution = x;
571
572 // Doc time for solver
573 double t_end = TimingHelpers::timer();
574 Solution_time = t_end - t_start;
575
576 if (Doc_time)
577 {
578 oomph_info << "Time for solve with BiCGStab [sec]: " << Solution_time
579 << std::endl;
580 }
581
582 if (Throw_error_after_max_iter)
583 {
584 std::string err = "Solver failed to converge and you requested an error";
585 err += " on convergence failures.";
586 throw OomphLibError(
587 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
588 }
589 }
590
591
592 /// ///////////////////////////////////////////////////////////////////
593 /// ///////////////////////////////////////////////////////////////////
594 /// ///////////////////////////////////////////////////////////////////
595
596
597 //==================================================================
598 /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
599 /// and returns the solution of the linear system.
600 /// Algorithm and variable names based on "Matrix Computations,
601 /// 2nd Ed." Golub & van Loan, John Hopkins University Press(1989),
602 /// page 529.
603 //==================================================================
604 template<typename MATRIX>
606 const DoubleVector& rhs,
607 DoubleVector& solution)
608 {
609#ifdef PARANOID
610 // check that the rhs vector is setup
611 if (!rhs.built())
612 {
613 std::ostringstream error_message_stream;
614 error_message_stream << "The vectors rhs must be setup";
615 throw OomphLibError(error_message_stream.str(),
616 OOMPH_CURRENT_FUNCTION,
617 OOMPH_EXCEPTION_LOCATION);
618 }
619
620 // check that the matrix is square
621 if (matrix_pt->nrow() != matrix_pt->ncol())
622 {
623 std::ostringstream error_message_stream;
624 error_message_stream << "The matrix at matrix_pt must be square.";
625 throw OomphLibError(error_message_stream.str(),
626 OOMPH_CURRENT_FUNCTION,
627 OOMPH_EXCEPTION_LOCATION);
628 }
629
630 // check that the matrix and the rhs vector have the same nrow()
631 if (matrix_pt->nrow() != rhs.nrow())
632 {
633 std::ostringstream error_message_stream;
634 error_message_stream
635 << "The matrix and the rhs vector must have the same number of rows.";
636 throw OomphLibError(error_message_stream.str(),
637 OOMPH_CURRENT_FUNCTION,
638 OOMPH_EXCEPTION_LOCATION);
639 }
640
641 // if the matrix is distributable then it too should have the same
642 // communicator as the rhs vector
643 DistributableLinearAlgebraObject* dist_matrix_pt =
644 dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
645 if (dist_matrix_pt != 0)
646 {
647 if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
648 {
649 std::ostringstream error_message_stream;
650 error_message_stream
651 << "The matrix matrix_pt must have the same communicator as the "
652 "vectors"
653 << " rhs and result must have the same communicator";
654 throw OomphLibError(error_message_stream.str(),
655 OOMPH_CURRENT_FUNCTION,
656 OOMPH_EXCEPTION_LOCATION);
657 }
658 }
659 // if the matrix is not distributable then it the rhs vector should not be
660 // distributed
661 else
662 {
663 if (rhs.distribution_pt()->distributed())
664 {
665 std::ostringstream error_message_stream;
666 error_message_stream
667 << "The matrix (matrix_pt) is not distributable and therefore the rhs"
668 << " vector must not be distributed";
669 throw OomphLibError(error_message_stream.str(),
670 OOMPH_CURRENT_FUNCTION,
671 OOMPH_EXCEPTION_LOCATION);
672 }
673 }
674 // if the result vector is setup then check it has the same distribution
675 // as the rhs
676 if (solution.built())
677 {
678 if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
679 {
680 std::ostringstream error_message_stream;
681 error_message_stream << "The solution vector distribution has been "
682 "setup; it must have the "
683 << "same distribution as the rhs vector.";
684 throw OomphLibError(error_message_stream.str(),
685 OOMPH_CURRENT_FUNCTION,
686 OOMPH_EXCEPTION_LOCATION);
687 }
688 }
689#endif
690
691 // setup the solution if it is not
692 if (!solution.distribution_pt()->built())
693 {
694 solution.build(this->distribution_pt(), 0.0);
695 }
696 // zero
697 else
698 {
699 solution.initialise(0.0);
700 }
701
702 // Get number of dofs
703 // unsigned n_dof=rhs.size();
704 unsigned nrow_local = this->nrow_local();
705
706 // Initialise counter
707 unsigned counter = 0;
708
709 // Time solver
710 double t_start = TimingHelpers::timer();
711
712 // Initialise: Zero initial guess so the initial residual is
713 // equal to the RHS
714 DoubleVector x(this->distribution_pt(), 0.0);
715 DoubleVector residual(rhs);
716 double residual_norm = residual.norm();
717 double rhs_norm = residual_norm;
718 if (rhs_norm == 0.0) rhs_norm = 1.0;
719
720 // Normalised residual
721 double normalised_residual_norm = residual_norm / rhs_norm;
722
723 // if required will document convergence history to screen or file (if
724 // stream open)
725 if (Doc_convergence_history)
726 {
727 if (!Output_file_stream.is_open())
728 {
729 oomph_info << 0 << " " << normalised_residual_norm << std::endl;
730 }
731 else
732 {
733 Output_file_stream << 0 << " " << normalised_residual_norm << std::endl;
734 }
735 }
736
737 // Check immediate convergence
738 if (normalised_residual_norm < Tolerance)
739 {
740 if (Doc_time)
741 {
742 oomph_info << "CG converged immediately" << std::endl;
743 }
744 solution = x;
745
746 // Doc time for solver
747 double t_end = TimingHelpers::timer();
748 Solution_time = t_end - t_start;
749
750 if (Doc_time)
751 {
752 oomph_info << "Time for solve with CG [sec]: " << Solution_time
753 << std::endl;
754 }
755 return;
756 }
757
758
759 // Setup preconditioner only if we're not re-solving
760 if (!Resolving)
761 {
762 // only setup the preconditioner if required
763 if (Setup_preconditioner_before_solve)
764 {
765 // Setup preconditioner from the Jacobian matrix
766 double t_start_prec = TimingHelpers::timer();
767
768 preconditioner_pt()->setup(matrix_pt);
769
770 // Doc time for setup of preconditioner
771 double t_end_prec = TimingHelpers::timer();
772 Preconditioner_setup_time = t_end_prec - t_start_prec;
773
774 if (Doc_time)
775 {
776 oomph_info << "Time for setup of preconditioner [sec]: "
777 << Preconditioner_setup_time << std::endl;
778 }
779 }
780 }
781 else
782 {
783 if (Doc_time)
784 {
785 oomph_info << "Setup of preconditioner is bypassed in resolve mode"
786 << std::endl;
787 }
788 }
789
790
791 // Auxiliary vectors
792 // Vector<double> z(n_dof),p(n_dof),jacobian_times_p(n_dof,0.0);
793 DoubleVector z(this->distribution_pt(), 0.0),
794 p(this->distribution_pt(), 0.0),
795 jacobian_times_p(this->distribution_pt(), 0.0);
796
797 // Auxiliary values
798 double alpha, beta, rz;
799 double prev_rz = 0.0;
800
801 // Main iteration
802 while ((normalised_residual_norm > Tolerance) && (counter != Max_iter))
803 {
804 // Apply precondtitioner: z=P^-1*r
805 preconditioner_pt()->preconditioner_solve(residual, z);
806
807 // P vector is computed differently for first and subsequent steps
808 if (counter == 0)
809 {
810 p = z;
811 rz = residual.dot(z);
812 }
813 // Subsequent steps
814 else
815 {
816 rz = residual.dot(z);
817 beta = rz / prev_rz;
818 for (unsigned i = 0; i < nrow_local; i++)
819 {
820 p[i] = z[i] + beta * p[i];
821 }
822 }
823
824
825 // Matrix vector product
826 matrix_pt->multiply(p, jacobian_times_p);
827 double pq = p.dot(jacobian_times_p);
828 alpha = rz / pq;
829
830 // Update
831 prev_rz = rz;
832 for (unsigned i = 0; i < nrow_local; i++)
833 {
834 x[i] += alpha * p[i];
835 residual[i] -= alpha * jacobian_times_p[i];
836 }
837
838
839 // Calculate the 2norm
840 residual_norm = residual.norm();
841
842 // Difference between the initial and current 2norm residual
843 normalised_residual_norm = residual_norm / rhs_norm;
844
845
846 // if required will document convergence history to screen or file (if
847 // stream open)
848 if (Doc_convergence_history)
849 {
850 if (!Output_file_stream.is_open())
851 {
852 oomph_info << counter << " " << normalised_residual_norm << std::endl;
853 }
854 else
855 {
856 Output_file_stream << counter << " " << normalised_residual_norm
857 << std::endl;
858 }
859 }
860
861 counter = counter + 1;
862
863 } // end while
864
865
866 if (counter >= Max_iter)
867 {
868 oomph_info << std::endl;
869 oomph_info << "CG did not converge to required tolerance! " << std::endl;
870 oomph_info << "Returning with normalised residual norm: "
871 << normalised_residual_norm << std::endl;
872 oomph_info << "after " << counter << " iterations." << std::endl;
873 oomph_info << std::endl;
874 }
875 else
876 {
877 if (Doc_time)
878 {
879 oomph_info << std::endl;
880 oomph_info << "CG converged. Normalised residual norm: "
881 << normalised_residual_norm << std::endl;
882 oomph_info << "Number of iterations to convergence: " << counter
883 << std::endl;
884 oomph_info << std::endl;
885 }
886 }
887
888
889 // Store number if iterations taken
890 Iterations = counter;
891
892 // Copy result back
893 solution = x;
894
895 // Doc time for solver
896 double t_end = TimingHelpers::timer();
897 Solution_time = t_end - t_start;
898
899 if (Doc_time)
900 {
901 oomph_info << "Time for solve with CG [sec]: " << Solution_time
902 << std::endl;
903 }
904
905 if ((counter >= Max_iter) && (Throw_error_after_max_iter))
906 {
907 std::string err = "Solver failed to converge and you requested an error";
908 err += " on convergence failures.";
909 throw OomphLibError(
910 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
911 }
912
913 } // end CG
914
915
916 //==================================================================
917 /// Re-solve the system defined by the last assembled Jacobian
918 /// and the rhs vector specified here. Solution is returned in
919 /// the vector result.
920 //==================================================================
921 template<typename MATRIX>
923 {
924 // We are re-solving
925 Resolving = true;
926
927#ifdef PARANOID
928 if (Matrix_pt == 0)
929 {
930 throw OomphLibError("No matrix was stored -- cannot re-solve",
931 OOMPH_CURRENT_FUNCTION,
932 OOMPH_EXCEPTION_LOCATION);
933 }
934#endif
935
936 // Call linear algebra-style solver
937 this->solve(Matrix_pt, rhs, result);
938
939 // Reset re-solving flag
940 Resolving = false;
941 }
942
943
944 //==================================================================
945 /// Solver: Takes pointer to problem and returns the results vector
946 /// which contains the solution of the linear system defined by
947 /// the problem's fully assembled Jacobian and residual vector.
948 //==================================================================
949 template<typename MATRIX>
950 void CG<MATRIX>::solve(Problem* const& problem_pt, DoubleVector& result)
951 {
952 // Initialise timer
953 double t_start = TimingHelpers::timer();
954
955 // We're not re-solving
956 Resolving = false;
957
958 // Get rid of any previously stored data
960
961 // Get Jacobian matrix in format specified by template parameter
962 // and nonlinear residual vector
963 Matrix_pt = new MATRIX;
964 DoubleVector f;
965 problem_pt->get_jacobian(f, *Matrix_pt);
966
967 // We've made the matrix, we can delete it...
968 Matrix_can_be_deleted = true;
969
970 // Doc time for setup
971 double t_end = TimingHelpers::timer();
972 Jacobian_setup_time = t_end - t_start;
973
974 if (Doc_time)
975 {
976 oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
977 << std::endl;
978 }
979
980 // set the distribution
981 if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
982 {
983 // the solver has the same distribution as the matrix if possible
984 this->build_distribution(
985 dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt)
986 ->distribution_pt());
987 }
988 else
989 {
990 // the solver has the same distribution as the RHS
991 this->build_distribution(f.distribution_pt());
992 }
993
994 // if the result vector is not setup
995 if (!result.distribution_pt()->built())
996 {
997 result.build(this->distribution_pt(), 0.0);
998 }
999
1000 // Call linear algebra-style solver
1001 if (!(*result.distribution_pt() == *this->distribution_pt()))
1002 {
1003 LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
1004 result.build(this->distribution_pt(), 0.0);
1005 this->solve_helper(Matrix_pt, f, result);
1006 result.redistribute(&temp_global_dist);
1007 }
1008 else
1009 {
1010 this->solve_helper(Matrix_pt, f, result);
1011 }
1012
1013 // Kill matrix unless it's still required for resolve
1014 if (!Enable_resolve) clean_up_memory();
1015 };
1016
1017
1018 /// ///////////////////////////////////////////////////////////////////
1019 /// ///////////////////////////////////////////////////////////////////
1020 /// ///////////////////////////////////////////////////////////////////
1021
1022
1023 //==================================================================
1024 /// Self-test to be called inside solve_helper to ensure
1025 /// that all inputs are consistent and everything that needs to
1026 /// be built, is.
1027 //==================================================================
1028 template<typename MATRIX>
1030 const DoubleVector& rhs,
1031 DoubleVector& solution,
1032 const double& n_dof)
1033 {
1034 // Check that if the matrix is not distributable then it should not
1035 // be distributed
1036 if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
1037 {
1038 if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
1039 ->distributed())
1040 {
1041 std::ostringstream error_message_stream;
1042 error_message_stream << "The matrix must not be distributed.";
1043 throw OomphLibError(error_message_stream.str(),
1044 OOMPH_CURRENT_FUNCTION,
1045 OOMPH_EXCEPTION_LOCATION);
1046 }
1047 } // if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)!=0)
1048 // Check that this rhs distribution is setup
1049 if (!rhs.built())
1050 {
1051 std::ostringstream error_message_stream;
1052 error_message_stream << "The rhs vector distribution must be setup.";
1053 throw OomphLibError(error_message_stream.str(),
1054 OOMPH_CURRENT_FUNCTION,
1055 OOMPH_EXCEPTION_LOCATION);
1056 }
1057 // Check that the rhs has the right number of global rows
1058 if (rhs.nrow() != n_dof)
1059 {
1060 std::ostringstream error_message_stream;
1061 error_message_stream << "RHS does not have the same dimension as "
1062 << "the linear system";
1063 throw OomphLibError(error_message_stream.str(),
1064 OOMPH_CURRENT_FUNCTION,
1065 OOMPH_EXCEPTION_LOCATION);
1066 }
1067 // Check that the rhs is not distributed
1068 if (rhs.distribution_pt()->distributed())
1069 {
1070 std::ostringstream error_message_stream;
1071 error_message_stream << "The rhs vector must not be distributed.";
1072 throw OomphLibError(error_message_stream.str(),
1073 OOMPH_CURRENT_FUNCTION,
1074 OOMPH_EXCEPTION_LOCATION);
1075 }
1076 // Check that if the result is setup it matches the distribution
1077 // of the rhs
1078 if (solution.built())
1079 {
1080 if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
1081 {
1082 std::ostringstream error_message_stream;
1083 error_message_stream << "If the result distribution is setup then it "
1084 << "must be the same as the rhs distribution";
1085 throw OomphLibError(error_message_stream.str(),
1086 OOMPH_CURRENT_FUNCTION,
1087 OOMPH_EXCEPTION_LOCATION);
1088 }
1089 } // if (solution.built())
1090 } // End of check_validity_of_solve_helper_inputs
1091
1092
1093 /// ////////////////////////////////////////////////////////////////////
1094 /// ////////////////////////////////////////////////////////////////////
1095 /// ////////////////////////////////////////////////////////////////////
1096
1097
1098 //==================================================================
1099 /// Solver: Takes pointer to problem and returns the results
1100 /// vector which contains the solution of the linear system defined
1101 /// by the problem's fully assembled Jacobian and residual vector.
1102 //==================================================================
1103 template<typename MATRIX>
1104 void GS<MATRIX>::solve(Problem* const& problem_pt, DoubleVector& result)
1105 {
1106 // Reset the Use_as_smoother_flag as the solver is not being used
1107 // as a smoother
1108 Use_as_smoother = false;
1109
1110 // Find the # of degrees of freedom (variables)
1111 unsigned n_dof = problem_pt->ndof();
1112
1113 // Initialise timer
1114 double t_start = TimingHelpers::timer();
1115
1116 // We're not re-solving
1117 Resolving = false;
1118
1119 // Get rid of any previously stored data
1121
1122 // Set up the distribution
1123 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
1124
1125 // Assign the distribution to the LinearSolver
1126 this->build_distribution(dist);
1127
1128 // Allocate space for the Jacobian matrix in format specified
1129 // by template parameter
1130 Matrix_pt = new MATRIX;
1131
1132 // Get the nonlinear residuals vector
1133 DoubleVector f;
1134
1135 // Assign the Jacobian and the residuals vector
1136 problem_pt->get_jacobian(f, *Matrix_pt);
1137
1138 // We've made the matrix, we can delete it...
1139 Matrix_can_be_deleted = true;
1140
1141 // Doc time for setup
1142 double t_end = TimingHelpers::timer();
1143 Jacobian_setup_time = t_end - t_start;
1144
1145 // If time documentation is enabled
1146 if (Doc_time)
1147 {
1148 oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
1149 << std::endl;
1150 }
1151
1152 // Call linear algebra-style solver
1153 this->solve_helper(Matrix_pt, f, result);
1154
1155 // Kill matrix unless it's still required for resolve
1156 if (!Enable_resolve) clean_up_memory();
1157 } // End of solve
1158
1159 //==================================================================
1160 /// Linear-algebra-type solver: Takes pointer to a matrix and
1161 /// rhs vector and returns the solution of the linear system.
1162 //==================================================================
1163 template<typename MATRIX>
1165 const DoubleVector& rhs,
1166 DoubleVector& solution)
1167 {
1168 // Get number of dofs
1169 unsigned n_dof = rhs.nrow();
1170
1171#ifdef PARANOID
1172 // Upcast the matrix to the appropriate type
1173 MATRIX* tmp_matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
1174
1175 // PARANOID Run the self-tests to check the inputs are correct
1176 this->check_validity_of_solve_helper_inputs<MATRIX>(
1177 tmp_matrix_pt, rhs, solution, n_dof);
1178
1179 // We don't need the pointer any more but we do still need the matrix
1180 // so just make tmp_matrix_pt a null pointer
1181 tmp_matrix_pt = 0;
1182#endif
1183
1184 // Set up the solution if it is not
1185 if (!solution.distribution_pt()->built())
1186 {
1187 // Build it!
1188 solution.build(this->distribution_pt(), 0.0);
1189 }
1190 // If the solution is already set up
1191 else
1192 {
1193 // If we're inside the multigrid solver then as we traverse up the
1194 // hierarchy we use the smoother on the updated approximate solution.
1195 // As such, we should ONLY be resetting all the values to zero if
1196 // we're NOT inside the multigrid solver
1197 if (!Use_as_smoother)
1198 {
1199 // Initialise the vector with all entries set to zero
1200 solution.initialise(0.0);
1201 }
1202 } // if (!solution.distribution_pt()->built())
1203
1204 // Initialise timer
1205 double t_start = TimingHelpers::timer();
1206
1207 // Copy the solution vector into x
1208 DoubleVector x(solution);
1209
1210 // Create a vector to hold the residual. This will only be built if
1211 // we're not inside the multigrid solver
1212 DoubleVector current_residual;
1213
1214 // Variable to hold the current residual norm. Only used if we're
1215 // not inside the multigrid solver
1216 double norm_res = 0.0;
1217
1218 // Variables to hold the initial residual norm. Only used if we're
1219 // not inside the multigrid solver
1220 double norm_f = 0.0;
1221
1222 // Initialise the value of Iterations
1223 Iterations = 0;
1224
1225 // Calculate the residual only if we're not inside the multigrid solver
1226 if (!Use_as_smoother)
1227 {
1228 // Build the residual vector
1229 current_residual.build(this->distribution_pt(), 0.0);
1230
1231 // Calculate the residual r=b-Ax
1232 matrix_pt->residual(x, rhs, current_residual);
1233
1234 // Calculate the 2-norm of the residual vector
1235 norm_res = current_residual.norm();
1236
1237 // Store the initial norm
1238 norm_f = norm_res;
1239
1240 // If required, document the convergence history to screen or file (if
1241 // the output stream open)
1242 if (Doc_convergence_history)
1243 {
1244 // If the output file stream isn't open
1245 if (!Output_file_stream.is_open())
1246 {
1247 // Output the result to screen
1248 oomph_info << Iterations << " " << norm_res << std::endl;
1249 }
1250 // If the output file stream is open
1251 else
1252 {
1253 // Document the result to file
1254 Output_file_stream << Iterations << " " << norm_res << std::endl;
1255 }
1256 } // if (Doc_convergence_history)
1257 } // if (!Use_as_smoother)
1258
1259 // Start of the main GS loop: run up to Max_iter times
1260 for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
1261 {
1262 // Loop over rows
1263 for (unsigned i = 0; i < n_dof; i++)
1264 {
1265 double dummy = rhs[i];
1266 for (unsigned j = 0; j < i; j++)
1267 {
1268 dummy -= (*(matrix_pt))(i, j) * x[j];
1269 }
1270 for (unsigned j = (i + 1); j < n_dof; j++)
1271 {
1272 dummy -= (*(matrix_pt))(i, j) * x[j];
1273 }
1274 x[i] = dummy / (*(matrix_pt))(i, i);
1275 } // for(unsigned i=0;i<n_dof;i++)
1276
1277 // Increment the value of Iterations
1278 Iterations++;
1279
1280 // Calculate the residual only if we're not inside the multigrid solver
1281 if (!Use_as_smoother)
1282 {
1283 // Get residual
1284 matrix_pt->residual(x, rhs, current_residual);
1285
1286 // Calculate the relative residual norm (i.e.
1287 // \frac{\|r_{i}\|}{\|r_{0}\|})
1288 norm_res = current_residual.norm() / norm_f;
1289
1290 // If required will document convergence history to screen or file (if
1291 // stream open)
1292 if (Doc_convergence_history)
1293 {
1294 if (!Output_file_stream.is_open())
1295 {
1296 oomph_info << Iterations << " " << norm_res << std::endl;
1297 }
1298 else
1299 {
1300 Output_file_stream << Iterations << " " << norm_res << std::endl;
1301 }
1302 } // if (Doc_convergence_history)
1303
1304 // Check the tolerance only if the residual norm is being computed
1305 if (norm_res < Tolerance)
1306 {
1307 // Break out of the for-loop
1308 break;
1309 }
1310 } // if (!Use_as_smoother)
1311 } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
1312
1313 // Calculate the residual only if we're not inside the multigrid solver
1314 if (!Use_as_smoother)
1315 {
1316 // If time documentation is enabled
1317 if (Doc_time)
1318 {
1319 oomph_info << "\nGS converged. Residual norm: " << norm_res
1320 << "\nNumber of iterations to convergence: " << Iterations
1321 << "\n"
1322 << std::endl;
1323 }
1324 } // if (!Use_as_smoother)
1325
1326 // Copy result into result
1327 solution = x;
1328
1329 // Doc. time for solver
1330 double t_end = TimingHelpers::timer();
1331 Solution_time = t_end - t_start;
1332 if (Doc_time)
1333 {
1334 oomph_info << "Time for solve with GS [sec]: " << Solution_time
1335 << std::endl;
1336 }
1337
1338 // If the solver failed to converge and the user asked for an error if
1339 // this happened
1340 if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
1341 {
1342 std::string error_message =
1343 "Solver failed to converge and you requested ";
1344 error_message += "an error on convergence failures.";
1345 throw OomphLibError(
1346 error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1347 }
1348 } // End of solve_helper
1349
1350
1351 /// ////////////////////////////////////////////////////////////////////
1352 /// ////////////////////////////////////////////////////////////////////
1353 /// ////////////////////////////////////////////////////////////////////
1354
1355
1356 //==================================================================
1357 /// Explicit template specialisation of the solver for CR
1358 /// matrices: Takes pointer to problem and returns the results
1359 /// vector which contains the solution of the linear system defined
1360 /// by the problem's fully assembled Jacobian and residual vector.
1361 //==================================================================
1362 void GS<CRDoubleMatrix>::solve(Problem* const& problem_pt,
1363 DoubleVector& result)
1364 {
1365 // Reset the Use_as_smoother_flag as the solver is not being used
1366 // as a smoother
1367 Use_as_smoother = false;
1368
1369 // Find # of degrees of freedom (variables)
1370 unsigned n_dof = problem_pt->ndof();
1371
1372 // Initialise timer
1373 double t_start = TimingHelpers::timer();
1374
1375 // We're not re-solving
1376 Resolving = false;
1377
1378 // Get rid of any previously stored data
1380
1381 // Set up the distribution
1382 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
1383
1384 // Build the distribution of the LinearSolver
1385 this->build_distribution(dist);
1386
1387 // Allocate space for the Jacobian matrix in CRDoubleMatrix format
1389
1390 // Build the matrix
1391 Matrix_pt->build(this->distribution_pt());
1392
1393 // Allocate space for the residuals vector
1394 DoubleVector f;
1395
1396 // Build it
1397 f.build(this->distribution_pt(), 0.0);
1398
1399 // Get the global Jacobian and the accompanying residuals vector
1400 problem_pt->get_jacobian(f, *Matrix_pt);
1401
1402 // Doc. time for setup
1403 double t_end = TimingHelpers::timer();
1404 Jacobian_setup_time = t_end - t_start;
1405
1406 // Run the generic setup function
1407 setup_helper(Matrix_pt);
1408
1409 // We've made the matrix, we can delete it...
1410 Matrix_can_be_deleted = true;
1411
1412 // If time documentation is enabled
1413 if (Doc_time)
1414 {
1415 // Output the time for the assembly of the Jacobian to screen
1416 oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
1417 << std::endl;
1418 }
1419
1420 // Call linear algebra-style solver
1421 solve_helper(Matrix_pt, f, result);
1422
1423 // Kill matrix unless it's still required for resolve
1425 } // End of solve
1426
1427 //==================================================================
1428 /// Explicit template specialisation of the smoother_setup
1429 /// function for CR matrices. Set up the smoother for the matrix
1430 /// specified by the pointer. This definition of the smoother_setup
1431 /// has the added functionality that it sorts the entries in the
1432 /// given matrix so that the CRDoubleMatrix implementation of the
1433 /// solver can be used.
1434 //==================================================================
1436 {
1437 // Assume the matrix has been passed in from the outside so we must
1438 // not delete it. This is needed to avoid pre- and post-smoothers
1439 // deleting the same matrix in the MG solver. If this was originally
1440 // set to TRUE then this will be sorted out in the other functions
1441 // from which this was called
1442 Matrix_can_be_deleted = false;
1443
1444 // Upcast the input matrix to system matrix to the type MATRIX
1445 Matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1446
1447 // The system matrix here is a CRDoubleMatrix. To make use of the
1448 // specific implementation of the solver for this type of matrix we
1449 // need to make sure the entries are arranged correctly
1450 Matrix_pt->sort_entries();
1451
1452 // Now get access to the vector Index_of_diagonal_entries
1453 Index_of_diagonal_entries = Matrix_pt->get_index_of_diagonal_entries();
1454
1455#ifdef PARANOID
1456 // Create a boolean variable to store the result of whether or not
1457 // the entries are in the correct order
1458 bool is_matrix_sorted = true;
1459
1460 // Check to make sure the entries have been sorted properly
1461 is_matrix_sorted = Matrix_pt->entries_are_sorted();
1462
1463 // If the entries are not sorted properly
1464 if (!is_matrix_sorted)
1465 {
1466 // Throw an error; the solver won't work unless the matrix is sorted
1467 // properly
1468 throw OomphLibError("Matrix is not sorted correctly",
1469 OOMPH_CURRENT_FUNCTION,
1470 OOMPH_EXCEPTION_LOCATION);
1471 }
1472#endif
1473 } // End of setup_helper
1474
1475 //==================================================================
1476 /// Explicit template specialisation of the solve_helper
1477 /// function for CR matrices. Exploiting the sparsity of the given
1478 /// matrix allows for a much faster iterative solver.
1479 //==================================================================
1481 const DoubleVector& rhs,
1482 DoubleVector& solution)
1483 {
1484 // Get number of dofs
1485 unsigned n_dof = rhs.nrow();
1486
1487#ifdef PARANOID
1488 // Upcast the matrix to the appropriate type
1489 CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1490
1491 // PARANOID Run the self-tests to check the inputs are correct
1492 this->check_validity_of_solve_helper_inputs<CRDoubleMatrix>(
1493 cr_matrix_pt, rhs, solution, n_dof);
1494
1495 // We don't need the pointer any more but we do still need the matrix
1496 // so just make tmp_matrix_pt a null pointer
1497 cr_matrix_pt = 0;
1498#endif
1499
1500 // Setup the solution if it is not
1501 if (!solution.distribution_pt()->built())
1502 {
1503 solution.build(this->distribution_pt(), 0.0);
1504 }
1505 // If the solution is already set up
1506 else
1507 {
1508 // If we're inside the multigrid solver then as we traverse up the
1509 // hierarchy we use the smoother on the updated approximate solution.
1510 // As such, we should ONLY be resetting all the values to zero if
1511 // we're NOT inside the multigrid solver
1512 if (!Use_as_smoother)
1513 {
1514 // Initialise the vector with all entries set to zero
1515 solution.initialise(0.0);
1516 }
1517 } // if (!solution.distribution_pt()->built())
1518
1519 // Initialise timer
1520 double t_start = TimingHelpers::timer();
1521
1522 // Copy the solution vector into x
1523 DoubleVector x(solution);
1524
1525 // Create a vector to hold the residual. This will only be built if
1526 // we're not inside the multigrid solver
1527 DoubleVector current_residual;
1528
1529 // Variable to hold the current residual norm. Only used if we're
1530 // not inside the multigrid solver
1531 double norm_res = 0.0;
1532
1533 // Variables to hold the initial residual norm. Only used if we're
1534 // not inside the multigrid solver
1535 double norm_f = 0.0;
1536
1537 // Initialise the value of Iterations
1538 Iterations = 0;
1539
1540 // Calculate the residual only if we're not inside the multigrid solver
1541 if (!Use_as_smoother)
1542 {
1543 // Build the residual vector
1544 current_residual.build(this->distribution_pt(), 0.0);
1545
1546 // Calculate the residual r=b-Ax
1547 matrix_pt->residual(x, rhs, current_residual);
1548
1549 // Calculate the 2-norm of the residual vector
1550 norm_res = current_residual.norm();
1551
1552 // Store the initial norm
1553 norm_f = norm_res;
1554
1555 // If required will document convergence history to screen
1556 // or file (if stream is open)
1558 {
1559 if (!Output_file_stream.is_open())
1560 {
1561 oomph_info << Iterations << " " << norm_res << std::endl;
1562 }
1563 else
1564 {
1565 Output_file_stream << Iterations << " " << norm_res << std::endl;
1566 }
1567 } // if (!Use_as_smoother)
1568 } // if (Doc_convergence_history)
1569
1570 // In each iteration of Gauss-Seidel we need to calculate
1571 // x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
1572 // where L represents the strict lower triangular portion of A,
1573 // D represents the diagonal of A and U represents the strict upper
1574 // triangular portion of A. Since the value of inv(L+D)*rhs remains
1575 // constant throughout we shall calculate it now and store it to
1576 // avoid calculating it again with each iteration.
1577
1578 // Create a temporary matrix pointer to allow for the extraction of the
1579 // row start, column index and value pointers
1580 CRDoubleMatrix* tmp_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1581
1582 // First acquire access to the value, row_start and column_index arrays
1583 // from the compressed row matrix
1584 const double* value_pt = tmp_matrix_pt->value();
1585 const int* row_start_pt = tmp_matrix_pt->row_start();
1586 const int* column_index_pt = tmp_matrix_pt->column_index();
1587
1588 // We've finished using the temporary matrix pointer so make it a null
1589 // pointer
1590 tmp_matrix_pt = 0;
1591
1592 // We can only calculate this constant term if the diagonal entries
1593 // of the matrix are nonzero so we check this first
1594 for (unsigned i = 0; i < n_dof; i++)
1595 {
1596 // Get the index of the last entry below or on the diagonal in row i
1597 unsigned diag_index = Index_of_diagonal_entries[i];
1598
1599 if (unsigned(*(column_index_pt + diag_index)) != i)
1600 {
1601 std::string err_strng = "Gauss-Seidel cannot operate on a matrix with ";
1602 err_strng += "zero diagonal entries.";
1603 throw OomphLibError(
1604 err_strng, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1605 }
1606 } // for (unsigned i=0;i<n_dof;i++)
1607
1608 // If there were no zero diagonal entries we can calculate the vector
1609 // constant_term=inv(L+D)*rhs.
1610
1611 // Create a vector to hold the constant term in the iteration
1612 DoubleVector constant_term(this->distribution_pt(), 0.0);
1613
1614 // Copy the entries of rhs into the vector, constant_term
1615 constant_term = rhs;
1616
1617 // Start by looping over the rows of constant_term
1618 for (unsigned i = 0; i < n_dof; i++)
1619 {
1620 // Get the index of the last entry below or on the diagonal in row i
1621 unsigned diag_index = Index_of_diagonal_entries[i];
1622
1623 // Get the index of the first entry in row i
1624 unsigned row_i_start = *(row_start_pt + i);
1625
1626 // If there are entries strictly below the diagonal then the
1627 // column index of the first entry in the i-th row cannot be
1628 // i (since we have ensured that nonzero entries exist on the
1629 // diagonal so this is the largest the column index of the
1630 // first entry in this row could be)
1631 if (unsigned(*(column_index_pt + row_i_start)) != i)
1632 {
1633 // Initialise the column index variable
1634 unsigned column_index = 0;
1635
1636 // Loop over the entries below the diagonal on row i
1637 for (unsigned j = row_i_start; j < diag_index; j++)
1638 {
1639 // Find the column index of this nonzero entry
1640 column_index = *(column_index_pt + j);
1641
1642 // Subtract the value of A(i,j)*constant_term(j) (where j<i)
1643 constant_term[i] -= (*(value_pt + j)) * constant_term[column_index];
1644 }
1645 } // if (*(column_index_pt+row_i_start)!=i)
1646
1647 // Finish off by calculating constant_term(i)/A(i,i)
1648 constant_term[i] /= *(value_pt + diag_index);
1649 } // for (unsigned i=0;i<n_dof;i++)
1650
1651 // Build a temporary vector to store the value of A*x with each iteration
1652 DoubleVector temp_vec(this->distribution_pt(), 0.0);
1653
1654 // Outermost loop: Run max_iter times (the iteration number)
1655 for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
1656 {
1657 // With each iteration we need to calculate
1658 // x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
1659 // =-inv(L+D)*U*x_old+constant_term.
1660 // Since we've already calculated the vector constant_term, it remains
1661 // for us to calculate inv(L+D)*U*x. To do this we break the calculation
1662 // into two parts:
1663 // (1) The matrix-vector multiplication temp_vec=U*x, and;
1664 // (2) The matrix-vector multiplication inv(L+D)*temp_vec.
1665
1666 //------------------------------
1667 // (1) Calculating temp_vec=U*x:
1668 //------------------------------
1669
1670 // Auxiliary variable to store the index of the first entry on the
1671 // upper triangular portion of the matrix
1672 unsigned upper_tri_start = 0;
1673
1674 // Auxiliary variable to store the index of the first entry in the
1675 // next row
1676 unsigned next_row_start = 0;
1677
1678 // Set the temporary vector temp_vec to be initially be the zero vector
1679 temp_vec.initialise(0.0);
1680
1681 // Loop over the rows of temp_vec (note: we do not loop over the final
1682 // row since all the entries on the last row of a strict upper triangular
1683 // matrix are zero)
1684 for (unsigned i = 0; i < n_dof - 1; i++)
1685 {
1686 // Get the index of the first entry on the upper triangular portion of
1687 // the matrix noting that the i-th entry in Index_of_diagonal_entries
1688 // will store the index of the diagonal entry of the i-th row
1689 upper_tri_start = Index_of_diagonal_entries[i] + 1;
1690
1691 // Get the index of the first entry in the next row
1692 next_row_start = *(row_start_pt + i + 1);
1693
1694 // Variable to store the column index of each entry
1695 unsigned column_index = 0;
1696
1697 // Loop over all of the entries above the diagonal
1698 for (unsigned j = upper_tri_start; j < next_row_start; j++)
1699 {
1700 // Get the column index of this entry
1701 column_index = *(column_index_pt + j);
1702
1703 // Update the value of temp_vec[i] by adding A(i,j)*x(j) (where j>i)
1704 temp_vec[i] += (*(value_pt + j)) * x[column_index];
1705 }
1706 } // for (unsigned i=0;i<n_dof;i++)
1707
1708 //-----------------------------------
1709 // (2) Calculating inv(L+D)*temp_vec:
1710 //-----------------------------------
1711
1712 // Now U*x=temp_vec has been calculated we need to calculate the vector
1713 // y=inv(L+D)*U*x=inv(L+D)*temp_vec.
1714 // Note: the i-th entry in y is given by
1715 // y(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*y(j))/A(i,i).
1716 // We can see from this that to calculate the current entry in y we
1717 // use its previously calculated entries and the current entry in
1718 // temp_vec. As a result, we do not need two separate vectors y and
1719 // temp_vec for this calculation. Instead, we can just update the
1720 // entries in temp_vec (carefully!) using:
1721 // temp_vec(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*temp_vec(j))/A(i,i).
1722
1723 // Start by looping over the rows of rhs
1724 for (unsigned i = 0; i < n_dof; i++)
1725 {
1726 // Get the index of the last entry below or on the diagonal in row i
1727 unsigned diag_index = Index_of_diagonal_entries[i];
1728
1729 // Get the index of the first entry in row i
1730 unsigned row_i_start = *(row_start_pt + i);
1731
1732 // If there are no entries strictly below the diagonal then the
1733 // column index of the first entry in the i-th row will be greater
1734 // than or equal to i
1735 if (unsigned(*(column_index_pt + row_i_start)) != i)
1736 {
1737 // Initialise the column index variable
1738 unsigned column_index = 0;
1739
1740 // Loop over the entries below the diagonal
1741 for (unsigned j = row_i_start; j < diag_index; j++)
1742 {
1743 // Find the column index of this nonzero entry
1744 column_index = *(column_index_pt + j);
1745
1746 // Subtract the value of A(i,j)*y(j) (where j<i)
1747 temp_vec[i] -= (*(value_pt + j)) * temp_vec[column_index];
1748 }
1749 } // if (*(column_index_pt+row_i_start)!=i)
1750
1751 // Finish off by calculating temp_vec(i)/A(i,i)
1752 temp_vec[i] /= *(value_pt + diag_index);
1753 } // for (unsigned i=0;i<n_dof;i++)
1754
1755 //-------------------------
1756 // Updating the solution x:
1757 //-------------------------
1758
1759 // The updated solution is given by
1760 // x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
1761 // =-temp_vec+constant_term,
1762 // We have calculated both constant_term and temp_vec so we simply
1763 // need to calculate x now
1764
1765 // Set x to be the vector, constant_term
1766 x = constant_term;
1767
1768 // Update x to give the next iterate
1769 x -= temp_vec;
1770
1771 // Increment the value of Iterations
1772 Iterations++;
1773
1774 // Calculate the residual only if we're not inside the multigrid solver
1775 if (!Use_as_smoother)
1776 {
1777 // Get residual
1778 matrix_pt->residual(x, rhs, current_residual);
1779
1780 // Calculate the relative residual norm (i.e.
1781 // \frac{\|r_{i}\|}{\|r_{0}\|})
1782 norm_res = current_residual.norm() / norm_f;
1783
1784 // If required will document convergence history to screen or file (if
1785 // stream open)
1787 {
1788 if (!Output_file_stream.is_open())
1789 {
1790 oomph_info << iter_num << " " << norm_res << std::endl;
1791 }
1792 else
1793 {
1794 Output_file_stream << iter_num << " " << norm_res << std::endl;
1795 }
1796 } // if (Doc_convergence_history)
1797
1798 // Check the tolerance only if the residual norm is being computed
1799 if (norm_res < Tolerance)
1800 {
1801 // Break out of the for-loop
1802 break;
1803 }
1804 } // if (!Use_as_smoother)
1805 } // for (unsigned iter_num=0;iter_num<max_iter;iter_num++)
1806
1807 // Calculate the residual only if we're not inside the multigrid solver
1808 if (!Use_as_smoother)
1809 {
1810 // If time documentation is enabled
1811 if (Doc_time)
1812 {
1813 oomph_info << "\nGS converged. Residual norm: " << norm_res
1814 << "\nNumber of iterations to convergence: " << Iterations
1815 << "\n"
1816 << std::endl;
1817 }
1818 } // if (!Use_as_smoother)
1819
1820 // Copy the solution into the solution vector
1821 solution = x;
1822
1823 // Document the time for the solver
1824 double t_end = TimingHelpers::timer();
1825 Solution_time = t_end - t_start;
1826
1827 // If time documentation is enabled
1828 if (Doc_time)
1829 {
1830 oomph_info << "Time for solve with Gauss-Seidel [sec]: " << Solution_time
1831 << std::endl;
1832 }
1833
1834 // If the solver failed to converge and the user asked for an error if
1835 // this happened
1837 {
1838 std::string error_message =
1839 "Solver failed to converge and you requested ";
1840 error_message += "an error on convergence failures.";
1841 throw OomphLibError(
1842 error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1843 }
1844 } // End of solve_helper
1845
1846
1847 /// ///////////////////////////////////////////////////////////////////
1848 /// ///////////////////////////////////////////////////////////////////
1849 /// ///////////////////////////////////////////////////////////////////
1850
1851
1852 //==================================================================
1853 /// Solver: Takes pointer to problem and returns the results
1854 /// vector which contains the solution of the linear system defined
1855 /// by the problem's fully assembled Jacobian and residual vector.
1856 //==================================================================
1857 template<typename MATRIX>
1858 void DampedJacobi<MATRIX>::solve(Problem* const& problem_pt,
1859 DoubleVector& result)
1860 {
1861 // Reset the Use_as_smoother_flag as the solver is not being used
1862 // as a smoother
1863 Use_as_smoother = false;
1864
1865 // Find the # of degrees of freedom (variables)
1866 unsigned n_dof = problem_pt->ndof();
1867
1868 // Initialise timer
1869 double t_start = TimingHelpers::timer();
1870
1871 // We're not re-solving
1872 Resolving = false;
1873
1874 // Get rid of any previously stored data
1876
1877 // Set up the distribution
1878 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
1879
1880 // Assign the distribution to the LinearSolver
1881 this->build_distribution(dist);
1882
1883 // Allocate space for the Jacobian matrix in format specified
1884 // by template parameter
1885 Matrix_pt = new MATRIX;
1886
1887 // Get the nonlinear residuals vector
1888 DoubleVector f;
1889
1890 // Assign the Jacobian and the residuals vector
1891 problem_pt->get_jacobian(f, *Matrix_pt);
1892
1893 // Extract the diagonal entries of the matrix and store them
1894 extract_diagonal_entries(Matrix_pt);
1895
1896 // We've made the matrix, we can delete it...
1897 Matrix_can_be_deleted = true;
1898
1899 // Doc time for setup
1900 double t_end = TimingHelpers::timer();
1901 Jacobian_setup_time = t_end - t_start;
1902
1903 // If time documentation is enabled
1904 if (Doc_time)
1905 {
1906 oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
1907 << std::endl;
1908 }
1909
1910 // Call linear algebra-style solver
1911 solve_helper(Matrix_pt, f, result);
1912
1913 // Kill matrix unless it's still required for resolve
1914 if (!Enable_resolve) clean_up_memory();
1915 } // End of solve
1916
1917 //==================================================================
1918 /// Linear-algebra-type solver: Takes pointer to a matrix and
1919 /// rhs vector and returns the solution of the linear system.
1920 //==================================================================
1921 template<typename MATRIX>
1923 const DoubleVector& rhs,
1924 DoubleVector& solution)
1925 {
1926 // Get number of dofs
1927 unsigned n_dof = rhs.nrow();
1928
1929#ifdef PARANOID
1930 // Upcast the matrix to the appropriate type
1931 MATRIX* tmp_matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
1932
1933 // PARANOID Run the self-tests to check the inputs are correct
1934 this->check_validity_of_solve_helper_inputs<MATRIX>(
1935 tmp_matrix_pt, rhs, solution, n_dof);
1936
1937 // We don't need the pointer any more but we do still need the matrix
1938 // so just make tmp_matrix_pt a null pointer
1939 tmp_matrix_pt = 0;
1940#endif
1941
1942 // Setup the solution if it is not
1943 if (!solution.distribution_pt()->built())
1944 {
1945 // Build the distribution of the solution vector if it hasn't been done
1946 // yet
1947 solution.build(this->distribution_pt(), 0.0);
1948 }
1949 // If the solution has already been set up
1950 else
1951 {
1952 // If we're inside the multigrid solver then as we traverse up the
1953 // hierarchy we use the smoother on the updated approximate solution.
1954 // As such, we should ONLY be resetting all the values to zero if
1955 // we're NOT inside the multigrid solver
1956 if (!Use_as_smoother)
1957 {
1958 // Initialise the vector with all entries set to zero
1959 solution.initialise(0.0);
1960 }
1961 } // if (!solution.distribution_pt()->built())
1962
1963 // Initialise timer
1964 double t_start = TimingHelpers::timer();
1965
1966 // Initial guess isn't necessarily zero (restricted solution from finer
1967 // grids) therefore x needs to be assigned values from the input.
1968 DoubleVector x(solution);
1969
1970 // Create a vector to store the value of the constant term, omega*inv(D)*r
1971 DoubleVector constant_term(this->distribution_pt(), 0.0);
1972
1973 // Calculate the constant term vector
1974 for (unsigned i = 0; i < n_dof; i++)
1975 {
1976 // Assign the i-th entry of constant_term
1977 constant_term[i] = Omega * Matrix_diagonal[i] * rhs[i];
1978 }
1979
1980 // Create a vector to hold the residual. This will only be built if
1981 // we're not inside the multigrid solver
1982 DoubleVector local_residual;
1983
1984 // Variable to store the 2-norm of the residual vector. Only used
1985 // if we are not working inside the MG solver
1986 double norm_res = 0.0;
1987
1988 // Variables to hold the initial residual norm. Only used if we're
1989 // not inside the multigrid solver
1990 double norm_f = 0.0;
1991
1992 // Initialise the value of Iterations
1993 Iterations = 0;
1994
1995 // Calculate the residual only if we're not inside the multigrid solver
1996 if (!Use_as_smoother)
1997 {
1998 // Build the local residual vector
1999 local_residual.build(this->distribution_pt(), 0.0);
2000
2001 // Calculate the residual vector
2002 matrix_pt->residual(x, rhs, local_residual);
2003
2004 // Calculate the 2-norm
2005 norm_res = local_residual.norm();
2006
2007 // Store the initial norm
2008 norm_f = norm_res;
2009
2010 // If required will document convergence history to screen
2011 // or file (if stream is open)
2012 if (Doc_convergence_history)
2013 {
2014 if (!Output_file_stream.is_open())
2015 {
2016 oomph_info << Iterations << " " << norm_res << std::endl;
2017 }
2018 else
2019 {
2020 Output_file_stream << Iterations << " " << norm_res << std::endl;
2021 }
2022 } // if (Doc_convergence_history)
2023 } // if (!Use_as_smoother)
2024
2025 // Create a temporary vector to store the value of A*e
2026 // on each iteration and another to store the residual
2027 // at the end of each iteration
2028 DoubleVector temp_vec(this->distribution_pt(), 0.0);
2029
2030 // Outermost loop: Run up to Max_iter times (the iteration number)
2031 for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
2032 {
2033 // Calculate A*e
2034 matrix_pt->multiply(x, temp_vec);
2035
2036 // Loop over each degree of freedom and update
2037 // the current approximation
2038 for (unsigned idof = 0; idof < n_dof; idof++)
2039 {
2040 // Scale the idof'th entry of temp_vec
2041 // by omega/A(i,i)
2042 temp_vec[idof] *= Omega * Matrix_diagonal[idof];
2043 }
2044
2045 // Update the value of e (in the system Ae=r)
2046 x += constant_term;
2047 x -= temp_vec;
2048
2049 // Increment the value of Iterations
2050 Iterations++;
2051
2052 // Calculate the residual only if we're not inside the multigrid solver
2053 if (!Use_as_smoother)
2054 {
2055 // Get residual
2056 matrix_pt->residual(x, rhs, local_residual);
2057
2058 // Calculate the 2-norm of the residual r=b-Ax
2059 norm_res = local_residual.norm() / norm_f;
2060
2061 // If required, this will document convergence history to
2062 // screen or file (if the stream is open)
2063 if (Doc_convergence_history)
2064 {
2065 if (!Output_file_stream.is_open())
2066 {
2067 oomph_info << Iterations << " " << norm_res << std::endl;
2068 }
2069 else
2070 {
2071 Output_file_stream << Iterations << " " << norm_res << std::endl;
2072 }
2073 } // if (Doc_convergence_history)
2074
2075 // Check the tolerance only if the residual norm is being computed
2076 if (norm_res < Tolerance)
2077 {
2078 // Break out of the for-loop
2079 break;
2080 }
2081 } // if (!Use_as_smoother)
2082 } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
2083
2084 // Calculate the residual only if we're not inside the multigrid solver
2085 if (!Use_as_smoother)
2086 {
2087 // If time documentation is enabled
2088 if (Doc_time)
2089 {
2090 oomph_info << "\nDamped Jacobi converged. Residual norm: " << norm_res
2091 << "\nNumber of iterations to convergence: " << Iterations
2092 << "\n"
2093 << std::endl;
2094 }
2095 } // if (!Use_as_smoother)
2096
2097 // Copy the solution into the solution vector
2098 solution = x;
2099
2100 // Doc. time for solver
2101 double t_end = TimingHelpers::timer();
2102 Solution_time = t_end - t_start;
2103 if (Doc_time)
2104 {
2105 oomph_info << "Time for solve with damped Jacobi [sec]: " << Solution_time
2106 << std::endl;
2107 }
2108
2109 // If the solver failed to converge and the user asked for an error if
2110 // this happened
2111 if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
2112 {
2113 std::string error_message =
2114 "Solver failed to converge and you requested ";
2115 error_message += "an error on convergence failures.";
2116 throw OomphLibError(
2117 error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
2118 }
2119 } // End of solve_helper
2120
2121
2122 /// ////////////////////////////////////////////////////////////////////
2123 /// ////////////////////////////////////////////////////////////////////
2124 /// ////////////////////////////////////////////////////////////////////
2125
2126
2127 //==================================================================
2128 /// \Short Re-solve the system defined by the last assembled Jacobian
2129 /// and the rhs vector specified here. Solution is returned in
2130 /// the vector result.
2131 //==================================================================
2132 template<typename MATRIX>
2134 {
2135 // We are re-solving
2136 Resolving = true;
2137
2138#ifdef PARANOID
2139 if (Matrix_pt == 0)
2140 {
2141 throw OomphLibError("No matrix was stored -- cannot re-solve",
2142 OOMPH_CURRENT_FUNCTION,
2143 OOMPH_EXCEPTION_LOCATION);
2144 }
2145#endif
2146
2147 // Call linear algebra-style solver
2148 this->solve(Matrix_pt, rhs, result);
2149
2150 // Reset re-solving flag
2151 Resolving = false;
2152 }
2153
2154
2155 //==================================================================
2156 /// Solver: Takes pointer to problem and returns the results vector
2157 /// which contains the solution of the linear system defined by
2158 /// the problem's fully assembled Jacobian and residual vector.
2159 //==================================================================
2160 template<typename MATRIX>
2161 void GMRES<MATRIX>::solve(Problem* const& problem_pt, DoubleVector& result)
2162 {
2163 // Find # of degrees of freedom (variables)
2164 unsigned n_dof = problem_pt->ndof();
2165
2166 // Initialise timer
2167 double t_start = TimingHelpers::timer();
2168
2169 // We're not re-solving
2170 Resolving = false;
2171
2172 // Get rid of any previously stored data
2174
2175 // setup the distribution
2176 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
2177 this->build_distribution(dist);
2178
2179 // Get Jacobian matrix in format specified by template parameter
2180 // and nonlinear residual vector
2181 Matrix_pt = new MATRIX;
2182 DoubleVector f;
2183 if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt) != 0)
2184 {
2185 if (dynamic_cast<CRDoubleMatrix*>(Matrix_pt) != 0)
2186 {
2187 dynamic_cast<CRDoubleMatrix*>(Matrix_pt)->build(
2188 this->distribution_pt());
2189 f.build(this->distribution_pt(), 0.0);
2190 }
2191 }
2192 problem_pt->get_jacobian(f, *Matrix_pt);
2193
2194 // We've made the matrix, we can delete it...
2195 Matrix_can_be_deleted = true;
2196
2197 // Doc time for setup
2198 double t_end = TimingHelpers::timer();
2199 Jacobian_setup_time = t_end - t_start;
2200
2201 if (Doc_time)
2202 {
2203 oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
2204 << std::endl;
2205 }
2206
2207 // If we want to compute the gradient for the globally convergent
2208 // Newton method, then do it here
2209 if (Compute_gradient)
2210 {
2211 // Compute it
2212 Matrix_pt->multiply_transpose(f, Gradient_for_glob_conv_newton_solve);
2213 // Set the flag
2214 Gradient_has_been_computed = true;
2215 }
2216
2217 // Call linear algebra-style solver
2218 // If the result distribution is wrong, then redistribute
2219 // before the solve and return to original distribution
2220 // afterwards
2221 if ((result.built()) &&
2222 (!(*result.distribution_pt() == *this->distribution_pt())))
2223 {
2224 LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
2225 result.build(this->distribution_pt(), 0.0);
2226 this->solve_helper(Matrix_pt, f, result);
2227 result.redistribute(&temp_global_dist);
2228 }
2229 // Otherwise just solve
2230 else
2231 {
2232 this->solve_helper(Matrix_pt, f, result);
2233 }
2234
2235 // Kill matrix unless it's still required for resolve
2236 if (!Enable_resolve) clean_up_memory();
2237 };
2238
2239
2240 //==========================================================================
2241 /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
2242 /// and returns the solution of the linear system.
2243 /// based on the algorithm presented in Templates for the
2244 /// Solution of Linear Systems: Building Blocks for Iterative
2245 /// Methods, Barrett,
2246 /// Berry et al, SIAM, 2006 and the implementation in the IML++ library :
2247 /// http://math.nist.gov/iml++/
2248 //==========================================================================
2249 template<typename MATRIX>
2251 const DoubleVector& rhs,
2252 DoubleVector& solution)
2253 {
2254 // Get number of dofs
2255 unsigned n_dof = rhs.nrow();
2256
2257#ifdef PARANOID
2258 // PARANOID check that if the matrix is distributable then it should not be
2259 // then it should not be distributed
2260 if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
2261 {
2262 if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
2263 ->distributed())
2264 {
2265 std::ostringstream error_message_stream;
2266 error_message_stream << "The matrix must not be distributed.";
2267 throw OomphLibError(error_message_stream.str(),
2268 OOMPH_CURRENT_FUNCTION,
2269 OOMPH_EXCEPTION_LOCATION);
2270 }
2271 }
2272 // PARANOID check that this rhs distribution is setup
2273 if (!rhs.built())
2274 {
2275 std::ostringstream error_message_stream;
2276 error_message_stream << "The rhs vector distribution must be setup.";
2277 throw OomphLibError(error_message_stream.str(),
2278 OOMPH_CURRENT_FUNCTION,
2279 OOMPH_EXCEPTION_LOCATION);
2280 }
2281 // PARANOID check that the rhs has the right number of global rows
2282 if (rhs.nrow() != n_dof)
2283 {
2284 throw OomphLibError(
2285 "RHS does not have the same dimension as the linear system",
2286 OOMPH_CURRENT_FUNCTION,
2287 OOMPH_EXCEPTION_LOCATION);
2288 }
2289 // PARANOID check that the rhs is not distributed
2290 if (rhs.distribution_pt()->distributed())
2291 {
2292 std::ostringstream error_message_stream;
2293 error_message_stream << "The rhs vector must not be distributed.";
2294 throw OomphLibError(error_message_stream.str(),
2295 OOMPH_CURRENT_FUNCTION,
2296 OOMPH_EXCEPTION_LOCATION);
2297 }
2298 // PARANOID check that if the result is setup it matches the distribution
2299 // of the rhs
2300 if (solution.built())
2301 {
2302 if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
2303 {
2304 std::ostringstream error_message_stream;
2305 error_message_stream << "If the result distribution is setup then it "
2306 "must be the same as the "
2307 << "rhs distribution";
2308 throw OomphLibError(error_message_stream.str(),
2309 OOMPH_CURRENT_FUNCTION,
2310 OOMPH_EXCEPTION_LOCATION);
2311 }
2312 }
2313#endif
2314
2315 // Reset the time spent applying the preconditioner
2316 Preconditioner_application_time = 0.0;
2317
2318 // Set up the solution if it is not
2319 if (!solution.built())
2320 {
2321 solution.build(this->distribution_pt(), 0.0);
2322 }
2323 // Otherwise initialise to zero
2324 else
2325 {
2326 solution.initialise(0.0);
2327 }
2328
2329 // Time solver
2330 double t_start = TimingHelpers::timer();
2331
2332 // Relative residual
2333 double resid;
2334
2335 // iteration counter
2336 unsigned iter = 1;
2337
2338 // if not using iteration restart set Restart to n_dof
2339 if (!Iteration_restart)
2340 {
2341 Restart = n_dof;
2342 }
2343
2344 // initialise vectors
2345 Vector<double> s(Restart + 1, 0);
2346 Vector<double> cs(Restart + 1);
2347 Vector<double> sn(Restart + 1);
2348 DoubleVector w(this->distribution_pt(), 0.0);
2349
2350 // Setup preconditioner only if we're not re-solving
2351 if (!Resolving)
2352 {
2353 // only setup the preconditioner before solve if require
2354 if (Setup_preconditioner_before_solve)
2355 {
2356 // Setup preconditioner from the Jacobian matrix
2357 double t_start_prec = TimingHelpers::timer();
2358
2359 // do not setup
2360 preconditioner_pt()->setup(matrix_pt);
2361
2362 // Doc time for setup of preconditioner
2363 double t_end_prec = TimingHelpers::timer();
2364 Preconditioner_setup_time = t_end_prec - t_start_prec;
2365
2366 if (Doc_time)
2367 {
2368 oomph_info << "Time for setup of preconditioner [sec]: "
2369 << Preconditioner_setup_time << std::endl;
2370 }
2371 }
2372 }
2373 else
2374 {
2375 if (Doc_time)
2376 {
2377 oomph_info << "Setup of preconditioner is bypassed in resolve mode"
2378 << std::endl;
2379 }
2380 }
2381
2382 // solve b-Jx = Mr for r (assumes x = 0);
2383 DoubleVector r(this->distribution_pt(), 0.0);
2384 if (Preconditioner_LHS)
2385 {
2386 // Start the timer
2387 double t_start_prec = TimingHelpers::timer();
2388
2389 // Apply the preconditioner
2390 preconditioner_pt()->preconditioner_solve(rhs, r);
2391
2392 // Calculate the time taken for the preconditioner solve
2393 Preconditioner_application_time +=
2394 (TimingHelpers::timer() - t_start_prec);
2395 }
2396 else
2397 {
2398 r = rhs;
2399 }
2400 double normb = 0;
2401
2402 // compute norm(r)
2403 double* r_pt = r.values_pt();
2404 for (unsigned i = 0; i < n_dof; i++)
2405 {
2406 normb += r_pt[i] * r_pt[i];
2407 }
2408 normb = sqrt(normb);
2409
2410 // set beta (the initial residual)
2411 double beta = normb;
2412
2413 // compute initial relative residual
2414 if (normb == 0.0) normb = 1;
2415 resid = beta / normb;
2416
2417 // if required will document convergence history to screen or file (if
2418 // stream open)
2419 if (Doc_convergence_history)
2420 {
2421 if (!Output_file_stream.is_open())
2422 {
2423 oomph_info << 0 << " " << resid << std::endl;
2424 }
2425 else
2426 {
2427 Output_file_stream << 0 << " " << resid << std::endl;
2428 }
2429 }
2430
2431 // if GMRES converges immediately
2432 if (resid <= Tolerance)
2433 {
2434 if (Doc_time)
2435 {
2436 oomph_info << "GMRES converged immediately. Normalised residual norm: "
2437 << resid << std::endl;
2438 }
2439 // Doc time for solver
2440 double t_end = TimingHelpers::timer();
2441 Solution_time = t_end - t_start;
2442
2443 if (Doc_time)
2444 {
2445 // Doc the time taken for the preconditioner applications
2446 oomph_info << "Time for all preconditioner applications [sec]: "
2447 << Preconditioner_application_time
2448 << "\n\nTime for solve with GMRES [sec]: " << Solution_time
2449 << std::endl;
2450 }
2451 return;
2452 }
2453
2454
2455 // initialise vector of orthogonal basis vectors (v) and upper hessenberg
2456 // matrix H
2457 // NOTE: for implementation purpose the upper hessenberg matrix indexes are
2458 // are swapped so the matrix is effectively transposed
2460 v.resize(Restart + 1);
2461 Vector<Vector<double>> H(Restart + 1);
2462
2463 // while...
2464 while (iter <= Max_iter)
2465 {
2466 // set zeroth basis vector v[0] to r/beta
2467 v[0].build(this->distribution_pt(), 0.0);
2468 double* v0_pt = v[0].values_pt();
2469 for (unsigned i = 0; i < n_dof; i++)
2470 {
2471 v0_pt[i] = r_pt[i] / beta;
2472 }
2473
2474 //
2475 s[0] = beta;
2476
2477 // inner iteration counter for restarted version
2478 unsigned iter_restart;
2479
2480 // perform iterations
2481 for (iter_restart = 0; iter_restart < Restart && iter <= Max_iter;
2482 iter_restart++, iter++)
2483 {
2484 // resize next column of upper hessenberg matrix
2485 H[iter_restart].resize(iter_restart + 2);
2486
2487 // solve Jv[i] = Mw for w
2488 {
2489 DoubleVector temp(this->distribution_pt(), 0.0);
2490 if (Preconditioner_LHS)
2491 {
2492 // solve Jv[i] = Mw for w
2493 matrix_pt->multiply(v[iter_restart], temp);
2494
2495 // Start the timer
2496 double t_start_prec = TimingHelpers::timer();
2497
2498 // Apply the preconditioner
2499 preconditioner_pt()->preconditioner_solve(temp, w);
2500
2501 // Calculate the time taken for the preconditioner solve
2502 Preconditioner_application_time +=
2503 (TimingHelpers::timer() - t_start_prec);
2504 }
2505 else
2506 {
2507 // Start the timer
2508 double t_start_prec = TimingHelpers::timer();
2509
2510 // w=JM^{-1}v by saad p270
2511 preconditioner_pt()->preconditioner_solve(v[iter_restart], temp);
2512
2513 // Calculate the time taken for the preconditioner solve
2514 Preconditioner_application_time +=
2515 (TimingHelpers::timer() - t_start_prec);
2516
2517 // Do a matrix multiplication
2518 matrix_pt->multiply(temp, w);
2519 }
2520 }
2521
2522 //
2523 double* w_pt = w.values_pt();
2524 for (unsigned k = 0; k <= iter_restart; k++)
2525 {
2526 //
2527 H[iter_restart][k] = 0.0;
2528 double* vk_pt = v[k].values_pt();
2529 for (unsigned i = 0; i < n_dof; i++)
2530 {
2531 H[iter_restart][k] += w_pt[i] * vk_pt[i];
2532 }
2533
2534 //
2535 for (unsigned i = 0; i < n_dof; i++)
2536 {
2537 w_pt[i] -= H[iter_restart][k] * vk_pt[i];
2538 }
2539 }
2540
2541 //
2542 {
2543 double temp_norm_w = 0.0;
2544 for (unsigned i = 0; i < n_dof; i++)
2545 {
2546 temp_norm_w += w_pt[i] * w_pt[i];
2547 }
2548 temp_norm_w = sqrt(temp_norm_w);
2549 H[iter_restart][iter_restart + 1] = temp_norm_w;
2550 }
2551
2552 //
2553 v[iter_restart + 1].build(this->distribution_pt(), 0.0);
2554 double* v_pt = v[iter_restart + 1].values_pt();
2555 for (unsigned i = 0; i < n_dof; i++)
2556 {
2557 v_pt[i] = w_pt[i] / H[iter_restart][iter_restart + 1];
2558 }
2559
2560 //
2561 for (unsigned k = 0; k < iter_restart; k++)
2562 {
2563 apply_plane_rotation(
2564 H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
2565 }
2566 generate_plane_rotation(H[iter_restart][iter_restart],
2567 H[iter_restart][iter_restart + 1],
2568 cs[iter_restart],
2569 sn[iter_restart]);
2570 apply_plane_rotation(H[iter_restart][iter_restart],
2571 H[iter_restart][iter_restart + 1],
2572 cs[iter_restart],
2573 sn[iter_restart]);
2574 apply_plane_rotation(s[iter_restart],
2575 s[iter_restart + 1],
2576 cs[iter_restart],
2577 sn[iter_restart]);
2578
2579 // compute current residual
2580 beta = std::fabs(s[iter_restart + 1]);
2581
2582 // compute relative residual
2583 resid = beta / normb;
2584
2585 // if required will document convergence history to screen or file (if
2586 // stream open)
2587 if (Doc_convergence_history)
2588 {
2589 if (!Output_file_stream.is_open())
2590 {
2591 oomph_info << iter << " " << resid << std::endl;
2592 }
2593 else
2594 {
2595 Output_file_stream << iter << " " << resid << std::endl;
2596 }
2597 }
2598
2599 // if required tolerance found
2600 if (resid < Tolerance)
2601 {
2602 // update result vector
2603 update(iter_restart, H, s, v, solution);
2604
2605 // document convergence
2606 if (Doc_time)
2607 {
2608 oomph_info << std::endl;
2609 oomph_info << "GMRES converged (1). Normalised residual norm: "
2610 << resid << std::endl;
2611 oomph_info << "Number of iterations to convergence: " << iter
2612 << std::endl;
2613 oomph_info << std::endl;
2614 }
2615
2616 // Doc time for solver
2617 double t_end = TimingHelpers::timer();
2618 Solution_time = t_end - t_start;
2619
2620 Iterations = iter;
2621
2622 if (Doc_time)
2623 {
2624 // Doc the time taken for the preconditioner applications
2625 oomph_info << "Time for all preconditioner applications [sec]: "
2626 << Preconditioner_application_time
2627 << "\n\nTime for solve with GMRES [sec]: "
2628 << Solution_time << std::endl;
2629 }
2630 return;
2631 }
2632 }
2633
2634 // update
2635 if (iter_restart > 0) update((iter_restart - 1), H, s, v, solution);
2636
2637 // solve Mr = (b-Jx) for r
2638 {
2639 DoubleVector temp(this->distribution_pt(), 0.0);
2640 matrix_pt->multiply(solution, temp);
2641 double* temp_pt = temp.values_pt();
2642 const double* rhs_pt = rhs.values_pt();
2643 for (unsigned i = 0; i < n_dof; i++)
2644 {
2645 temp_pt[i] = rhs_pt[i] - temp_pt[i];
2646 }
2647
2648 if (Preconditioner_LHS)
2649 {
2650 // Start the timer
2651 double t_start_prec = TimingHelpers::timer();
2652
2653 preconditioner_pt()->preconditioner_solve(temp, r);
2654
2655 // Calculate the time taken for the preconditioner solve
2656 Preconditioner_application_time +=
2657 (TimingHelpers::timer() - t_start_prec);
2658 }
2659 }
2660
2661 // compute current residual
2662 beta = 0.0;
2663 r_pt = r.values_pt();
2664 for (unsigned i = 0; i < n_dof; i++)
2665 {
2666 beta += r_pt[i] * r_pt[i];
2667 }
2668 beta = sqrt(beta);
2669
2670 // if relative residual within tolerance
2671 resid = beta / normb;
2672 if (resid < Tolerance)
2673 {
2674 if (Doc_time)
2675 {
2676 oomph_info << std::endl;
2677 oomph_info << "GMRES converged (2). Normalised residual norm: "
2678 << resid << std::endl;
2679 oomph_info << "Number of iterations to convergence: " << iter
2680 << std::endl;
2681 oomph_info << std::endl;
2682 }
2683
2684 // Doc time for solver
2685 double t_end = TimingHelpers::timer();
2686 Solution_time = t_end - t_start;
2687
2688 Iterations = iter;
2689
2690 if (Doc_time)
2691 {
2692 // Doc the time taken for the preconditioner applications
2693 oomph_info << "Time for all preconditioner applications [sec]: "
2694 << Preconditioner_application_time
2695 << "\n\nTime for solve with GMRES [sec]: "
2696 << Solution_time << std::endl;
2697 }
2698 return;
2699 }
2700 }
2701
2702 // DRAIG: I want to delete this but GMRES doesn't compute the residuals
2703 // properly after this scope when it doesn't converge within the specified
2704 // number of iterations. Not sure why. If you're bored, try to figure
2705 // out why.
2706 {
2707 DoubleVector temp(this->distribution_pt(), 0.0);
2708 matrix_pt->multiply(solution, temp);
2709 temp *= -1.0;
2710 temp += rhs;
2711 resid = temp.norm() / normb;
2712 }
2713
2714 // otherwise GMRES failed convergence
2715 oomph_info << std::endl;
2716 oomph_info << "GMRES did not converge to required tolerance! " << std::endl;
2717 oomph_info << "Returning with normalised residual norm: " << resid
2718 << std::endl;
2719 oomph_info << "after " << Max_iter << " iterations." << std::endl;
2720 oomph_info << std::endl;
2721
2722
2723 if (Throw_error_after_max_iter)
2724 {
2725 std::string err = "Solver failed to converge and you requested an error";
2726 err += " on convergence failures.";
2727 throw OomphLibError(
2728 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
2729 }
2730
2731 return;
2732
2733 } // End GMRES
2734
2735
2736 /// ////////////////////////////////////////////////////////////////////
2737 /// ////////////////////////////////////////////////////////////////////
2738 /// ////////////////////////////////////////////////////////////////////
2739
2740
2741 //==================================================================
2742 /// Solver: Takes pointer to problem and returns the results vector
2743 /// which contains the solution of the linear system defined by
2744 /// the problem's fully assembled Jacobian and residual vector.
2745 //==================================================================
2746 void AugmentedProblemGMRES::solve(Problem* const& problem_pt,
2747 DoubleVector& result)
2748 {
2749 // Find # of degrees of freedom (in the non-augmented problem)
2750 unsigned n_dof = problem_pt->ndof();
2751
2752 // Initialise timer
2753 double t_start = TimingHelpers::timer();
2754
2755 // We're not re-solving
2756 Resolving = false;
2757
2758 // Get rid of any previously stored data
2760
2761 // Create a distribution
2762 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
2763
2764 // ...now build it
2765 this->build_distribution(dist);
2766
2767 // Get the Jacobian matrix and the nonlinear residual vector
2769 DoubleVector f;
2770 if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt) != 0)
2771 {
2772 if (dynamic_cast<CRDoubleMatrix*>(Matrix_pt) != 0)
2773 {
2774 dynamic_cast<CRDoubleMatrix*>(Matrix_pt)->build(
2775 this->distribution_pt());
2776 f.build(this->distribution_pt(), 0.0);
2777 }
2778 }
2779 problem_pt->get_jacobian(f, *Matrix_pt);
2780
2781 // We've made the matrix, we can delete it...
2782 Matrix_can_be_deleted = true;
2783
2784 // Doc time for setup
2785 double t_end = TimingHelpers::timer();
2786 Jacobian_setup_time = t_end - t_start;
2787
2788 if (Doc_time)
2789 {
2790 oomph_info << "Time for setup of Jacobian [sec]: " << Jacobian_setup_time
2791 << std::endl;
2792 }
2793
2794 // Reset the Schur complement scalar entry value
2796
2797 // Call linear algebra-style solver. If the result distribution is wrong,
2798 // then redistribute before the solve and return to original distribution
2799 // afterwards
2800 if ((!(*result.distribution_pt() == *this->distribution_pt())) &&
2801 result.built())
2802 {
2803 // Create a temporary copy of the current result distribution
2804 LinearAlgebraDistribution temp_global_dist(result.distribution_pt());
2805
2806 // Re-build the result vector as an augmented vector
2807 result.build(dist, 0.0);
2808
2809 // Call the auxilliary helper function to do the actual solve
2810 this->solve_helper(Matrix_pt, f, result);
2811
2812 // Re-distribute result vector
2813 result.redistribute(&temp_global_dist);
2814 }
2815 // Otherwise just solve
2816 else
2817 {
2818 // Call the auxilliary helper function to do the actual solve
2819 this->solve_helper(Matrix_pt, f, result);
2820 }
2821
2822 // Kill matrix unless it's still required for resolve
2824 };
2825
2826 //==================================================================
2827 /// \Short Re-solve the system defined by the last assembled Jacobian
2828 /// and the rhs vector specified here. Solution is returned in
2829 /// the vector result.
2830 //==================================================================
2832 DoubleVector& result)
2833 {
2834 // We are re-solving
2835 Resolving = true;
2836
2837#ifdef PARANOID
2838 if (Matrix_pt == 0)
2839 {
2840 throw OomphLibError("No matrix was stored -- cannot re-solve",
2841 OOMPH_CURRENT_FUNCTION,
2842 OOMPH_EXCEPTION_LOCATION);
2843 }
2844#endif
2845
2846 // Call linear algebra-style solver
2847 this->solve(Matrix_pt, rhs, result);
2848
2849 // Reset re-solving flag
2850 Resolving = false;
2851 }
2852
2853
2854 //=============================================================================
2855 /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
2856 /// and returns the solution of the linear system.
2857 /// based on the algorithm presented in Templates for the
2858 /// Solution of Linear Systems: Building Blocks for Iterative Methods,
2859 /// Barrett, Berry et al, SIAM, 2006 and the implementation in the IML++
2860 /// library : http://math.nist.gov/iml++/
2861 //=============================================================================
2863 DoubleMatrixBase* const& input_matrix_pt,
2864 const DoubleVector& rhs,
2865 DoubleVector& lhs)
2866 {
2867 // Get number of dofs (in the non-augmented problem)
2868 unsigned n_dof = input_matrix_pt->nrow();
2869
2870 // How many dofs in the augmented problem?
2871 unsigned n_aug_dof = n_dof + 1;
2872
2873 // Get a handle to this Problems communicator
2875
2876 // Set up the (non-augmented) distribution
2877 LinearAlgebraDistribution dist(comm_pt, n_dof, false);
2878
2879 // Set up the (augmented) distribution
2880 LinearAlgebraDistribution aug_dist(comm_pt, n_aug_dof, false);
2881
2882 // If the input vector doesn't have size N+1 then something has gone wrong
2883 if ((rhs.nrow() != n_dof) || (lhs.nrow() != n_dof))
2884 {
2885 // Create an output stream
2886 std::ostringstream error_message_stream;
2887
2888 // Create the error message
2889 error_message_stream << "RHS vector has " << rhs.nrow()
2890 << " rows and the "
2891 << "LHS vector has " << rhs.nrow()
2892 << " rows but\nthey "
2893 << "should both have " << n_dof << " rows!"
2894 << std::endl;
2895
2896 // Throw an error
2897 throw OomphLibError(error_message_stream.str(),
2898 OOMPH_CURRENT_FUNCTION,
2899 OOMPH_EXCEPTION_LOCATION);
2900 }
2901
2902#ifdef PARANOID
2903 // PARANOID check that if the matrix is distributable then it should not be
2904 // then it should not be distributed
2905 if (dynamic_cast<DistributableLinearAlgebraObject*>(input_matrix_pt) != 0)
2906 {
2907 if (dynamic_cast<DistributableLinearAlgebraObject*>(input_matrix_pt)
2908 ->distributed())
2909 {
2910 std::ostringstream error_message_stream;
2911 error_message_stream << "The matrix must not be distributed.";
2912 throw OomphLibError(error_message_stream.str(),
2913 OOMPH_CURRENT_FUNCTION,
2914 OOMPH_EXCEPTION_LOCATION);
2915 }
2916 }
2917 // PARANOID check that this rhs distribution is setup
2918 if (!rhs.built())
2919 {
2920 std::ostringstream error_message_stream;
2921 error_message_stream << "The rhs vector distribution must be setup.";
2922 throw OomphLibError(error_message_stream.str(),
2923 OOMPH_CURRENT_FUNCTION,
2924 OOMPH_EXCEPTION_LOCATION);
2925 }
2926 // PARANOID check that the rhs is not distributed
2927 if (rhs.distribution_pt()->distributed())
2928 {
2929 std::ostringstream error_message_stream;
2930 error_message_stream << "The rhs vector must not be distributed.";
2931 throw OomphLibError(error_message_stream.str(),
2932 OOMPH_CURRENT_FUNCTION,
2933 OOMPH_EXCEPTION_LOCATION);
2934 }
2935 // PARANOID check that if the result is setup it matches the distribution
2936 // of the rhs
2937 if (lhs.built())
2938 {
2939 if (!(*rhs.distribution_pt() == *lhs.distribution_pt()))
2940 {
2941 std::ostringstream error_message_stream;
2942 error_message_stream << "If the result distribution is setup then it "
2943 "must be the same as the "
2944 << "rhs distribution";
2945 throw OomphLibError(error_message_stream.str(),
2946 OOMPH_CURRENT_FUNCTION,
2947 OOMPH_EXCEPTION_LOCATION);
2948 }
2949 }
2950#endif
2951
2952 // Set up the solution if it is not
2953 if (!lhs.built())
2954 {
2955 // Use the augmented distribution which has size N+1 (not the solvers
2956 // distribution which has size N)
2957 lhs.build(dist, 0.0);
2958 }
2959 // Otherwise initialise to zero
2960 else
2961 {
2962 lhs.initialise(0.0);
2963 }
2964
2965 // Create a vector to store the augmented LHS vector entries
2966 DoubleVector solution(aug_dist, 0.0);
2967
2968 // Time solver
2969 double t_start = TimingHelpers::timer();
2970
2971 // Relative residual
2972 double resid;
2973
2974 // Iteration counter
2975 unsigned iter = 1;
2976
2977 // if not using iteration restart set Restart to n_aug_dof
2978 if (!Iteration_restart)
2979 {
2980 Restart = n_aug_dof;
2981 }
2982
2983 // initialise vectors
2984 Vector<double> s(Restart + 1, 0);
2985 Vector<double> cs(Restart + 1);
2986 Vector<double> sn(Restart + 1);
2987 DoubleVector w(aug_dist, 0.0);
2988
2989 // Set up the preconditioner only if we're not re-solving
2990 if (!Resolving)
2991 {
2992 // Only setup the preconditioner before solve if require
2994 {
2995 // Setup preconditioner from the Jacobian matrix
2996 double t_start_prec = TimingHelpers::timer();
2997
2998 // Do not setup
2999 preconditioner_pt()->setup(input_matrix_pt);
3000
3001 // Doc time for setup of preconditioner
3002 double t_end_prec = TimingHelpers::timer();
3003 Preconditioner_setup_time = t_end_prec - t_start_prec;
3004
3005 // If we're meant to document timings
3006 if (Doc_time)
3007 {
3008 oomph_info << "Time for setup of preconditioner [sec]: "
3009 << Preconditioner_setup_time << std::endl;
3010 }
3011 } // if (Setup_preconditioner_before_solve)
3012 }
3013 else
3014 {
3015 // If we're meant to document timings
3016 if (Doc_time)
3017 {
3018 oomph_info << "Setup of preconditioner is bypassed in resolve mode"
3019 << std::endl;
3020 }
3021 } // if (!Resolving)
3022
3023 // The default system matrix is simply the one that was passed in
3024 CRDoubleMatrix* matrix_pt = dynamic_cast<CRDoubleMatrix*>(input_matrix_pt);
3025
3026 // Storage for the time taken for all preconditioner applications
3027 double t_prec_application_total = 0.0;
3028
3029 //---------------------------------------------------------------
3030 // Compute the scalar component of the diagonal Schur complement
3031 // preconditioner. NOTE: Only compute it if it is currently set
3032 // to the default value of unity. It's quite unlikely that it
3033 // will actually have value one (up to machine precision) in
3034 // practice. This check ensures that we only compute the value
3035 // during the solve function and not during the resolve.
3036 //---------------------------------------------------------------
3037 // If the Schur complement value has its default value of one
3038 if (std::fabs(Schur_complement_scalar - 1.0) < 1.0e-14)
3039 {
3040 // If the bordering vectors haven't actually been properly created
3041 if ((B_pt == 0) || (C_pt == 0))
3042 {
3043 // Throw an error
3044 throw OomphLibError("Bordering vectors have not been created!",
3045 OOMPH_CURRENT_FUNCTION,
3046 OOMPH_EXCEPTION_LOCATION);
3047 }
3048 // Make sure the bordering vectors have the right size
3049 else if ((B_pt->nrow() != n_dof) || (C_pt->nrow() != n_dof))
3050 {
3051 // Throw an error
3052 throw OomphLibError("Bordering vectors do not have the right size!",
3053 OOMPH_CURRENT_FUNCTION,
3054 OOMPH_EXCEPTION_LOCATION);
3055 }
3056
3057 // Allocate storage for A^{-1}b
3058 DoubleVector a_inv_b(this->distribution_pt(), 0.0);
3059
3060 // Compute A^{-1}b with the preconditioner approximation to A^{-1}
3062
3063 // Now compute the scalar component of the Schur complement preconditioner
3064 Schur_complement_scalar = C_pt->dot(a_inv_b);
3065
3066 // If the Schur complement scalar value is pratically zero or is NaN
3067 if ((std::fabs(Schur_complement_scalar) < 1.0e-14) ||
3068 std::isnan(Schur_complement_scalar))
3069 {
3070 // Reset it to have value 2.0 (not 1.0 otherwise we would end up
3071 // recomputing it every resolve) as we "should" be able to use
3072 // any scalar value and end up with the same number of distinct
3073 // eigenvalues for the preconditioned system).
3075 }
3076 } // if (std::fabs(Schur_complement_scalar-1.0)<1.0e-14)
3077
3078 // Storage for the augmented RHS vector
3079 DoubleVector rhs_aug(aug_dist, 0.0);
3080
3081 // Loop over the first n_dof entries of the RHS vector
3082 for (unsigned i = 0; i < n_dof; i++)
3083 {
3084 // Copy the i-th entry over
3085 rhs_aug[i] = rhs[i];
3086 }
3087
3088 // Update the RHS entry associated with the augmented equation
3089 rhs_aug[n_dof] = (*Rhs_pt);
3090
3091 // Storage for b-Jx=Mr for r (assumes x = 0);
3092 DoubleVector r(aug_dist, 0.0);
3093
3094 // If we're using LHS preconditioning
3096 {
3097 // Initialise timer
3098 double t_prec_application_start = TimingHelpers::timer();
3099
3100 // Use the helper function to apply the Schur complement preconditioner
3102
3103 // Doc time for setup of preconditioner
3104 double t_prec_application_end = TimingHelpers::timer();
3105
3106 // Update the preconditioner application time total
3107 t_prec_application_total +=
3108 (t_prec_application_end - t_prec_application_start);
3109 }
3110 // If we're using RHS preconditioning
3111 else
3112 {
3113 // We don't need to do any preconditioning here; just copy the vector over
3114 r = rhs_aug;
3115 } // if (Preconditioner_LHS)
3116
3117 // Compute the norm of the RHS vector
3118 double normb = r.norm();
3119
3120 // Set beta (the initial residual)
3121 double beta = normb;
3122
3123 // If the RHS vector is (numerically equivalent to) the zero vector
3124 if (normb == 0.0)
3125 {
3126 // Initialise the norm of b as unity
3127 normb = 1.0;
3128 }
3129
3130 // Compute the relative residual
3131 resid = beta / normb;
3132
3133 // If required will document convergence history to screen or file (if
3134 // stream open)
3136 {
3137 if (!Output_file_stream.is_open())
3138 {
3139 oomph_info << 0 << " " << resid << std::endl;
3140 }
3141 else
3142 {
3143 Output_file_stream << 0 << " " << resid << std::endl;
3144 }
3145 } // if (Doc_convergence_history)
3146
3147 // If GMRES converges immediately
3148 if (resid <= Tolerance)
3149 {
3150 if (Doc_time)
3151 {
3152 oomph_info << "AugmentedProblemGMRES converged immediately. Normalised "
3153 << "residual norm: " << resid << std::endl;
3154 }
3155
3156 // Doc time for solver
3157 double t_end = TimingHelpers::timer();
3158 Solution_time = t_end - t_start;
3159
3160 if (Doc_time)
3161 {
3163 << "Time for all preconditioner applications [sec]: "
3164 << t_prec_application_total
3165 << "\nTotal time for solve with AugmentedProblemGMRES [sec]: "
3166 << Solution_time << std::endl;
3167 }
3168
3169 // We're done so finish here
3170 return;
3171 }
3172
3173 // initialise vector of orthogonal basis vectors (v) and upper hessenberg
3174 // matrix H
3175 // NOTE: for implementation purpose the upper hessenberg matrix indexes are
3176 // are swapped so the matrix is effectively transposed
3178 v.resize(Restart + 1);
3180
3181 // while...
3182 while (iter <= Max_iter)
3183 {
3184 // Copy the values of r into the zeroth basis vector v[0]
3185 v[0] = r;
3186
3187 // Now scale it so the zeroth basis vector v[0] is r/beta
3188 v[0] /= beta;
3189
3190 // Set the first entry of s
3191 s[0] = beta;
3192
3193 // Counter for the inner iteration counter (restarted version)
3194 unsigned iter_restart = 0;
3195
3196 // Start performing iterations
3197 for (iter_restart = 0; iter_restart < Restart && iter <= Max_iter;
3198 iter_restart++, iter++)
3199 {
3200 // Resize next column of upper hessenberg matrix
3201 H[iter_restart].resize(iter_restart + 2);
3202
3203 // solve Jv[i] = Mw for w
3204 {
3205 // Allocate space for the product J*v[i]
3206 DoubleVector temp(aug_dist, 0.0);
3207
3208 // If we're using LHS preconditioning
3210 {
3211 // Compute J*v[i] where J is now the augmented system matrix
3212 augmented_matrix_multiply(matrix_pt, v[iter_restart], temp);
3213
3214 // Initialise timer
3215 double t_prec_application_start = TimingHelpers::timer();
3216
3217 // Use the helper function to apply the Schur complement
3218 // preconditioner
3220
3221 // End timer
3222 double t_prec_application_end = TimingHelpers::timer();
3223
3224 // Update the preconditioner application time total
3225 t_prec_application_total +=
3226 (t_prec_application_end - t_prec_application_start);
3227 }
3228 else
3229 {
3230 // Initialise timer
3231 double t_prec_application_start = TimingHelpers::timer();
3232
3233 // Compute w=JM^{-1}v by saad p270
3234 apply_schur_complement_preconditioner(v[iter_restart], temp);
3235
3236 // Doc time for setup of preconditioner
3237 double t_prec_application_end = TimingHelpers::timer();
3238
3239 // Update the preconditioner application time total
3240 t_prec_application_total +=
3241 (t_prec_application_end - t_prec_application_start);
3242
3243 // Calculate w=Jv_m where v_m=M^{-1}v
3244 augmented_matrix_multiply(matrix_pt, temp, w);
3245 }
3246 } // Solve Jv[i]=Mw for w
3247
3248 // Get a pointer to the entries of w
3249 double* w_pt = w.values_pt();
3250
3251 // Loop over the rows of the Hessenberg matrix
3252 for (unsigned k = 0; k <= iter_restart; k++)
3253 {
3254 // Compute the (k,iter_restart)-th entry of the Hessenberg matrix
3255 // (which is stored in its transposed form)
3256 H[iter_restart][k] = w.dot(v[k]);
3257
3258 // Get a pointer to the entries of the k-th basis vector
3259 double* vk_pt = v[k].values_pt();
3260
3261 // Loop over the entries of w
3262 for (unsigned i = 0; i < n_aug_dof; i++)
3263 {
3264 // Update the i-th entry of w
3265 w_pt[i] -= H[iter_restart][k] * vk_pt[i];
3266 }
3267 } // for (unsigned k=0;k<=iter_restart;k++)
3268
3269 // Calculate the subdiagonal Hessenberg entry
3270 H[iter_restart][iter_restart + 1] = w.norm();
3271
3272 // Copy the entries of w into the next basis vector
3273 v[iter_restart + 1] = w;
3274
3275 // Now normalise the basis vector with the appropriate Hessenberg entry
3276 v[iter_restart + 1] /= H[iter_restart][iter_restart + 1];
3277
3278 // Loop over the rows of the Hessenberg matrix
3279 for (unsigned k = 0; k < iter_restart; k++)
3280 {
3281 // Apply a Givens rotation to the appropriate entries of H
3283 H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
3284 }
3285
3286 generate_plane_rotation(H[iter_restart][iter_restart],
3287 H[iter_restart][iter_restart + 1],
3288 cs[iter_restart],
3289 sn[iter_restart]);
3290 apply_plane_rotation(H[iter_restart][iter_restart],
3291 H[iter_restart][iter_restart + 1],
3292 cs[iter_restart],
3293 sn[iter_restart]);
3294 apply_plane_rotation(s[iter_restart],
3295 s[iter_restart + 1],
3296 cs[iter_restart],
3297 sn[iter_restart]);
3298
3299 // Compute the current residual
3300 beta = std::fabs(s[iter_restart + 1]);
3301
3302 // compute relative residual
3303 resid = beta / normb;
3304
3305 // DRAIG: Delete...
3306 // oomph_info << "Residual at iteration "
3307 // << iter << ": " << resid << std::endl;
3308
3309 // if required will document convergence history to screen or file (if
3310 // stream open)
3312 {
3313 if (!Output_file_stream.is_open())
3314 {
3315 oomph_info << iter << " " << resid << std::endl;
3316 }
3317 else
3318 {
3319 Output_file_stream << iter << " " << resid << std::endl;
3320 }
3321 } // if (Doc_convergence_history)
3322
3323 // If required tolerance found
3324 if (resid < Tolerance)
3325 {
3326 // update result vector
3327 update(iter_restart, H, s, v, solution);
3328
3329 // Loop over the first n_dof entries of the RHS vector
3330 for (unsigned i = 0; i < n_dof; i++)
3331 {
3332 // Copy the i-th entry over
3333 lhs[i] = solution[i];
3334 }
3335
3336 // Update the bifurcation parameter
3337 (*X_pt) -= solution[n_dof];
3338
3339 // If we're meant to document convergence
3340 if (Doc_time)
3341 {
3342 // Output the convergence to screen
3343 oomph_info << "\nAugmentedProblemGMRES converged (1). Normalised "
3344 << "residual norm: " << resid
3345 << "\nNumber of iterations "
3346 << "to convergence: " << iter << "\n"
3347 << std::endl;
3348 }
3349
3350 // Doc time for solver
3351 double t_end = TimingHelpers::timer();
3352 Solution_time = t_end - t_start;
3353
3354 Iterations = iter;
3355
3356 // If we're meant to document convergence
3357 if (Doc_time)
3358 {
3360 << "Time for all preconditioner applications [sec]: "
3361 << t_prec_application_total
3362 << "\nTotal time for solve with AugmentedProblemGMRES [sec]: "
3363 << Solution_time << std::endl;
3364 }
3365
3366 // We're done now so jump out here
3367 return;
3368 }
3369#ifdef PARANOID
3370 // If the tolerance wasn't met
3371 else
3372 {
3373 // Warn the user if the sub-diagonal Hessenberg entry was almost
3374 // zero...
3375 if (std::fabs(H[iter_restart][iter_restart + 1]) < 1.0e-12)
3376 {
3377 // Create an output stream
3378 std::ostringstream warning_message_stream;
3379
3380 // Create the warning message
3381 warning_message_stream << "Hessenberg subdiagonal entry (="
3382 << H[iter_restart][iter_restart + 1]
3383 << ") almost zero at iteration "
3384 << iter_restart << "\nbut residual "
3385 << "is only " << resid << " which does "
3386 << "not meet the tolerance (" << Tolerance
3387 << ")!" << std::endl;
3388
3389 // Warn the user
3390 OomphLibWarning(warning_message_stream.str(),
3391 OOMPH_CURRENT_FUNCTION,
3392 OOMPH_EXCEPTION_LOCATION);
3393 }
3394 } // if (resid<Tolerance)
3395#endif
3396 } // for (iter_restart=0;iter_restart<Restart&&iter<=Max_iter;...
3397
3398 // If we're not on the first iteration
3399 if (iter_restart > 0)
3400 {
3401 // Do an update
3402 update((iter_restart - 1), H, s, v, solution);
3403 }
3404
3405 // Solve Mr = (b-Jx) for r
3406 {
3407 DoubleVector temp(aug_dist, 0.0);
3408 augmented_matrix_multiply(matrix_pt, solution, temp);
3409 double* temp_pt = temp.values_pt();
3410 const double* rhs_pt = rhs.values_pt();
3411 for (unsigned i = 0; i < n_dof; i++)
3412 {
3413 // Update the i-th entry of temp to be b-Jx
3414 temp_pt[i] = rhs_pt[i] - temp_pt[i];
3415 }
3416
3417 // If we're using LHS preconditioning
3419 {
3420 // Initialise timer
3421 double t_prec_application_start = TimingHelpers::timer();
3422
3423 // Use the helper function to apply the Schur complement
3424 // preconditioner
3426
3427 // Doc time for setup of preconditioner
3428 double t_prec_application_end = TimingHelpers::timer();
3429
3430 // Update the preconditioner application time total
3431 t_prec_application_total +=
3432 (t_prec_application_end - t_prec_application_start);
3433 }
3434 } // Solve Mr = (b-Jx) for r
3435
3436 // Compute the current residual
3437 beta = r.norm();
3438
3439 // if relative residual within tolerance
3440 resid = beta / normb;
3441 if (resid < Tolerance)
3442 {
3443 // If we're meant to document convergence
3444 if (Doc_time)
3445 {
3446 oomph_info << "\nAugmentedProblemGMRES converged (2). Normalised "
3447 << "residual norm: " << resid << "\nNumber of iterations "
3448 << "to convergence: " << iter << "\n"
3449 << std::endl;
3450 }
3451
3452 // Loop over the first n_dof entries of the RHS vector
3453 for (unsigned i = 0; i < n_dof; i++)
3454 {
3455 // Copy the i-th entry over
3456 lhs[i] = solution[i];
3457 }
3458
3459 // Update the bifurcation parameter
3460 (*X_pt) -= solution[n_dof];
3461
3462 // Doc time for solver
3463 double t_end = TimingHelpers::timer();
3464 Solution_time = t_end - t_start;
3465
3466 Iterations = iter;
3467
3468 // If we're meant to document convergence
3469 if (Doc_time)
3470 {
3472 << "Time for all preconditioner applications [sec]: "
3473 << t_prec_application_total
3474 << "\nTotal time for solve with AugmentedProblemGMRES [sec]: "
3475 << Solution_time << std::endl;
3476 }
3477
3478 return;
3479 } // if (resid<Tolerance)
3480 } // while (iter<=Max_iter)
3481
3482 // DRAIG: Delete...
3483 {
3484 DoubleVector temp(aug_dist, 0.0);
3485 augmented_matrix_multiply(matrix_pt, solution, temp);
3486 temp *= -1.0;
3487 temp += rhs;
3488 resid = temp.norm() / normb;
3489 }
3490
3491 // Loop over the first n_dof entries of the RHS vector
3492 for (unsigned i = 0; i < n_dof; i++)
3493 {
3494 // Copy the i-th entry over
3495 lhs[i] = solution[i];
3496 }
3497
3498 // Update the bifurcation parameter
3499 (*X_pt) -= solution[n_dof];
3500
3501 // Otherwise AugmentedProblemGMRES failed convergence
3502 oomph_info << "\nAugmentedProblemGMRES did not converge to required "
3503 << "tolerance!\nReturning with normalised residual norm: "
3504 << resid << "\nafter " << Max_iter << " iterations.\n"
3505 << std::endl;
3506
3508 {
3509 std::string error_message_string = "Solver failed to converge and you ";
3510 error_message_string += "requested an error on convergence failures.";
3511 throw OomphLibError(
3512 error_message_string, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
3513 }
3514
3515 return;
3516 } // End AugmentedProblemGMRES
3517
3518
3519 // Ensure build of required objects
3520 template class BiCGStab<CCDoubleMatrix>;
3521 template class BiCGStab<CRDoubleMatrix>;
3522 template class BiCGStab<DenseDoubleMatrix>;
3523
3524 template class CG<CCDoubleMatrix>;
3525 template class CG<CRDoubleMatrix>;
3526 template class CG<DenseDoubleMatrix>;
3527
3528 template class GS<CCDoubleMatrix>;
3529 template class GS<CRDoubleMatrix>;
3530 template class GS<DenseDoubleMatrix>;
3531
3532 template class DampedJacobi<CCDoubleMatrix>;
3533 template class DampedJacobi<CRDoubleMatrix>;
3534 template class DampedJacobi<DenseDoubleMatrix>;
3535
3536 template class GMRES<CCDoubleMatrix>;
3537 template class GMRES<CRDoubleMatrix>;
3538 template class GMRES<DenseDoubleMatrix>;
3539
3540 // Solvers for SumOfMatrices class
3541 template class BiCGStab<SumOfMatrices>;
3542 template class CG<SumOfMatrices>;
3543 template class GS<SumOfMatrices>;
3544 template class GMRES<SumOfMatrices>;
3545} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
void apply_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Apply plane rotation. This is done using the update:
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
DoubleVector * C_pt
Pointer to the row vector in the bordered system.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
CRDoubleMatrix * Matrix_pt
Pointer to matrix.
bool Iteration_restart
boolean indicating if iteration restarting is used
unsigned Restart
The number of iterations before the iteration proceedure is restarted if iteration restart is used.
void apply_schur_complement_preconditioner(const DoubleVector &rhs, DoubleVector &soln)
Apply the block-diagonal Schur complement preconditioner to compute the LHS which has size N+1 (the a...
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 generate_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Generate a plane rotation. This is done by finding the values of (i....
double Schur_complement_scalar
The scalar component of the Schur complement preconditioner.
void augmented_matrix_multiply(CRDoubleMatrix *matrix_pt, const DoubleVector &x, DoubleVector &soln)
Multiply the vector x by the augmented system matrix.
unsigned Iterations
Number of iterations taken.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
DoubleVector * B_pt
Pointer to the column vector in the bordered system.
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
void update(const unsigned &k, const Vector< Vector< double > > &H, const Vector< double > &s, const Vector< DoubleVector > &v, DoubleVector &x)
Helper function to update the result vector using the result, x=x_0+V_m*y.
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1060
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1072
double * value()
Access to C-style value array.
Definition: matrices.h:1084
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
This is where the actual work is done – different implementations for different matrix types.
void solve(Problem *const &problem_pt, DoubleVector &result)
Use damped Jacobi iteration as an IterativeLinearSolver: This obtains the Jacobian matrix J and the r...
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
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
virtual void multiply(const DoubleVector &x, DoubleVector &soln) const =0
Multiply the matrix by the vector x: soln=Ax.
virtual void residual(const DoubleVector &x, const DoubleVector &b, DoubleVector &residual_)
Find the residual, i.e. r=b-Ax the residual.
Definition: matrices.h:326
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void initialise(const double &v)
initialise the whole vector with value v
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
double norm() const
compute the 2 norm of this vector
bool built() const
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
double * values_pt()
access function to the underlying values
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
unsigned Iterations
Number of iterations taken.
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 clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
MATRIX * Matrix_pt
System matrix pointer in the format specified by the template argument.
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...
double Tolerance
Convergence tolerance.
bool Throw_error_after_max_iter
Should we throw an error instead of just returning when we hit the max iterations?
double Preconditioner_setup_time
Preconditioner setup time.
double Jacobian_setup_time
Jacobian setup time.
static IdentityPreconditioner Default_preconditioner
Default preconditioner: The base class for preconditioners is a fully functional (if trivial!...
unsigned Max_iter
Maximum number of iterations.
double Solution_time
linear solver solution time
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
std::ofstream Output_file_stream
Output file stream for convergence history.
bool Doc_convergence_history
Flag indicating if the convergence history is to be documented.
bool Setup_preconditioner_before_solve
indicates whether the preconditioner should be setup before solve. Default = true;
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
bool built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:77
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:73
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
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....
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1674
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1246
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:3890
void check_validity_of_solve_helper_inputs(MATRIX *const &matrix_pt, const DoubleVector &rhs, 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 solver (or elsewhere) the residual ...
unsigned Max_iter
Max. # of Newton iterations.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void clean_up_memory()
Clean up function that deletes anything dynamically allocated in this namespace.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...