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// Version 1.0; svn revision $LastChangedRevision: 1282 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //The actual solve functions for Iterative solvers.
31 
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35 #include <oomph-lib-config.h>
36 #endif
37 
38 // Oomph-lib includes
40 
41 // Required to force_ get templated builds of iterative solvers for
42 // sumofmatrices class.
43 #include "sum_of_matrices.h"
44 
45 
46 namespace oomph
47 {
48 
49 
50 //==================================================================
51 /// \short Default preconditioner for iterative solvers: The base
52 /// class for preconditioners is a fully functional (if trivial!)
53 /// preconditioner.
54 //==================================================================
55  IdentityPreconditioner IterativeLinearSolver::Default_preconditioner;
56 
57 
58 ///////////////////////////////////////////////////////////////////
59 ///////////////////////////////////////////////////////////////////
60 ///////////////////////////////////////////////////////////////////
61 
62 
63 //==================================================================
64 /// \short Re-solve the system defined by the last assembled Jacobian
65 /// and the rhs vector specified here. Solution is returned in
66 /// the vector result.
67 //==================================================================
68  template<typename MATRIX>
70  DoubleVector &result)
71  {
72 
73  // We are re-solving
74  Resolving=true;
75 
76 #ifdef PARANOID
77  if (Matrix_pt==0)
78  {
79  throw OomphLibError(
80  "No matrix was stored -- cannot re-solve",
81  OOMPH_CURRENT_FUNCTION,
82  OOMPH_EXCEPTION_LOCATION);
83  }
84 #endif
85 
86  // Call linear algebra-style solver
87  this->solve(Matrix_pt,rhs,result);
88 
89  // Reset re-solving flag
90  Resolving=false;
91 
92  }
93 
94 //==================================================================
95 /// Solver: Takes pointer to problem and returns the results vector
96 /// which contains the solution of the linear system defined by
97 /// the problem's fully assembled Jacobian and residual vector.
98 //==================================================================
99  template<typename MATRIX>
100  void BiCGStab<MATRIX>::solve(Problem* const &problem_pt,
101  DoubleVector &result)
102  {
103  // Initialise timer
104 #ifdef OOMPH_HAS_MPI
105  double t_start = TimingHelpers::timer();
106 #else
107  clock_t t_start = clock();
108 #endif
109 
110  // We're not re-solving
111  Resolving=false;
112 
113  // Get rid of any previously stored data
114  clean_up_memory();
115 
116  // Get Jacobian matrix in format specified by template parameter
117  // and nonlinear residual vector
118  Matrix_pt=new MATRIX;
119  DoubleVector f;
120  problem_pt->get_jacobian(f,*Matrix_pt);
121 
122 
123  // We've made the matrix, we can delete it...
124  Matrix_can_be_deleted=true;
125 
126  // Doc time for setup
127 #ifdef OOMPH_HAS_MPI
128  double t_end = TimingHelpers::timer();
129  Jacobian_setup_time= t_end-t_start;
130 #else
131  clock_t t_end = clock();
132  Jacobian_setup_time=double(t_end-t_start)/CLOCKS_PER_SEC;
133 #endif
134 
135  if(Doc_time)
136  {
137  oomph_info << "Time for setup of Jacobian [sec]: "
138  << Jacobian_setup_time << std::endl;
139  }
140 
141  // set the distribution
142  if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
143  {
144  // the solver has the same distribution as the matrix if possible
145  this->build_distribution(dynamic_cast<DistributableLinearAlgebraObject*>
146  (Matrix_pt)->distribution_pt());
147  }
148  else
149  {
150  // the solver has the same distribution as the RHS
151  this->build_distribution(f.distribution_pt());
152  }
153 
154  // Call linear algebra-style solver
155  if((result.built()) &&
156  (!(*result.distribution_pt() == *this->distribution_pt())))
157  {
159  temp_global_dist(result.distribution_pt());
160  result.build(this->distribution_pt(),0.0);
161  this->solve_helper(Matrix_pt,f,result);
162  result.redistribute(&temp_global_dist);
163  }
164  else
165  {
166  this->solve_helper(Matrix_pt,f,result);
167  }
168 
169  // Kill matrix unless it's still required for resolve
170  if (!Enable_resolve) clean_up_memory();
171 
172  };
173 
174 
175 
176 
177 //==================================================================
178 /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
179 /// and returns the solution of the linear system.
180 /// Algorithm and variable names based on "Numerical Linear Algebra
181 /// for High-Performance Computers" by Dongarra, Duff, Sorensen & van
182 /// der Vorst. SIAM (1998), page 185.
183 //==================================================================
184  template<typename MATRIX>
186  const DoubleVector &rhs,
187  DoubleVector &solution)
188  {
189 #ifdef PARANOID
190  // check that the rhs vector is setup
191  if (!rhs.built())
192  {
193  std::ostringstream error_message_stream;
194  error_message_stream
195  << "The vectors rhs must be setup";
196  throw OomphLibError(error_message_stream.str(),
197  OOMPH_CURRENT_FUNCTION,
198  OOMPH_EXCEPTION_LOCATION);
199  }
200 
201  // check that the matrix is square
202  if (matrix_pt->nrow() != matrix_pt->ncol())
203  {
204  std::ostringstream error_message_stream;
205  error_message_stream
206  << "The matrix at matrix_pt must be square.";
207  throw OomphLibError(error_message_stream.str(),
208  OOMPH_CURRENT_FUNCTION,
209  OOMPH_EXCEPTION_LOCATION);
210  }
211 
212  // check that the matrix and the rhs vector have the same nrow()
213  if (matrix_pt->nrow() != rhs.nrow())
214  {
215  std::ostringstream error_message_stream;
216  error_message_stream
217  << "The matrix and the rhs vector must have the same number of rows.";
218  throw OomphLibError(error_message_stream.str(),
219  OOMPH_CURRENT_FUNCTION,
220  OOMPH_EXCEPTION_LOCATION);
221  }
222 
223  // if the matrix is distributable then it too should have the same
224  // communicator as the rhs vector
225  DistributableLinearAlgebraObject* dist_matrix_pt =
226  dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
227  if (dist_matrix_pt != 0)
228  {
229  if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
230  {
231  std::ostringstream error_message_stream;
232  error_message_stream
233  << "The matrix matrix_pt must have the same communicator as the vectors"
234  << " rhs and result must have the same communicator";
235  throw OomphLibError(error_message_stream.str(),
236  OOMPH_CURRENT_FUNCTION,
237  OOMPH_EXCEPTION_LOCATION);
238  }
239  }
240  // if the matrix is not distributable then it the rhs vector should not be
241  // distributed
242  else
243  {
244  if (rhs.distribution_pt()->distributed())
245  {
246  std::ostringstream error_message_stream;
247  error_message_stream
248  << "The matrix (matrix_pt) is not distributable and therefore the rhs"
249  << " vector must not be distributed";
250  throw OomphLibError(error_message_stream.str(),
251  OOMPH_CURRENT_FUNCTION,
252  OOMPH_EXCEPTION_LOCATION);
253  }
254  }
255  // if the result vector is setup then check it has the same distribution
256  // as the rhs
257  if (solution.built())
258  {
259  if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
260  {
261  std::ostringstream error_message_stream;
262  error_message_stream
263  << "The solution vector distribution has been setup; it must have the "
264  << "same distribution as the rhs vector.";
265  throw OomphLibError(error_message_stream.str(),
266  OOMPH_CURRENT_FUNCTION,
267  OOMPH_EXCEPTION_LOCATION);
268  }
269  }
270 #endif
271 
272  // setup the solution if it is not
273  if (!solution.distribution_pt()->built())
274  {
275  solution.build(this->distribution_pt(),0.0);
276  }
277  // zero
278  else
279  {
280  solution.initialise(0.0);
281  }
282 
283  // Get number of dofs
284 // unsigned n_dof=rhs.size();
285  unsigned nrow_local = this->nrow_local();
286 
287  // Time solver
288 #ifdef OOMPH_HAS_MPI
289  double t_start = TimingHelpers::timer();
290 #else
291  clock_t t_start = clock();
292 #endif
293 
294  // Initialise: Zero initial guess so the initial residual is
295  // equal to the RHS, i.e. the nonlinear residual
296  DoubleVector residual(rhs);
297  double residual_norm = residual.norm();
298  double rhs_norm = residual_norm;
299  if (rhs_norm==0.0) rhs_norm=1.0;
300  DoubleVector x(rhs.distribution_pt(),0.0);
301 
302  // Hat residual by copy operation
303  DoubleVector r_hat(residual);
304 
305  // Normalised residual
306  double normalised_residual_norm=residual_norm/rhs_norm;
307 
308  // if required will document convergence history to screen or file (if
309  // stream open)
310  if (Doc_convergence_history)
311  {
312  if (!Output_file_stream.is_open())
313  {
314  oomph_info << 0 << " "
315  << normalised_residual_norm << std::endl;
316  }
317  else
318  {
319  Output_file_stream << 0 << " "
320  << normalised_residual_norm << std::endl;
321  }
322  }
323 
324  // Check immediate convergence
325  if (normalised_residual_norm<Tolerance)
326  {
327  if(Doc_time)
328  {
329  oomph_info << "BiCGStab converged immediately" << std::endl;
330  }
331  solution=x;
332 
333  // Doc time for solver
334  double t_end = TimingHelpers::timer();
335  Solution_time = t_end-t_start;
336 
337  if(Doc_time)
338  {
339  oomph_info << "Time for solve with BiCGStab [sec]: "
340  << Solution_time << std::endl;
341  }
342  return;
343  }
344 
345  // Setup preconditioner only if we're not re-solving
346  if (!Resolving)
347  {
348  // only setup the preconditioner if required
349  if (Setup_preconditioner_before_solve)
350  {
351  //Setup preconditioner from the Jacobian matrix
352  double t_start_prec = TimingHelpers::timer();
353 
354  preconditioner_pt()->setup(matrix_pt);
355 
356  // Doc time for setup of preconditioner
357  double t_end_prec = TimingHelpers::timer();
358  Preconditioner_setup_time = t_end_prec-t_start_prec;
359 
360  if(Doc_time)
361  {
362  oomph_info << "Time for setup of preconditioner [sec]: "
363  << Preconditioner_setup_time << std::endl;
364  }
365  }
366  }
367  else
368  {
369  if(Doc_time)
370  {
371  oomph_info << "Setup of preconditioner is bypassed in resolve mode"
372  << std::endl;
373  }
374  }
375 
376  // Some auxiliary variables
377  double rho_prev=0.0;
378  double alpha=0.0;
379  double omega=0.0;
380  double rho,beta,dot_prod,dot_prod_tt,dot_prod_ts;
381  double s_norm,r_norm;
382 
383  // Some vectors
384  DoubleVector p(this->distribution_pt(),0.0),p_hat(this->distribution_pt(),0.0),
385  v(this->distribution_pt(),0.0),z(this->distribution_pt(),0.0),t(this->distribution_pt(),0.0),
386  s(this->distribution_pt(),0.0);
387 
388  // Loop over max. number of iterations
389  for (unsigned iter=1;iter<=Max_iter;iter++)
390  {
391  // Dot product for rho
392  rho = r_hat.dot(residual);
393 
394  // Breakdown?
395  if (rho==0.0)
396  {
397  oomph_info << "BiCGStab has broken down after " << iter
398  << " iterations" << std::endl;
399  oomph_info << "Returning with current normalised residual of "
400  << normalised_residual_norm << std::endl;
401  }
402 
403  // First step is different
404  if (iter==1)
405  {
406  p=residual;
407  }
408  else
409  {
410  beta=(rho/rho_prev)*(alpha/omega);
411  for (unsigned i=0;i<nrow_local;i++)
412  {
413  p[i] = residual[i] + beta*(p[i]-omega*v[i]);
414  }
415  }
416 
417  // Apply precondtitioner: p_hat=P^-1*p
418  preconditioner_pt()->preconditioner_solve(p,p_hat);
419 
420  // Matrix vector product: v=A*p_hat
421  matrix_pt->multiply(p_hat,v);
422  dot_prod = r_hat.dot(v);
423  alpha=rho/dot_prod;
424  for (unsigned i=0;i<nrow_local;i++)
425  {
426  s[i]=residual[i]-alpha*v[i];
427  }
428  s_norm = s.norm();
429 
430  // Normalised residual
431  normalised_residual_norm=s_norm/rhs_norm;
432 
433  // if required will document convergence history to screen or file (if
434  // stream open)
435  if (Doc_convergence_history)
436  {
437  if (!Output_file_stream.is_open())
438  {
439  oomph_info << double(iter-0.5) << " "
440  << normalised_residual_norm << std::endl;
441  }
442  else
443  {
444  Output_file_stream << double(iter-0.5) << " "
445  << normalised_residual_norm <<std::endl;
446  }
447  }
448 
449  // Converged?
450  if (normalised_residual_norm<Tolerance)
451  {
452  for (unsigned i=0;i<nrow_local;i++)
453  {
454  solution[i]=x[i]+alpha*p_hat[i];
455  }
456 
457  if(Doc_time)
458  {
459  oomph_info << std::endl;
460  oomph_info << "BiCGStab converged. Normalised residual norm: "
461  << normalised_residual_norm << std::endl;
462  oomph_info << "Number of iterations to convergence: "
463  << iter << std::endl;
464  oomph_info << std::endl;
465  }
466 
467  // Store number of iterations taken
468  Iterations = iter;
469 
470  // Doc time for solver
471  double t_end = TimingHelpers::timer();
472  Solution_time = t_end-t_start;
473 
474  if(Doc_time)
475  {
476  oomph_info << "Time for solve with BiCGStab [sec]: "
477  << Solution_time << std::endl;
478  }
479 
480  return;
481  }
482 
483  //Apply precondtitioner: z=P^-1*s
484  preconditioner_pt()->preconditioner_solve(s,z);
485 
486  // Matrix vector product: t=A*z
487  matrix_pt->multiply(z,t);
488  dot_prod_ts=t.dot(s);
489  dot_prod_tt=t.dot(t);
490  omega=dot_prod_ts/dot_prod_tt;
491  for (unsigned i=0;i<nrow_local;i++)
492  {
493  x[i]+=alpha*p_hat[i]+omega*z[i];
494  residual[i]=s[i]-omega*t[i];
495  }
496  r_norm = residual.norm();
497  rho_prev=rho;
498 
499  // Check convergence again
500  normalised_residual_norm=r_norm/rhs_norm;
501 
502  // if required will document convergence history to screen or file (if
503  // stream open)
504  if (Doc_convergence_history)
505  {
506  if (!Output_file_stream.is_open())
507  {
508  oomph_info << iter << " "
509  << normalised_residual_norm << std::endl;
510  }
511  else
512  {
513  Output_file_stream << iter << " "
514  << normalised_residual_norm << std::endl;
515  }
516  }
517 
518 
519  if (normalised_residual_norm<Tolerance)
520  {
521  if(Doc_time)
522  {
523  oomph_info << std::endl;
524  oomph_info << "BiCGStab converged. Normalised residual norm: "
525  << normalised_residual_norm << std::endl;
526  oomph_info << "Number of iterations to convergence: "
527  << iter << std::endl;
528  oomph_info << std::endl;
529  }
530  solution=x;
531 
532  // Store the number of itertions taken.
533  Iterations = iter;
534 
535  // Doc time for solver
536  double t_end = TimingHelpers::timer();
537  Solution_time = t_end-t_start;
538 
539  if(Doc_time)
540  {
541  oomph_info << "Time for solve with BiCGStab [sec]: "
542  << Solution_time << std::endl;
543  }
544  return;
545  }
546 
547 
548  // Breakdown: Omega has to be >0 for to be able to continue
549  if (omega == 0.0)
550  {
551  oomph_info << std::endl;
552  oomph_info << "BiCGStab breakdown with omega=0.0. "
553  << "Normalised residual norm: "
554  << normalised_residual_norm << std::endl;
555  oomph_info << "Number of iterations so far: " << iter << std::endl;
556  oomph_info << std::endl;
557  solution=x;
558 
559  // Store the number of itertions taken.
560  Iterations = iter;
561 
562  // Doc time for solver
563  double t_end = TimingHelpers::timer();
564  Solution_time = t_end-t_start;
565 
566  if(Doc_time)
567  {
568  oomph_info << "Time for solve with BiCGStab [sec]: "
569  << Solution_time << std::endl;
570  }
571  return;
572  }
573 
574 
575 
576  } // end of iteration loop
577 
578 
579  // No convergence
580  oomph_info << std::endl;
581  oomph_info << "BiCGStab did not converge to required tolerance! "
582  << std::endl;
583  oomph_info << "Returning with normalised residual norm: "
584  << normalised_residual_norm << std::endl;
585  oomph_info << "after " << Max_iter<< " iterations." << std::endl;
586  oomph_info << std::endl;
587 
588  solution=x;
589 
590  // Doc time for solver
591  double t_end = TimingHelpers::timer();
592  Solution_time = t_end-t_start;
593 
594  if(Doc_time)
595  {
596  oomph_info << "Time for solve with BiCGStab [sec]: "
597  << Solution_time << std::endl;
598  }
599 
600  if(Throw_error_after_max_iter)
601  {
602  std::string err = "Solver failed to converge and you requested an error";
603  err += " on convergence failures.";
604  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
605  OOMPH_CURRENT_FUNCTION);
606  }
607 
608  }
609 
610 
611 //////////////////////////////////////////////////////////////////////
612 //////////////////////////////////////////////////////////////////////
613 //////////////////////////////////////////////////////////////////////
614 
615 
616 
617 //==================================================================
618 /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
619 /// and returns the solution of the linear system.
620 /// Algorithm and variable names based on "Matrix Computations,
621 /// 2nd Ed." Golub & van Loan, John Hopkins University Press(1989),
622 /// page 529.
623 //==================================================================
624  template<typename MATRIX>
626  const DoubleVector& rhs,
627  DoubleVector& solution)
628  {
629 #ifdef PARANOID
630  // check that the rhs vector is setup
631  if (!rhs.built())
632  {
633  std::ostringstream error_message_stream;
634  error_message_stream
635  << "The vectors rhs must be setup";
636  throw OomphLibError(error_message_stream.str(),
637  OOMPH_CURRENT_FUNCTION,
638  OOMPH_EXCEPTION_LOCATION);
639  }
640 
641  // check that the matrix is square
642  if (matrix_pt->nrow() != matrix_pt->ncol())
643  {
644  std::ostringstream error_message_stream;
645  error_message_stream
646  << "The matrix at matrix_pt must be square.";
647  throw OomphLibError(error_message_stream.str(),
648  OOMPH_CURRENT_FUNCTION,
649  OOMPH_EXCEPTION_LOCATION);
650  }
651 
652  // check that the matrix and the rhs vector have the same nrow()
653  if (matrix_pt->nrow() != rhs.nrow())
654  {
655  std::ostringstream error_message_stream;
656  error_message_stream
657  << "The matrix and the rhs vector must have the same number of rows.";
658  throw OomphLibError(error_message_stream.str(),
659  OOMPH_CURRENT_FUNCTION,
660  OOMPH_EXCEPTION_LOCATION);
661  }
662 
663  // if the matrix is distributable then it too should have the same
664  // communicator as the rhs vector
665  DistributableLinearAlgebraObject* dist_matrix_pt =
666  dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
667  if (dist_matrix_pt != 0)
668  {
669  if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
670  {
671  std::ostringstream error_message_stream;
672  error_message_stream
673  << "The matrix matrix_pt must have the same communicator as the vectors"
674  << " rhs and result must have the same communicator";
675  throw OomphLibError(error_message_stream.str(),
676  OOMPH_CURRENT_FUNCTION,
677  OOMPH_EXCEPTION_LOCATION);
678  }
679  }
680  // if the matrix is not distributable then it the rhs vector should not be
681  // distributed
682  else
683  {
684  if (rhs.distribution_pt()->distributed())
685  {
686  std::ostringstream error_message_stream;
687  error_message_stream
688  << "The matrix (matrix_pt) is not distributable and therefore the rhs"
689  << " vector must not be distributed";
690  throw OomphLibError(error_message_stream.str(),
691  OOMPH_CURRENT_FUNCTION,
692  OOMPH_EXCEPTION_LOCATION);
693  }
694  }
695  // if the result vector is setup then check it has the same distribution
696  // as the rhs
697  if (solution.built())
698  {
699  if (!(*solution.distribution_pt() == *rhs.distribution_pt()))
700  {
701  std::ostringstream error_message_stream;
702  error_message_stream
703  << "The solution vector distribution has been setup; it must have the "
704  << "same distribution as the rhs vector.";
705  throw OomphLibError(error_message_stream.str(),
706  OOMPH_CURRENT_FUNCTION,
707  OOMPH_EXCEPTION_LOCATION);
708  }
709  }
710 #endif
711 
712  // setup the solution if it is not
713  if (!solution.distribution_pt()->built())
714  {
715  solution.build(this->distribution_pt(),0.0);
716  }
717  // zero
718  else
719  {
720  solution.initialise(0.0);
721  }
722 
723  // Get number of dofs
724 // unsigned n_dof=rhs.size();
725  unsigned nrow_local = this->nrow_local();
726 
727  // Initialise counter
728  unsigned counter = 0;
729 
730  // Time solver
731  double t_start = TimingHelpers::timer();
732 
733  // Initialise: Zero initial guess so the initial residual is
734  // equal to the RHS
735  DoubleVector x(this->distribution_pt(),0.0);
736  DoubleVector residual(rhs);
737  double residual_norm = residual.norm();
738  double rhs_norm=residual_norm;
739  if (rhs_norm==0.0) rhs_norm=1.0;
740 
741  // Normalised residual
742  double normalised_residual_norm=residual_norm/rhs_norm;
743 
744  // if required will document convergence history to screen or file (if
745  // stream open)
746  if (Doc_convergence_history)
747  {
748  if (!Output_file_stream.is_open())
749  {
750  oomph_info << 0 << " "
751  << normalised_residual_norm <<std::endl;
752  }
753  else
754  {
755  Output_file_stream << 0 << " "
756  << normalised_residual_norm <<std::endl;
757  }
758  }
759 
760  // Check immediate convergence
761  if (normalised_residual_norm<Tolerance)
762  {
763  if(Doc_time)
764  {
765  oomph_info << "CG converged immediately" << std::endl;
766  }
767  solution=x;
768 
769  // Doc time for solver
770  double t_end = TimingHelpers::timer();
771  Solution_time = t_end-t_start;
772 
773  if(Doc_time)
774  {
775  oomph_info << "Time for solve with CG [sec]: "
776  << Solution_time << std::endl;
777  }
778  return;
779  }
780 
781 
782  // Setup preconditioner only if we're not re-solving
783  if (!Resolving)
784  {
785  // only setup the preconditioner if required
786  if (Setup_preconditioner_before_solve)
787  {
788  //Setup preconditioner from the Jacobian matrix
789  double t_start_prec = TimingHelpers::timer();
790 
791  preconditioner_pt()->setup(matrix_pt);
792 
793  // Doc time for setup of preconditioner
794  double t_end_prec = TimingHelpers::timer();
795  Preconditioner_setup_time = t_end_prec-t_start_prec;
796 
797  if(Doc_time)
798  {
799  oomph_info << "Time for setup of preconditioner [sec]: "
800  << Preconditioner_setup_time << std::endl;
801  }
802  }
803  }
804  else
805  {
806  if(Doc_time)
807  {
808  oomph_info << "Setup of preconditioner is bypassed in resolve mode"
809  << std::endl;
810  }
811  }
812 
813 
814  // Auxiliary vectors
815 // Vector<double> z(n_dof),p(n_dof),jacobian_times_p(n_dof,0.0);
816  DoubleVector z(this->distribution_pt(),0.0), p(this->distribution_pt(),0.0),
817  jacobian_times_p(this->distribution_pt(),0.0);
818 
819  // Auxiliary values
820  double alpha,beta,rz;
821  double prev_rz=0.0;
822 
823  // Main iteration
824  while((normalised_residual_norm>Tolerance)&&(counter!=Max_iter))
825  {
826 
827  //Apply precondtitioner: z=P^-1*r
828  preconditioner_pt()->preconditioner_solve(residual,z);
829 
830  // P vector is computed differently for first and subsequent steps
831  if(counter==0)
832  {
833  p=z;
834  rz=residual.dot(z);
835  }
836  // Subsequent steps
837  else
838  {
839  rz=residual.dot(z);
840  beta=rz/prev_rz;
841  for (unsigned i=0;i<nrow_local;i++)
842  {
843  p[i]=z[i]+beta*p[i];
844  }
845  }
846 
847 
848  // Matrix vector product
849  matrix_pt->multiply(p,jacobian_times_p);
850  double pq = p.dot(jacobian_times_p);
851  alpha=rz/pq;
852 
853  // Update
854  prev_rz=rz;
855  for(unsigned i=0;i<nrow_local;i++)
856  {
857  x[i]+=alpha*p[i];
858  residual[i]-=alpha*jacobian_times_p[i];
859  }
860 
861 
862  //Calculate the 2norm
863  residual_norm = residual.norm();
864 
865  //Difference between the initial and current 2norm residual
866  normalised_residual_norm=residual_norm/rhs_norm;
867 
868 
869  // if required will document convergence history to screen or file (if
870  // stream open)
871  if (Doc_convergence_history)
872  {
873  if (!Output_file_stream.is_open())
874  {
875  oomph_info << counter << " "
876  << normalised_residual_norm << std::endl;
877  }
878  else
879  {
880  Output_file_stream << counter << " "
881  << normalised_residual_norm << std::endl;
882  }
883  }
884 
885  counter=counter+1;
886 
887  }//end while
888 
889 
890  if (counter >= Max_iter)
891  {
892  oomph_info << std::endl;
893  oomph_info << "CG did not converge to required tolerance! " << std::endl;
894  oomph_info << "Returning with normalised residual norm: "
895  << normalised_residual_norm << std::endl;
896  oomph_info << "after " << counter << " iterations." << std::endl;
897  oomph_info << std::endl;
898  }
899  else
900  {
901  if(Doc_time)
902  {
903  oomph_info << std::endl;
904  oomph_info << "CG converged. Normalised residual norm: "
905  << normalised_residual_norm << std::endl;
906  oomph_info << "Number of iterations to convergence: "
907  << counter << std::endl;
908  oomph_info << std::endl;
909  }
910  }
911 
912 
913  // Store number if iterations taken
914  Iterations = counter;
915 
916  // Copy result back
917  solution=x;
918 
919  // Doc time for solver
920  double t_end = TimingHelpers::timer();
921  Solution_time = t_end-t_start;
922 
923  if(Doc_time)
924  {
925  oomph_info << "Time for solve with CG [sec]: "
926  << Solution_time << std::endl;
927  }
928 
929  if((counter >= Max_iter) && (Throw_error_after_max_iter))
930  {
931  std::string err = "Solver failed to converge and you requested an error";
932  err += " on convergence failures.";
933  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
934  OOMPH_CURRENT_FUNCTION);
935  }
936 
937  }//end CG
938 
939 
940 //==================================================================
941 /// \short Re-solve the system defined by the last assembled Jacobian
942 /// and the rhs vector specified here. Solution is returned in
943 /// the vector result.
944 //==================================================================
945  template<typename MATRIX>
947  DoubleVector &result)
948  {
949 
950  // We are re-solving
951  Resolving=true;
952 
953 #ifdef PARANOID
954  if (Matrix_pt==0)
955  {
956  throw OomphLibError(
957  "No matrix was stored -- cannot re-solve",
958  OOMPH_CURRENT_FUNCTION,
959  OOMPH_EXCEPTION_LOCATION);
960  }
961 #endif
962 
963  // Call linear algebra-style solver
964  this->solve(Matrix_pt,rhs,result);
965 
966  // Reset re-solving flag
967  Resolving=false;
968 
969  }
970 
971 
972 //==================================================================
973 /// Solver: Takes pointer to problem and returns the results vector
974 /// which contains the solution of the linear system defined by
975 /// the problem's fully assembled Jacobian and residual vector.
976 //==================================================================
977  template<typename MATRIX>
978  void CG<MATRIX>::solve(Problem* const &problem_pt, DoubleVector &result)
979  {
980  // Initialise timer
981  double t_start = TimingHelpers::timer();
982 
983  // We're not re-solving
984  Resolving=false;
985 
986  // Get rid of any previously stored data
987  clean_up_memory();
988 
989  // Get Jacobian matrix in format specified by template parameter
990  // and nonlinear residual vector
991  Matrix_pt=new MATRIX;
992  DoubleVector f;
993  problem_pt->get_jacobian(f,*Matrix_pt);
994 
995  // We've made the matrix, we can delete it...
996  Matrix_can_be_deleted=true;
997 
998  // Doc time for setup
999  double t_end = TimingHelpers::timer();
1000  Jacobian_setup_time= t_end-t_start;
1001 
1002  if(Doc_time)
1003  {
1004  oomph_info << "Time for setup of Jacobian [sec]: "
1005  << Jacobian_setup_time << std::endl;
1006  }
1007 
1008  // set the distribution
1009  if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt))
1010  {
1011  // the solver has the same distribution as the matrix if possible
1012  this->build_distribution(dynamic_cast<DistributableLinearAlgebraObject*>
1013  (Matrix_pt)->distribution_pt());
1014  }
1015  else
1016  {
1017  // the solver has the same distribution as the RHS
1018  this->build_distribution(f.distribution_pt());
1019  }
1020 
1021  // if the result vector is not setup
1022  if (!result.distribution_pt()->built())
1023  {
1024  result.build(this->distribution_pt(),0.0);
1025  }
1026 
1027  // Call linear algebra-style solver
1028  if (!(*result.distribution_pt() == *this->distribution_pt()))
1029  {
1031  temp_global_dist(result.distribution_pt());
1032  result.build(this->distribution_pt(),0.0);
1033  this->solve_helper(Matrix_pt,f,result);
1034  result.redistribute(&temp_global_dist);
1035  }
1036  else
1037  {
1038  this->solve_helper(Matrix_pt,f,result);
1039  }
1040 
1041  // Kill matrix unless it's still required for resolve
1042  if (!Enable_resolve) clean_up_memory();
1043 
1044  };
1045 
1046 
1047 //////////////////////////////////////////////////////////////////////
1048 //////////////////////////////////////////////////////////////////////
1049 //////////////////////////////////////////////////////////////////////
1050 
1051 
1052 //==================================================================
1053 /// \short Self-test to be called inside solve_helper to ensure
1054 /// that all inputs are consistent and everything that needs to
1055 /// be built, is.
1056 //==================================================================
1057  template<typename MATRIX>
1059  const DoubleVector& rhs,
1060  DoubleVector& solution,
1061  const double& n_dof)
1062  {
1063  // Check that if the matrix is not distributable then it should not
1064  // be distributed
1065  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
1066  {
1067  if (dynamic_cast<DistributableLinearAlgebraObject*>
1068  (matrix_pt)->distributed())
1069  {
1070  std::ostringstream error_message_stream;
1071  error_message_stream << "The matrix must not be distributed.";
1072  throw OomphLibError(error_message_stream.str(),
1073  OOMPH_CURRENT_FUNCTION,
1074  OOMPH_EXCEPTION_LOCATION);
1075  }
1076  } // if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)!=0)
1077  // Check that this rhs distribution is setup
1078  if (!rhs.built())
1079  {
1080  std::ostringstream error_message_stream;
1081  error_message_stream << "The rhs vector distribution must be setup.";
1082  throw OomphLibError(error_message_stream.str(),
1083  OOMPH_CURRENT_FUNCTION,
1084  OOMPH_EXCEPTION_LOCATION);
1085  }
1086  // Check that the rhs has the right number of global rows
1087  if(rhs.nrow()!=n_dof)
1088  {
1089  std::ostringstream error_message_stream;
1090  error_message_stream << "RHS does not have the same dimension as "
1091  << "the linear system";
1092  throw OomphLibError(error_message_stream.str(),
1093  OOMPH_CURRENT_FUNCTION,
1094  OOMPH_EXCEPTION_LOCATION);
1095  }
1096  // Check that the rhs is not distributed
1097  if (rhs.distribution_pt()->distributed())
1098  {
1099  std::ostringstream error_message_stream;
1100  error_message_stream << "The rhs vector must not be distributed.";
1101  throw OomphLibError(error_message_stream.str(),
1102  OOMPH_CURRENT_FUNCTION,
1103  OOMPH_EXCEPTION_LOCATION);
1104  }
1105  // Check that if the result is setup it matches the distribution
1106  // of the rhs
1107  if (solution.built())
1108  {
1109  if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
1110  {
1111  std::ostringstream error_message_stream;
1112  error_message_stream << "If the result distribution is setup then it "
1113  << "must be the same as the rhs distribution";
1114  throw OomphLibError(error_message_stream.str(),
1115  OOMPH_CURRENT_FUNCTION,
1116  OOMPH_EXCEPTION_LOCATION);
1117  }
1118  } // if (solution.built())
1119  } // End of check_validity_of_solve_helper_inputs
1120 
1121 
1122 ///////////////////////////////////////////////////////////////////////
1123 ///////////////////////////////////////////////////////////////////////
1124 ///////////////////////////////////////////////////////////////////////
1125 
1126 
1127 //==================================================================
1128 /// \short Solver: Takes pointer to problem and returns the results
1129 /// vector which contains the solution of the linear system defined
1130 /// by the problem's fully assembled Jacobian and residual vector.
1131 //==================================================================
1132  template<typename MATRIX>
1133  void GS<MATRIX>::solve(Problem* const &problem_pt,DoubleVector &result)
1134  {
1135  // Reset the Use_as_smoother_flag as the solver is not being used
1136  // as a smoother
1137  Use_as_smoother=false;
1138 
1139  // Find the # of degrees of freedom (variables)
1140  unsigned n_dof=problem_pt->ndof();
1141 
1142  // Initialise timer
1143  double t_start=TimingHelpers::timer();
1144 
1145  // We're not re-solving
1146  Resolving=false;
1147 
1148  // Get rid of any previously stored data
1149  clean_up_memory();
1150 
1151  // Set up the distribution
1152  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),
1153  n_dof,false);
1154 
1155  // Assign the distribution to the LinearSolver
1156  this->build_distribution(dist);
1157 
1158  // Allocate space for the Jacobian matrix in format specified
1159  // by template parameter
1160  Matrix_pt=new MATRIX;
1161 
1162  // Get the nonlinear residuals vector
1163  DoubleVector f;
1164 
1165  // Assign the Jacobian and the residuals vector
1166  problem_pt->get_jacobian(f,*Matrix_pt);
1167 
1168  // We've made the matrix, we can delete it...
1169  Matrix_can_be_deleted=true;
1170 
1171  // Doc time for setup
1172  double t_end=TimingHelpers::timer();
1173  Jacobian_setup_time=t_end-t_start;
1174 
1175  // If time documentation is enabled
1176  if(Doc_time)
1177  {
1178  oomph_info << "Time for setup of Jacobian [sec]: "
1179  << Jacobian_setup_time << std::endl;
1180  }
1181 
1182  // Call linear algebra-style solver
1183  this->solve_helper(Matrix_pt,f,result);
1184 
1185  // Kill matrix unless it's still required for resolve
1186  if (!Enable_resolve) clean_up_memory();
1187  } // End of solve
1188 
1189 //==================================================================
1190 /// \short Linear-algebra-type solver: Takes pointer to a matrix and
1191 /// rhs vector and returns the solution of the linear system.
1192 //==================================================================
1193  template<typename MATRIX>
1195  const DoubleVector& rhs,
1196  DoubleVector& solution)
1197  {
1198  // Get number of dofs
1199  unsigned n_dof=rhs.nrow();
1200 
1201 #ifdef PARANOID
1202  // Upcast the matrix to the appropriate type
1203  MATRIX* tmp_matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
1204 
1205  // PARANOID Run the self-tests to check the inputs are correct
1206  this->check_validity_of_solve_helper_inputs<MATRIX>(tmp_matrix_pt,rhs,
1207  solution,n_dof);
1208 
1209  // We don't need the pointer any more but we do still need the matrix
1210  // so just make tmp_matrix_pt a null pointer
1211  tmp_matrix_pt=0;
1212 #endif
1213 
1214  // Set up the solution if it is not
1215  if (!solution.distribution_pt()->built())
1216  {
1217  // Build it!
1218  solution.build(this->distribution_pt(),0.0);
1219  }
1220  // If the solution is already set up
1221  else
1222  {
1223  // If we're inside the multigrid solver then as we traverse up the
1224  // hierarchy we use the smoother on the updated approximate solution.
1225  // As such, we should ONLY be resetting all the values to zero if
1226  // we're NOT inside the multigrid solver
1227  if (!Use_as_smoother)
1228  {
1229  // Initialise the vector with all entries set to zero
1230  solution.initialise(0.0);
1231  }
1232  } // if (!solution.distribution_pt()->built())
1233 
1234  // Initialise timer
1235  double t_start=TimingHelpers::timer();
1236 
1237  // Copy the solution vector into x
1238  DoubleVector x(solution);
1239 
1240  // Create a vector to hold the residual. This will only be built if
1241  // we're not inside the multigrid solver
1242  DoubleVector current_residual;
1243 
1244  // Variable to hold the current residual norm. Only used if we're
1245  // not inside the multigrid solver
1246  double norm_res=0.0;
1247 
1248  // Variables to hold the initial residual norm. Only used if we're
1249  // not inside the multigrid solver
1250  double norm_f=0.0;
1251 
1252  // Initialise the value of Iterations
1253  Iterations=0;
1254 
1255  // Calculate the residual only if we're not inside the multigrid solver
1256  if (!Use_as_smoother)
1257  {
1258  // Build the residual vector
1259  current_residual.build(this->distribution_pt(),0.0);
1260 
1261  // Calculate the residual r=b-Ax
1262  matrix_pt->residual(x,rhs,current_residual);
1263 
1264  // Calculate the 2-norm of the residual vector
1265  norm_res=current_residual.norm();
1266 
1267  // Store the initial norm
1268  norm_f=norm_res;
1269 
1270  // If required, document the convergence history to screen or file (if
1271  // the output stream open)
1272  if (Doc_convergence_history)
1273  {
1274  // If the output file stream isn't open
1275  if (!Output_file_stream.is_open())
1276  {
1277  // Output the result to screen
1278  oomph_info << Iterations << " " << norm_res << std::endl;
1279  }
1280  // If the output file stream is open
1281  else
1282  {
1283  // Document the result to file
1284  Output_file_stream << Iterations << " " << norm_res <<std::endl;
1285  }
1286  } // if (Doc_convergence_history)
1287  } // if (!Use_as_smoother)
1288 
1289  // Start of the main GS loop: run up to Max_iter times
1290  for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
1291  {
1292  // Loop over rows
1293  for(unsigned i=0;i<n_dof;i++)
1294  {
1295  double dummy=rhs[i];
1296  for (unsigned j=0;j<i;j++)
1297  {
1298  dummy-=(*(matrix_pt))(i,j)*x[j];
1299  }
1300  for (unsigned j=(i+1);j<n_dof;j++)
1301  {
1302  dummy-=(*(matrix_pt))(i,j)*x[j];
1303  }
1304  x[i]=dummy/(*(matrix_pt))(i,i);
1305  } // for(unsigned i=0;i<n_dof;i++)
1306 
1307  // Increment the value of Iterations
1308  Iterations++;
1309 
1310  // Calculate the residual only if we're not inside the multigrid solver
1311  if (!Use_as_smoother)
1312  {
1313  // Get residual
1314  matrix_pt->residual(x,rhs,current_residual);
1315 
1316  // Calculate the relative residual norm (i.e. \frac{\|r_{i}\|}{\|r_{0}\|})
1317  norm_res=current_residual.norm()/norm_f;
1318 
1319  // If required will document convergence history to screen or file (if
1320  // stream open)
1321  if (Doc_convergence_history)
1322  {
1323  if (!Output_file_stream.is_open())
1324  {
1325  oomph_info << Iterations << " " << norm_res << std::endl;
1326  }
1327  else
1328  {
1329  Output_file_stream << Iterations << " " << norm_res << std::endl;
1330  }
1331  } // if (Doc_convergence_history)
1332 
1333  // Check the tolerance only if the residual norm is being computed
1334  if (norm_res<Tolerance)
1335  {
1336  // Break out of the for-loop
1337  break;
1338  }
1339  } // if (!Use_as_smoother)
1340  } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
1341 
1342  // Calculate the residual only if we're not inside the multigrid solver
1343  if (!Use_as_smoother)
1344  {
1345  // If time documentation is enabled
1346  if(Doc_time)
1347  {
1348  oomph_info << "\nGS converged. Residual norm: " << norm_res
1349  << "\nNumber of iterations to convergence: " << Iterations
1350  << "\n" << std::endl;
1351  }
1352  } // if (!Use_as_smoother)
1353 
1354  // Copy result into result
1355  solution=x;
1356 
1357  // Doc. time for solver
1358  double t_end = TimingHelpers::timer();
1359  Solution_time = t_end-t_start;
1360  if(Doc_time)
1361  {
1362  oomph_info << "Time for solve with GS [sec]: "
1363  << Solution_time << std::endl;
1364  }
1365 
1366  // If the solver failed to converge and the user asked for an error if
1367  // this happened
1368  if((Iterations>Max_iter-1)&&(Throw_error_after_max_iter))
1369  {
1370  std::string error_message="Solver failed to converge and you requested ";
1371  error_message+="an error on convergence failures.";
1372  throw OomphLibError(error_message,
1373  OOMPH_EXCEPTION_LOCATION,
1374  OOMPH_CURRENT_FUNCTION);
1375  }
1376  } // End of solve_helper
1377 
1378 
1379 ///////////////////////////////////////////////////////////////////////
1380 ///////////////////////////////////////////////////////////////////////
1381 ///////////////////////////////////////////////////////////////////////
1382 
1383 
1384 //==================================================================
1385 /// \short Explicit template specialisation of the solver for CR
1386 /// matrices: Takes pointer to problem and returns the results
1387 /// vector which contains the solution of the linear system defined
1388 /// by the problem's fully assembled Jacobian and residual vector.
1389 //==================================================================
1390  void GS<CRDoubleMatrix>::solve(Problem* const &problem_pt,
1391  DoubleVector& result)
1392  {
1393  // Reset the Use_as_smoother_flag as the solver is not being used
1394  // as a smoother
1395  Use_as_smoother=false;
1396 
1397  // Find # of degrees of freedom (variables)
1398  unsigned n_dof=problem_pt->ndof();
1399 
1400  // Initialise timer
1401  double t_start=TimingHelpers::timer();
1402 
1403  // We're not re-solving
1404  Resolving=false;
1405 
1406  // Get rid of any previously stored data
1407  clean_up_memory();
1408 
1409  // Set up the distribution
1410  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n_dof,false);
1411 
1412  // Build the distribution of the LinearSolver
1413  this->build_distribution(dist);
1414 
1415  // Allocate space for the Jacobian matrix in CRDoubleMatrix format
1417 
1418  // Build the matrix
1419  Matrix_pt->build(this->distribution_pt());
1420 
1421  // Allocate space for the residuals vector
1422  DoubleVector f;
1423 
1424  // Build it
1425  f.build(this->distribution_pt(),0.0);
1426 
1427  // Get the global Jacobian and the accompanying residuals vector
1428  problem_pt->get_jacobian(f,*Matrix_pt);
1429 
1430  // Doc. time for setup
1431  double t_end=TimingHelpers::timer();
1432  Jacobian_setup_time=t_end-t_start;
1433 
1434  // Run the generic setup function
1435  setup_helper(Matrix_pt);
1436 
1437  // We've made the matrix, we can delete it...
1438  Matrix_can_be_deleted=true;
1439 
1440  // If time documentation is enabled
1441  if(Doc_time)
1442  {
1443  // Output the time for the assembly of the Jacobian to screen
1444  oomph_info << "Time for setup of Jacobian [sec]: "
1445  << Jacobian_setup_time << std::endl;
1446  }
1447 
1448  // Call linear algebra-style solver
1449  solve_helper(Matrix_pt,f,result);
1450 
1451  // Kill matrix unless it's still required for resolve
1453  } // End of solve
1454 
1455 //==================================================================
1456 /// \short Explicit template specialisation of the smoother_setup
1457 /// function for CR matrices. Set up the smoother for the matrix
1458 /// specified by the pointer. This definition of the smoother_setup
1459 /// has the added functionality that it sorts the entries in the
1460 /// given matrix so that the CRDoubleMatrix implementation of the
1461 /// solver can be used.
1462 //==================================================================
1464  {
1465  // Assume the matrix has been passed in from the outside so we must
1466  // not delete it. This is needed to avoid pre- and post-smoothers
1467  // deleting the same matrix in the MG solver. If this was originally
1468  // set to TRUE then this will be sorted out in the other functions
1469  // from which this was called
1470  Matrix_can_be_deleted=false;
1471 
1472  // Upcast the input matrix to system matrix to the type MATRIX
1473  Matrix_pt=dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1474 
1475  // The system matrix here is a CRDoubleMatrix. To make use of the
1476  // specific implementation of the solver for this type of matrix we
1477  // need to make sure the entries are arranged correctly
1478  Matrix_pt->sort_entries();
1479 
1480  // Now get access to the vector Index_of_diagonal_entries
1481  Index_of_diagonal_entries=Matrix_pt->get_index_of_diagonal_entries();
1482 
1483 #ifdef PARANOID
1484  // Create a boolean variable to store the result of whether or not
1485  // the entries are in the correct order
1486  bool is_matrix_sorted=true;
1487 
1488  // Check to make sure the entries have been sorted properly
1489  is_matrix_sorted=Matrix_pt->entries_are_sorted();
1490 
1491  // If the entries are not sorted properly
1492  if (!is_matrix_sorted)
1493  {
1494  // Throw an error; the solver won't work unless the matrix is sorted
1495  // properly
1496  throw OomphLibError("Matrix is not sorted correctly",
1497  OOMPH_CURRENT_FUNCTION,
1498  OOMPH_EXCEPTION_LOCATION);
1499  }
1500 #endif
1501  } // End of setup_helper
1502 
1503 //==================================================================
1504 /// \short Explicit template specialisation of the solve_helper
1505 /// function for CR matrices. Exploiting the sparsity of the given
1506 /// matrix allows for a much faster iterative solver.
1507 //==================================================================
1509  const DoubleVector &rhs,
1510  DoubleVector &solution)
1511  {
1512  // Get number of dofs
1513  unsigned n_dof=rhs.nrow();
1514 
1515 #ifdef PARANOID
1516  // Upcast the matrix to the appropriate type
1517  CRDoubleMatrix* cr_matrix_pt=dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1518 
1519  // PARANOID Run the self-tests to check the inputs are correct
1520  this->check_validity_of_solve_helper_inputs<CRDoubleMatrix>(cr_matrix_pt,rhs,
1521  solution,n_dof);
1522 
1523  // We don't need the pointer any more but we do still need the matrix
1524  // so just make tmp_matrix_pt a null pointer
1525  cr_matrix_pt=0;
1526 #endif
1527 
1528  // Setup the solution if it is not
1529  if (!solution.distribution_pt()->built())
1530  {
1531  solution.build(this->distribution_pt(),0.0);
1532  }
1533  // If the solution is already set up
1534  else
1535  {
1536  // If we're inside the multigrid solver then as we traverse up the
1537  // hierarchy we use the smoother on the updated approximate solution.
1538  // As such, we should ONLY be resetting all the values to zero if
1539  // we're NOT inside the multigrid solver
1540  if (!Use_as_smoother)
1541  {
1542  // Initialise the vector with all entries set to zero
1543  solution.initialise(0.0);
1544  }
1545  } // if (!solution.distribution_pt()->built())
1546 
1547  // Initialise timer
1548  double t_start=TimingHelpers::timer();
1549 
1550  // Copy the solution vector into x
1551  DoubleVector x(solution);
1552 
1553  // Create a vector to hold the residual. This will only be built if
1554  // we're not inside the multigrid solver
1555  DoubleVector current_residual;
1556 
1557  // Variable to hold the current residual norm. Only used if we're
1558  // not inside the multigrid solver
1559  double norm_res=0.0;
1560 
1561  // Variables to hold the initial residual norm. Only used if we're
1562  // not inside the multigrid solver
1563  double norm_f=0.0;
1564 
1565  // Initialise the value of Iterations
1566  Iterations=0;
1567 
1568  // Calculate the residual only if we're not inside the multigrid solver
1569  if (!Use_as_smoother)
1570  {
1571  // Build the residual vector
1572  current_residual.build(this->distribution_pt(),0.0);
1573 
1574  // Calculate the residual r=b-Ax
1575  matrix_pt->residual(x,rhs,current_residual);
1576 
1577  // Calculate the 2-norm of the residual vector
1578  norm_res=current_residual.norm();
1579 
1580  // Store the initial norm
1581  norm_f=norm_res;
1582 
1583  // If required will document convergence history to screen
1584  // or file (if stream is open)
1586  {
1587  if (!Output_file_stream.is_open())
1588  {
1589  oomph_info << Iterations << " " << norm_res << std::endl;
1590  }
1591  else
1592  {
1593  Output_file_stream << Iterations << " " << norm_res << std::endl;
1594  }
1595  } // if (!Use_as_smoother)
1596  } // if (Doc_convergence_history)
1597 
1598  // In each iteration of Gauss-Seidel we need to calculate
1599  // x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
1600  // where L represents the strict lower triangular portion of A,
1601  // D represents the diagonal of A and U represents the strict upper
1602  // triangular portion of A. Since the value of inv(L+D)*rhs remains
1603  // constant throughout we shall calculate it now and store it to
1604  // avoid calculating it again with each iteration.
1605 
1606  // Create a temporary matrix pointer to allow for the extraction of the
1607  // row start, column index and value pointers
1608  CRDoubleMatrix* tmp_matrix_pt=dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1609 
1610  // First acquire access to the value, row_start and column_index arrays
1611  // from the compressed row matrix
1612  const double* value_pt=tmp_matrix_pt->value();
1613  const int* row_start_pt=tmp_matrix_pt->row_start();
1614  const int* column_index_pt=tmp_matrix_pt->column_index();
1615 
1616  // We've finished using the temporary matrix pointer so make it a null pointer
1617  tmp_matrix_pt=0;
1618 
1619  // We can only calculate this constant term if the diagonal entries
1620  // of the matrix are nonzero so we check this first
1621  for (unsigned i=0;i<n_dof;i++)
1622  {
1623  // Get the index of the last entry below or on the diagonal in row i
1624  unsigned diag_index=Index_of_diagonal_entries[i];
1625 
1626  if (unsigned(*(column_index_pt+diag_index))!=i)
1627  {
1628  std::string err_strng="Gauss-Seidel cannot operate on a matrix with ";
1629  err_strng+="zero diagonal entries.";
1630  throw OomphLibError(err_strng,
1631  OOMPH_CURRENT_FUNCTION,
1632  OOMPH_EXCEPTION_LOCATION);
1633  }
1634  } // for (unsigned i=0;i<n_dof;i++)
1635 
1636  // If there were no zero diagonal entries we can calculate the vector
1637  // constant_term=inv(L+D)*rhs.
1638 
1639  // Create a vector to hold the constant term in the iteration
1640  DoubleVector constant_term(this->distribution_pt(),0.0);
1641 
1642  // Copy the entries of rhs into the vector, constant_term
1643  constant_term=rhs;
1644 
1645  // Start by looping over the rows of constant_term
1646  for (unsigned i=0;i<n_dof;i++)
1647  {
1648  // Get the index of the last entry below or on the diagonal in row i
1649  unsigned diag_index=Index_of_diagonal_entries[i];
1650 
1651  // Get the index of the first entry in row i
1652  unsigned row_i_start=*(row_start_pt+i);
1653 
1654  // If there are entries strictly below the diagonal then the
1655  // column index of the first entry in the i-th row cannot be
1656  // i (since we have ensured that nonzero entries exist on the
1657  // diagonal so this is the largest the column index of the
1658  // first entry in this row could be)
1659  if (unsigned(*(column_index_pt+row_i_start))!=i)
1660  {
1661  // Initialise the column index variable
1662  unsigned column_index=0;
1663 
1664  // Loop over the entries below the diagonal on row i
1665  for (unsigned j=row_i_start;j<diag_index;j++)
1666  {
1667  // Find the column index of this nonzero entry
1668  column_index=*(column_index_pt+j);
1669 
1670  // Subtract the value of A(i,j)*constant_term(j) (where j<i)
1671  constant_term[i]-=(*(value_pt+j))*constant_term[column_index];
1672  }
1673  } // if (*(column_index_pt+row_i_start)!=i)
1674 
1675  // Finish off by calculating constant_term(i)/A(i,i)
1676  constant_term[i]/=*(value_pt+diag_index);
1677  } // for (unsigned i=0;i<n_dof;i++)
1678 
1679  // Build a temporary vector to store the value of A*x with each iteration
1680  DoubleVector temp_vec(this->distribution_pt(),0.0);
1681 
1682  // Outermost loop: Run max_iter times (the iteration number)
1683  for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
1684  {
1685  // With each iteration we need to calculate
1686  // x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
1687  // =-inv(L+D)*U*x_old+constant_term.
1688  // Since we've already calculated the vector constant_term, it remains
1689  // for us to calculate inv(L+D)*U*x. To do this we break the calculation
1690  // into two parts:
1691  // (1) The matrix-vector multiplication temp_vec=U*x, and;
1692  // (2) The matrix-vector multiplication inv(L+D)*temp_vec.
1693 
1694  //------------------------------
1695  // (1) Calculating temp_vec=U*x:
1696  //------------------------------
1697 
1698  // Auxiliary variable to store the index of the first entry on the
1699  // upper triangular portion of the matrix
1700  unsigned upper_tri_start=0;
1701 
1702  // Auxiliary variable to store the index of the first entry in the
1703  // next row
1704  unsigned next_row_start=0;
1705 
1706  // Set the temporary vector temp_vec to be initially be the zero vector
1707  temp_vec.initialise(0.0);
1708 
1709  // Loop over the rows of temp_vec (note: we do not loop over the final
1710  // row since all the entries on the last row of a strict upper triangular
1711  // matrix are zero)
1712  for (unsigned i=0;i<n_dof-1;i++)
1713  {
1714  // Get the index of the first entry on the upper triangular portion of
1715  // the matrix noting that the i-th entry in Index_of_diagonal_entries
1716  // will store the index of the diagonal entry of the i-th row
1717  upper_tri_start=Index_of_diagonal_entries[i]+1;
1718 
1719  // Get the index of the first entry in the next row
1720  next_row_start=*(row_start_pt+i+1);
1721 
1722  // Variable to store the column index of each entry
1723  unsigned column_index=0;
1724 
1725  // Loop over all of the entries above the diagonal
1726  for (unsigned j=upper_tri_start;j<next_row_start;j++)
1727  {
1728  // Get the column index of this entry
1729  column_index=*(column_index_pt+j);
1730 
1731  // Update the value of temp_vec[i] by adding A(i,j)*x(j) (where j>i)
1732  temp_vec[i]+=(*(value_pt+j))*x[column_index];
1733  }
1734  } // for (unsigned i=0;i<n_dof;i++)
1735 
1736  //-----------------------------------
1737  // (2) Calculating inv(L+D)*temp_vec:
1738  //-----------------------------------
1739 
1740  // Now U*x=temp_vec has been calculated we need to calculate the vector
1741  // y=inv(L+D)*U*x=inv(L+D)*temp_vec.
1742  // Note: the i-th entry in y is given by
1743  // y(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*y(j))/A(i,i).
1744  // We can see from this that to calculate the current entry in y we
1745  // use its previously calculated entries and the current entry in
1746  // temp_vec. As a result, we do not need two separate vectors y and
1747  // temp_vec for this calculation. Instead, we can just update the
1748  // entries in temp_vec (carefully!) using:
1749  // temp_vec(i)=(temp_vec(i)-\sum_{j=0}^{i-1}A(i,j)*temp_vec(j))/A(i,i).
1750 
1751  // Start by looping over the rows of rhs
1752  for (unsigned i=0;i<n_dof;i++)
1753  {
1754  // Get the index of the last entry below or on the diagonal in row i
1755  unsigned diag_index=Index_of_diagonal_entries[i];
1756 
1757  // Get the index of the first entry in row i
1758  unsigned row_i_start=*(row_start_pt+i);
1759 
1760  // If there are no entries strictly below the diagonal then the
1761  // column index of the first entry in the i-th row will be greater
1762  // than or equal to i
1763  if (unsigned(*(column_index_pt+row_i_start))!=i)
1764  {
1765  // Initialise the column index variable
1766  unsigned column_index=0;
1767 
1768  // Loop over the entries below the diagonal
1769  for (unsigned j=row_i_start;j<diag_index;j++)
1770  {
1771  // Find the column index of this nonzero entry
1772  column_index=*(column_index_pt+j);
1773 
1774  // Subtract the value of A(i,j)*y(j) (where j<i)
1775  temp_vec[i]-=(*(value_pt+j))*temp_vec[column_index];
1776  }
1777  } // if (*(column_index_pt+row_i_start)!=i)
1778 
1779  // Finish off by calculating temp_vec(i)/A(i,i)
1780  temp_vec[i]/=*(value_pt+diag_index);
1781  } // for (unsigned i=0;i<n_dof;i++)
1782 
1783  //-------------------------
1784  // Updating the solution x:
1785  //-------------------------
1786 
1787  // The updated solution is given by
1788  // x_new=-inv(L+D)*U*x_old+inv(L+D)*rhs,
1789  // =-temp_vec+constant_term,
1790  // We have calculated both constant_term and temp_vec so we simply
1791  // need to calculate x now
1792 
1793  // Set x to be the vector, constant_term
1794  x=constant_term;
1795 
1796  // Update x to give the next iterate
1797  x-=temp_vec;
1798 
1799  // Increment the value of Iterations
1800  Iterations++;
1801 
1802  // Calculate the residual only if we're not inside the multigrid solver
1803  if (!Use_as_smoother)
1804  {
1805  // Get residual
1806  matrix_pt->residual(x,rhs,current_residual);
1807 
1808  // Calculate the relative residual norm (i.e. \frac{\|r_{i}\|}{\|r_{0}\|})
1809  norm_res=current_residual.norm()/norm_f;
1810 
1811  // If required will document convergence history to screen or file (if
1812  // stream open)
1814  {
1815  if (!Output_file_stream.is_open())
1816  {
1817  oomph_info << iter_num << " "
1818  << norm_res << std::endl;
1819  }
1820  else
1821  {
1822  Output_file_stream << iter_num << " "
1823  << norm_res << std::endl;
1824  }
1825  } // if (Doc_convergence_history)
1826 
1827  // Check the tolerance only if the residual norm is being computed
1828  if (norm_res<Tolerance)
1829  {
1830  // Break out of the for-loop
1831  break;
1832  }
1833  } // if (!Use_as_smoother)
1834  } // for (unsigned iter_num=0;iter_num<max_iter;iter_num++)
1835 
1836  // Calculate the residual only if we're not inside the multigrid solver
1837  if (!Use_as_smoother)
1838  {
1839  // If time documentation is enabled
1840  if(Doc_time)
1841  {
1842  oomph_info << "\nGS converged. Residual norm: " << norm_res
1843  << "\nNumber of iterations to convergence: " << Iterations
1844  << "\n" << std::endl;
1845  }
1846  } // if (!Use_as_smoother)
1847 
1848  // Copy the solution into the solution vector
1849  solution=x;
1850 
1851  // Document the time for the solver
1852  double t_end=TimingHelpers::timer();
1853  Solution_time=t_end-t_start;
1854 
1855  // If time documentation is enabled
1856  if(Doc_time)
1857  {
1858  oomph_info << "Time for solve with Gauss-Seidel [sec]: "
1859  << Solution_time << std::endl;
1860  }
1861 
1862  // If the solver failed to converge and the user asked for an error if
1863  // this happened
1864  if((Iterations>Max_iter-1)&&(Throw_error_after_max_iter))
1865  {
1866  std::string error_message="Solver failed to converge and you requested ";
1867  error_message+="an error on convergence failures.";
1868  throw OomphLibError(error_message,
1869  OOMPH_EXCEPTION_LOCATION,
1870  OOMPH_CURRENT_FUNCTION);
1871  }
1872  } // End of solve_helper
1873 
1874 
1875 //////////////////////////////////////////////////////////////////////
1876 //////////////////////////////////////////////////////////////////////
1877 //////////////////////////////////////////////////////////////////////
1878 
1879 
1880 //==================================================================
1881 /// \short Solver: Takes pointer to problem and returns the results
1882 /// vector which contains the solution of the linear system defined
1883 /// by the problem's fully assembled Jacobian and residual vector.
1884 //==================================================================
1885  template<typename MATRIX>
1886  void DampedJacobi<MATRIX>::solve(Problem* const &problem_pt,
1887  DoubleVector &result)
1888  {
1889  // Reset the Use_as_smoother_flag as the solver is not being used
1890  // as a smoother
1891  Use_as_smoother=false;
1892 
1893  // Find the # of degrees of freedom (variables)
1894  unsigned n_dof=problem_pt->ndof();
1895 
1896  // Initialise timer
1897  double t_start=TimingHelpers::timer();
1898 
1899  // We're not re-solving
1900  Resolving=false;
1901 
1902  // Get rid of any previously stored data
1903  clean_up_memory();
1904 
1905  // Set up the distribution
1906  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),
1907  n_dof,false);
1908 
1909  // Assign the distribution to the LinearSolver
1910  this->build_distribution(dist);
1911 
1912  // Allocate space for the Jacobian matrix in format specified
1913  // by template parameter
1914  Matrix_pt=new MATRIX;
1915 
1916  // Get the nonlinear residuals vector
1917  DoubleVector f;
1918 
1919  // Assign the Jacobian and the residuals vector
1920  problem_pt->get_jacobian(f,*Matrix_pt);
1921 
1922  // Extract the diagonal entries of the matrix and store them
1923  extract_diagonal_entries(Matrix_pt);
1924 
1925  // We've made the matrix, we can delete it...
1926  Matrix_can_be_deleted=true;
1927 
1928  // Doc time for setup
1929  double t_end=TimingHelpers::timer();
1930  Jacobian_setup_time=t_end-t_start;
1931 
1932  // If time documentation is enabled
1933  if(Doc_time)
1934  {
1935  oomph_info << "Time for setup of Jacobian [sec]: "
1936  << Jacobian_setup_time << std::endl;
1937  }
1938 
1939  // Call linear algebra-style solver
1940  solve_helper(Matrix_pt,f,result);
1941 
1942  // Kill matrix unless it's still required for resolve
1943  if (!Enable_resolve) clean_up_memory();
1944  } // End of solve
1945 
1946 //==================================================================
1947 /// \short Linear-algebra-type solver: Takes pointer to a matrix and
1948 /// rhs vector and returns the solution of the linear system.
1949 //==================================================================
1950  template<typename MATRIX>
1952  const DoubleVector& rhs,
1953  DoubleVector& solution)
1954  {
1955  // Get number of dofs
1956  unsigned n_dof=rhs.nrow();
1957 
1958 #ifdef PARANOID
1959  // Upcast the matrix to the appropriate type
1960  MATRIX* tmp_matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
1961 
1962  // PARANOID Run the self-tests to check the inputs are correct
1963  this->check_validity_of_solve_helper_inputs<MATRIX>(tmp_matrix_pt,rhs,
1964  solution,n_dof);
1965 
1966  // We don't need the pointer any more but we do still need the matrix
1967  // so just make tmp_matrix_pt a null pointer
1968  tmp_matrix_pt=0;
1969 #endif
1970 
1971  // Setup the solution if it is not
1972  if (!solution.distribution_pt()->built())
1973  {
1974  // Build the distribution of the solution vector if it hasn't been done yet
1975  solution.build(this->distribution_pt(),0.0);
1976  }
1977  // If the solution has already been set up
1978  else
1979  {
1980  // If we're inside the multigrid solver then as we traverse up the
1981  // hierarchy we use the smoother on the updated approximate solution.
1982  // As such, we should ONLY be resetting all the values to zero if
1983  // we're NOT inside the multigrid solver
1984  if (!Use_as_smoother)
1985  {
1986  // Initialise the vector with all entries set to zero
1987  solution.initialise(0.0);
1988  }
1989  } // if (!solution.distribution_pt()->built())
1990 
1991  // Initialise timer
1992  double t_start = TimingHelpers::timer();
1993 
1994  // Initial guess isn't necessarily zero (restricted solution from finer
1995  // grids) therefore x needs to be assigned values from the input.
1996  DoubleVector x(solution);
1997 
1998  // Create a vector to store the value of the constant term, omega*inv(D)*r
1999  DoubleVector constant_term(this->distribution_pt(),0.0);
2000 
2001  // Calculate the constant term vector
2002  for(unsigned i=0;i<n_dof;i++)
2003  {
2004  // Assign the i-th entry of constant_term
2005  constant_term[i]=Omega*Matrix_diagonal[i]*rhs[i];
2006  }
2007 
2008  // Create a vector to hold the residual. This will only be built if
2009  // we're not inside the multigrid solver
2010  DoubleVector local_residual;
2011 
2012  // Variable to store the 2-norm of the residual vector. Only used
2013  // if we are not working inside the MG solver
2014  double norm_res=0.0;
2015 
2016  // Variables to hold the initial residual norm. Only used if we're
2017  // not inside the multigrid solver
2018  double norm_f=0.0;
2019 
2020  // Initialise the value of Iterations
2021  Iterations=0;
2022 
2023  // Calculate the residual only if we're not inside the multigrid solver
2024  if (!Use_as_smoother)
2025  {
2026  // Build the local residual vector
2027  local_residual.build(this->distribution_pt(),0.0);
2028 
2029  // Calculate the residual vector
2030  matrix_pt->residual(x,rhs,local_residual);
2031 
2032  // Calculate the 2-norm
2033  norm_res=local_residual.norm();
2034 
2035  // Store the initial norm
2036  norm_f=norm_res;
2037 
2038  // If required will document convergence history to screen
2039  // or file (if stream is open)
2040  if (Doc_convergence_history)
2041  {
2042  if (!Output_file_stream.is_open())
2043  {
2044  oomph_info << Iterations << " " << norm_res << std::endl;
2045  }
2046  else
2047  {
2048  Output_file_stream << Iterations << " " << norm_res <<std::endl;
2049  }
2050  } // if (Doc_convergence_history)
2051  } // if (!Use_as_smoother)
2052 
2053  // Create a temporary vector to store the value of A*e
2054  // on each iteration and another to store the residual
2055  // at the end of each iteration
2056  DoubleVector temp_vec(this->distribution_pt(),0.0);
2057 
2058  // Outermost loop: Run up to Max_iter times (the iteration number)
2059  for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
2060  {
2061  // Calculate A*e
2062  matrix_pt->multiply(x,temp_vec);
2063 
2064  // Loop over each degree of freedom and update
2065  // the current approximation
2066  for (unsigned idof=0;idof<n_dof;idof++)
2067  {
2068  // Scale the idof'th entry of temp_vec
2069  // by omega/A(i,i)
2070  temp_vec[idof]*=Omega*Matrix_diagonal[idof];
2071  }
2072 
2073  // Update the value of e (in the system Ae=r)
2074  x+=constant_term;
2075  x-=temp_vec;
2076 
2077  // Increment the value of Iterations
2078  Iterations++;
2079 
2080  // Calculate the residual only if we're not inside the multigrid solver
2081  if (!Use_as_smoother)
2082  {
2083  // Get residual
2084  matrix_pt->residual(x,rhs,local_residual);
2085 
2086  // Calculate the 2-norm of the residual r=b-Ax
2087  norm_res=local_residual.norm()/norm_f;
2088 
2089  // If required, this will document convergence history to
2090  // screen or file (if the stream is open)
2091  if (Doc_convergence_history)
2092  {
2093  if (!Output_file_stream.is_open())
2094  {
2095  oomph_info << Iterations << " "
2096  << norm_res << std::endl;
2097  }
2098  else
2099  {
2100  Output_file_stream << Iterations << " "
2101  << norm_res << std::endl;
2102  }
2103  } // if (Doc_convergence_history)
2104 
2105  // Check the tolerance only if the residual norm is being computed
2106  if (norm_res<Tolerance)
2107  {
2108  // Break out of the for-loop
2109  break;
2110  }
2111  } // if (!Use_as_smoother)
2112  } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
2113 
2114  // Calculate the residual only if we're not inside the multigrid solver
2115  if (!Use_as_smoother)
2116  {
2117  // If time documentation is enabled
2118  if(Doc_time)
2119  {
2120  oomph_info << "\nDamped Jacobi converged. Residual norm: " << norm_res
2121  << "\nNumber of iterations to convergence: " << Iterations
2122  << "\n" << std::endl;
2123  }
2124  } // if (!Use_as_smoother)
2125 
2126  // Copy the solution into the solution vector
2127  solution=x;
2128 
2129  // Doc. time for solver
2130  double t_end=TimingHelpers::timer();
2131  Solution_time=t_end-t_start;
2132  if(Doc_time)
2133  {
2134  oomph_info << "Time for solve with damped Jacobi [sec]: "
2135  << Solution_time << std::endl;
2136  }
2137 
2138  // If the solver failed to converge and the user asked for an error if
2139  // this happened
2140  if((Iterations>Max_iter-1)&&(Throw_error_after_max_iter))
2141  {
2142  std::string error_message="Solver failed to converge and you requested ";
2143  error_message+="an error on convergence failures.";
2144  throw OomphLibError(error_message,
2145  OOMPH_EXCEPTION_LOCATION,
2146  OOMPH_CURRENT_FUNCTION);
2147  }
2148  } // End of solve_helper
2149 
2150 
2151 ///////////////////////////////////////////////////////////////////////
2152 ///////////////////////////////////////////////////////////////////////
2153 ///////////////////////////////////////////////////////////////////////
2154 
2155 
2156 //==================================================================
2157 /// \Short Re-solve the system defined by the last assembled Jacobian
2158 /// and the rhs vector specified here. Solution is returned in
2159 /// the vector result.
2160 //==================================================================
2161  template<typename MATRIX>
2163  DoubleVector &result)
2164  {
2165  // We are re-solving
2166  Resolving=true;
2167 
2168 #ifdef PARANOID
2169  if (Matrix_pt==0)
2170  {
2171  throw OomphLibError("No matrix was stored -- cannot re-solve",
2172  OOMPH_CURRENT_FUNCTION,
2173  OOMPH_EXCEPTION_LOCATION);
2174  }
2175 #endif
2176 
2177  // Call linear algebra-style solver
2178  this->solve(Matrix_pt,rhs,result);
2179 
2180  // Reset re-solving flag
2181  Resolving=false;
2182  }
2183 
2184 
2185 //==================================================================
2186 /// Solver: Takes pointer to problem and returns the results vector
2187 /// which contains the solution of the linear system defined by
2188 /// the problem's fully assembled Jacobian and residual vector.
2189 //==================================================================
2190  template<typename MATRIX>
2191  void GMRES<MATRIX>::solve(Problem* const &problem_pt, DoubleVector &result)
2192  {
2193 
2194  //Find # of degrees of freedom (variables)
2195  unsigned n_dof = problem_pt->ndof();
2196 
2197  // Initialise timer
2198  double t_start = TimingHelpers::timer();
2199 
2200  // We're not re-solving
2201  Resolving=false;
2202 
2203  // Get rid of any previously stored data
2204  clean_up_memory();
2205 
2206  // setup the distribution
2207  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),
2208  n_dof,false);
2209  this->build_distribution(dist);
2210 
2211  // Get Jacobian matrix in format specified by template parameter
2212  // and nonlinear residual vector
2213  Matrix_pt=new MATRIX;
2214  DoubleVector f;
2215  if (dynamic_cast<DistributableLinearAlgebraObject*>(Matrix_pt) != 0)
2216  {
2217  if (dynamic_cast<CRDoubleMatrix*>(Matrix_pt) != 0)
2218  {
2219  dynamic_cast<CRDoubleMatrix* >(Matrix_pt)->build(this->distribution_pt());
2220  f.build(this->distribution_pt(),0.0);
2221  }
2222  }
2223  problem_pt->get_jacobian(f,*Matrix_pt);
2224 
2225  // We've made the matrix, we can delete it...
2226  Matrix_can_be_deleted=true;
2227 
2228  // Doc time for setup
2229  double t_end = TimingHelpers::timer();
2230  Jacobian_setup_time= t_end-t_start;
2231 
2232  if(Doc_time)
2233  {
2234  oomph_info << "Time for setup of Jacobian [sec]: "
2235  << Jacobian_setup_time << std::endl;
2236  }
2237 
2238  // Call linear algebra-style solver
2239  //If the result distribution is wrong, then redistribute
2240  //before the solve and return to original distribution
2241  //afterwards
2242  if((result.built()) &&
2243  (!(*result.distribution_pt() == *this->distribution_pt())))
2244  {
2246  temp_global_dist(result.distribution_pt());
2247  result.build(this->distribution_pt(),0.0);
2248  this->solve_helper(Matrix_pt,f,result);
2249  result.redistribute(&temp_global_dist);
2250  }
2251  //Otherwise just solve
2252  else
2253  {
2254  this->solve_helper(Matrix_pt,f,result);
2255  }
2256 
2257  // Kill matrix unless it's still required for resolve
2258  if (!Enable_resolve) clean_up_memory();
2259 
2260  };
2261 
2262 
2263 //=============================================================================
2264 /// Linear-algebra-type solver: Takes pointer to a matrix and rhs vector
2265 /// and returns the solution of the linear system.
2266 /// based on the algorithm presented in Templates for the
2267 /// Solution of Linear Systems: Building Blocks for Iterative Methods, Barrett,
2268 /// Berry et al, SIAM, 2006 and the implementation in the IML++ library :
2269 /// http://math.nist.gov/iml++/
2270 //=============================================================================
2271  template <typename MATRIX>
2273  const DoubleVector &rhs,
2274  DoubleVector &solution)
2275  {
2276 
2277  // Get number of dofs
2278  unsigned n_dof=rhs.nrow();
2279 
2280 #ifdef PARANOID
2281  // PARANOID check that if the matrix is distributable then it should not be
2282  // then it should not be distributed
2283  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
2284  {
2285  if (dynamic_cast<DistributableLinearAlgebraObject*>
2286  (matrix_pt)->distributed())
2287  {
2288  std::ostringstream error_message_stream;
2289  error_message_stream
2290  << "The matrix must not be distributed.";
2291  throw OomphLibError(error_message_stream.str(),
2292  OOMPH_CURRENT_FUNCTION,
2293  OOMPH_EXCEPTION_LOCATION);
2294  }
2295  }
2296  // PARANOID check that this rhs distribution is setup
2297  if (!rhs.built())
2298  {
2299  std::ostringstream error_message_stream;
2300  error_message_stream
2301  << "The rhs vector distribution must be setup.";
2302  throw OomphLibError(error_message_stream.str(),
2303  OOMPH_CURRENT_FUNCTION,
2304  OOMPH_EXCEPTION_LOCATION);
2305  }
2306  // PARANOID check that the rhs has the right number of global rows
2307  if(rhs.nrow() != n_dof)
2308  {
2309  throw OomphLibError(
2310  "RHS does not have the same dimension as the linear system",
2311  OOMPH_CURRENT_FUNCTION,
2312  OOMPH_EXCEPTION_LOCATION);
2313  }
2314  // PARANOID check that the rhs is not distributed
2315  if (rhs.distribution_pt()->distributed())
2316  {
2317  std::ostringstream error_message_stream;
2318  error_message_stream
2319  << "The rhs vector must not be distributed.";
2320  throw OomphLibError(error_message_stream.str(),
2321  OOMPH_CURRENT_FUNCTION,
2322  OOMPH_EXCEPTION_LOCATION);
2323  }
2324  // PARANOID check that if the result is setup it matches the distribution
2325  // of the rhs
2326  if (solution.built())
2327  {
2328  if (!(*rhs.distribution_pt() == *solution.distribution_pt()))
2329  {
2330  std::ostringstream error_message_stream;
2331  error_message_stream
2332  << "If the result distribution is setup then it must be the same as the "
2333  << "rhs distribution";
2334  throw OomphLibError(error_message_stream.str(),
2335  OOMPH_CURRENT_FUNCTION,
2336  OOMPH_EXCEPTION_LOCATION);
2337  }
2338  }
2339 #endif
2340 
2341  // Set up the solution if it is not
2342  if (!solution.built())
2343  {
2344  solution.build(this->distribution_pt(),0.0);
2345  }
2346  // Otherwise initialise to zero
2347  else
2348  {
2349  solution.initialise(0.0);
2350  }
2351 
2352  // Time solver
2353  double t_start=TimingHelpers::timer();
2354 
2355  // Relative residual
2356  double resid;
2357 
2358  // iteration counter
2359  unsigned iter = 1;
2360 
2361  // if not using iteration restart set Restart to n_dof
2362  if (!Iteration_restart)
2363  {
2364  Restart = n_dof;
2365  }
2366 
2367  // initialise vectors
2368  Vector<double> s(Restart + 1,0);
2369  Vector<double> cs(Restart + 1);
2370  Vector<double> sn(Restart + 1);
2371  DoubleVector w(this->distribution_pt(),0.0);
2372 
2373  // Setup preconditioner only if we're not re-solving
2374  if (!Resolving)
2375  {
2376  // only setup the preconditioner before solve if require
2377  if (Setup_preconditioner_before_solve)
2378  {
2379 
2380  //Setup preconditioner from the Jacobian matrix
2381  double t_start_prec = TimingHelpers::timer();
2382 
2383  // do not setup
2384  preconditioner_pt()->setup(matrix_pt);
2385 
2386  // Doc time for setup of preconditioner
2387  double t_end_prec = TimingHelpers::timer();
2388  Preconditioner_setup_time = t_end_prec-t_start_prec;
2389 
2390  if(Doc_time)
2391  {
2392  oomph_info << "Time for setup of preconditioner [sec]: "
2393  << Preconditioner_setup_time << std::endl;
2394  }
2395  }
2396  }
2397  else
2398  {
2399  if(Doc_time)
2400  {
2401  oomph_info << "Setup of preconditioner is bypassed in resolve mode"
2402  << std::endl;
2403  }
2404  }
2405 
2406  // solve b-Jx = Mr for r (assumes x = 0);
2407  DoubleVector r(this->distribution_pt(),0.0);
2408  if(Preconditioner_LHS)
2409  {
2410  preconditioner_pt()->preconditioner_solve(rhs,r);
2411  }
2412  else
2413  {
2414  r=rhs;
2415  }
2416  double normb = 0;
2417 
2418  // compute norm(r)
2419  double* r_pt = r.values_pt();
2420  for (unsigned i = 0; i < n_dof; i++)
2421  {
2422  normb += r_pt[i] * r_pt[i];
2423  }
2424  normb = sqrt(normb);
2425 
2426  // set beta (the initial residual)
2427  double beta = normb;
2428 
2429  // compute initial relative residual
2430  if (normb == 0.0) normb = 1;
2431  resid = beta / normb;
2432 
2433  // if required will document convergence history to screen or file (if
2434  // stream open)
2435  if (Doc_convergence_history)
2436  {
2437  if (!Output_file_stream.is_open())
2438  {
2439  oomph_info << 0 << " " << resid << std::endl;
2440  }
2441  else
2442  {
2443  Output_file_stream << 0 << " " << resid <<std::endl;
2444  }
2445  }
2446 
2447  // if GMRES converges immediately
2448  if (resid <= Tolerance)
2449  {
2450  if(Doc_time)
2451  {
2452  oomph_info << "GMRES converged immediately. Normalised residual norm: "
2453  << resid << std::endl;
2454  }
2455  // Doc time for solver
2456  double t_end = TimingHelpers::timer();
2457  Solution_time = t_end-t_start;
2458 
2459  if(Doc_time)
2460  {
2461  oomph_info << "Time for solve with GMRES [sec]: "
2462  << Solution_time << std::endl;
2463  }
2464  return;
2465  }
2466 
2467 
2468  // initialise vector of orthogonal basis vectors (v) and upper hessenberg
2469  // matrix H
2470  // NOTE: for implementation purpose the upper hessenberg matrix indexes are
2471  // are swapped so the matrix is effectively transposed
2473  v.resize(Restart + 1);
2474  Vector<Vector<double> > H(Restart + 1);
2475 
2476  // while...
2477  while (iter <= Max_iter)
2478  {
2479 
2480  // set zeroth basis vector v[0] to r/beta
2481  v[0].build(this->distribution_pt(),0.0);
2482  double* v0_pt = v[0].values_pt();
2483  for (unsigned i = 0; i < n_dof; i++)
2484  {
2485  v0_pt[i] = r_pt[i] / beta;
2486  }
2487 
2488  //
2489  s[0] = beta;
2490 
2491  // inner iteration counter for restarted version
2492  unsigned iter_restart;
2493 
2494  // perform iterations
2495  for (iter_restart = 0; iter_restart < Restart && iter <= Max_iter;
2496  iter_restart++, iter++)
2497  {
2498 
2499  // resize next column of upper hessenberg matrix
2500  H[iter_restart].resize(iter_restart+2);
2501 
2502  // solve Jv[i] = Mw for w
2503  {
2504  DoubleVector temp(this->distribution_pt(),0.0);
2505  if(Preconditioner_LHS)
2506  {
2507  // solve Jv[i] = Mw for w
2508  matrix_pt->multiply(v[iter_restart],temp);
2509  preconditioner_pt()->preconditioner_solve(temp,w);
2510  }
2511  else
2512  {
2513  //w=JM^{-1}v by saad p270
2514  preconditioner_pt()->preconditioner_solve(v[iter_restart],temp);
2515  matrix_pt->multiply(temp,w);
2516  }
2517  }
2518 
2519  //
2520  double* w_pt = w.values_pt();
2521  for (unsigned k = 0; k <= iter_restart; k++)
2522  {
2523 
2524  //
2525  H[iter_restart][k] = 0.0;
2526  double* vk_pt = v[k].values_pt();
2527  for (unsigned i = 0; i < n_dof; i++)
2528  {
2529  H[iter_restart][k] += w_pt[i]*vk_pt[i];
2530  }
2531 
2532  //
2533  for (unsigned i = 0; i < n_dof; i++)
2534  {
2535  w_pt[i] -= H[iter_restart][k] * vk_pt[i];
2536  }
2537  }
2538 
2539  //
2540  {
2541  double temp_norm_w = 0.0;
2542  for (unsigned i = 0; i < n_dof; i++)
2543  {
2544  temp_norm_w += w_pt[i]*w_pt[i];
2545  }
2546  temp_norm_w = sqrt(temp_norm_w);
2547  H[iter_restart][iter_restart+1] = temp_norm_w;
2548  }
2549 
2550  //
2551  v[iter_restart + 1].build(this->distribution_pt(),0.0);
2552  double* v_pt = v[iter_restart + 1].values_pt();
2553  for (unsigned i = 0; i < n_dof; i++)
2554  {
2555  v_pt[i] = w_pt[i] / H[iter_restart][iter_restart+1];
2556  }
2557 
2558  //
2559  for (unsigned k = 0; k < iter_restart; k++)
2560  {
2561  apply_plane_rotation(H[iter_restart][k], H[iter_restart][k+1], cs[k],
2562  sn[k]);
2563  }
2564  generate_plane_rotation(H[iter_restart][iter_restart],
2565  H[iter_restart][iter_restart+1],
2566  cs[iter_restart],
2567  sn[iter_restart]);
2568  apply_plane_rotation(H[iter_restart][iter_restart],
2569  H[iter_restart][iter_restart+1], cs[iter_restart],
2570  sn[iter_restart]);
2571  apply_plane_rotation(s[iter_restart],s[iter_restart+1],cs[iter_restart],
2572  sn[iter_restart]);
2573 
2574  // compute current residual
2575  beta = std::fabs(s[iter_restart+1]);
2576 
2577  // compute relative residual
2578  resid = beta/normb;
2579 
2580  // if required will document convergence history to screen or file (if
2581  // stream open)
2582  if (Doc_convergence_history)
2583  {
2584  if (!Output_file_stream.is_open())
2585  {
2586  oomph_info << iter << " "<< resid << std::endl;
2587  }
2588  else
2589  {
2590  Output_file_stream << iter << " " << resid <<std::endl;
2591  }
2592  }
2593 
2594  // if required tolerance found
2595  if (resid < Tolerance)
2596  {
2597  // update result vector
2598  update(iter_restart, H, s, v, solution);
2599 
2600  // document convergence
2601  if(Doc_time)
2602  {
2603  oomph_info << std::endl;
2604  oomph_info << "GMRES converged (1). Normalised residual norm: "
2605  << resid << std::endl;
2606  oomph_info << "Number of iterations to convergence: "
2607  << iter << std::endl;
2608  oomph_info << std::endl;
2609  }
2610 
2611  // Doc time for solver
2612  double t_end = TimingHelpers::timer();
2613  Solution_time = t_end-t_start;
2614 
2615  Iterations = iter;
2616 
2617  if(Doc_time)
2618  {
2619  oomph_info << "Time for solve with GMRES [sec]: "
2620  << Solution_time << std::endl;
2621  }
2622  return;
2623  }
2624  }
2625 
2626  // update
2627  if (iter_restart>0) update((iter_restart - 1), H, s, v, solution);
2628 
2629  // solve Mr = (b-Jx) for r
2630  {
2631  DoubleVector temp(this->distribution_pt(),0.0);
2632  matrix_pt->multiply(solution,temp);
2633  double* temp_pt = temp.values_pt();
2634  const double* rhs_pt = rhs.values_pt();
2635  for (unsigned i = 0; i < n_dof; i++)
2636  {
2637  temp_pt[i] = rhs_pt[i] - temp_pt[i];
2638  }
2639 
2640  if(Preconditioner_LHS)
2641  {
2642  preconditioner_pt()->preconditioner_solve(temp,r);
2643  }
2644  }
2645 
2646  // compute current residual
2647  beta = 0.0;
2648  r_pt = r.values_pt();
2649  for (unsigned i = 0; i < n_dof; i++)
2650  {
2651  beta += r_pt[i] * r_pt[i];
2652  }
2653  beta = sqrt(beta);
2654 
2655  // if relative residual within tolerance
2656  resid = beta / normb;
2657  if (resid < Tolerance)
2658  {
2659  if(Doc_time)
2660  {
2661  oomph_info << std::endl;
2662  oomph_info << "GMRES converged (2). Normalised residual norm: "
2663  << resid << std::endl;
2664  oomph_info << "Number of iterations to convergence: "
2665  << iter << std::endl;
2666  oomph_info << std::endl;
2667  }
2668 
2669  // Doc time for solver
2670  double t_end = TimingHelpers::timer();
2671  Solution_time = t_end-t_start;
2672 
2673  Iterations = iter;
2674 
2675  if(Doc_time)
2676  {
2677  oomph_info << "Time for solve with GMRES [sec]: "
2678  << Solution_time << std::endl;
2679  }
2680  return;
2681  }
2682  }
2683 
2684 
2685  // otherwise GMRES failed convergence
2686  oomph_info << std::endl;
2687  oomph_info << "GMRES did not converge to required tolerance! "
2688  << std::endl;
2689  oomph_info << "Returning with normalised residual norm: "
2690  << resid << std::endl;
2691  oomph_info << "after " << Max_iter<< " iterations." << std::endl;
2692  oomph_info << std::endl;
2693 
2694 
2695  if(Throw_error_after_max_iter)
2696  {
2697  std::string err = "Solver failed to converge and you requested an error";
2698  err += " on convergence failures.";
2699  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
2700  OOMPH_CURRENT_FUNCTION);
2701  }
2702 
2703  return;
2704 
2705  } // End GMRES
2706 
2707 
2708 
2709 
2710 //Ensure build of required objects
2711  template class BiCGStab<CCDoubleMatrix>;
2712  template class BiCGStab<CRDoubleMatrix>;
2713  template class BiCGStab<DenseDoubleMatrix>;
2714 
2715  template class CG<CCDoubleMatrix>;
2716  template class CG<CRDoubleMatrix>;
2717  template class CG<DenseDoubleMatrix>;
2718 
2719  template class GS<CCDoubleMatrix>;
2720  template class GS<CRDoubleMatrix>;
2721  template class GS<DenseDoubleMatrix>;
2722 
2723  template class DampedJacobi<CCDoubleMatrix>;
2724  template class DampedJacobi<CRDoubleMatrix>;
2725  template class DampedJacobi<DenseDoubleMatrix>;
2726 
2727  template class GMRES<CCDoubleMatrix>;
2728  template class GMRES<CRDoubleMatrix>;
2729  template class GMRES<DenseDoubleMatrix>;
2730 
2731 // Solvers for SumOfMatrices class
2732  template class BiCGStab<SumOfMatrices>;
2733  template class CG<SumOfMatrices>;
2734  template class GS<SumOfMatrices>;
2735  template class GMRES<SumOfMatrices>;
2736 
2737 
2738 }
bool Throw_error_after_max_iter
Should we throw an error instead of just returning when we hit the max iterations?
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1056
bool built() const
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here...
static IdentityPreconditioner Default_preconditioner
Default preconditioner: The base class for preconditioners is a fully functional (if trivial!) precon...
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 Tolerance
Convergence tolerance.
double * values_pt()
access function to the underlying values
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
cstr elem_len * i
Definition: cfortran.h:607
unsigned Max_iter
Maximum number of iterations.
unsigned nrow() const
access function to the number of global rows.
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:80
The Problem class.
Definition: problem.h:152
virtual void multiply(const DoubleVector &x, DoubleVector &soln) const =0
Multiply the matrix by the vector x: soln=Ax.
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...
char t
Definition: cfortran.h:572
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
bool distributed() const
distribution is serial or distributed
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here...
OomphInfo oomph_info
std::ofstream Output_file_stream
Output file stream for convergence history.
The GMRES method.
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed ...
double * value()
Access to C-style value array.
Definition: matrices.h:1062
The conjugate gradient method.
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.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
bool Use_as_smoother
When a derived class object is being used as a smoother in the MG solver (or elsewhere) the residual ...
double norm() const
compute the 2 norm of this vector
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)
Use damped Jacobi iteration as an IterativeLinearSolver: This obtains the Jacobian matrix J and the r...
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...
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
virtual void residual(const DoubleVector &x, const DoubleVector &b, DoubleVector &residual_)
Find the residual, i.e. r=b-Ax the residual.
Definition: matrices.h:343
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:84
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
double Jacobian_setup_time
Jacobian setup time.
void initialise(const double &v)
initialise the whole vector with value v
static char t char * s
Definition: cfortran.h:572
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1574
unsigned Max_iter
Max. # of Newton iterations.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here...
MATRIX * Matrix_pt
System matrix pointer in the format specified by the template argument.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
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...
double timer()
returns the time in seconds after some point in past
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
unsigned Iterations
Number of iterations taken.
double Solution_time
linear solver solution time
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
Definition: problem.cc:3851
The conjugate gradient method.
The Gauss Seidel method.
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
bool Doc_convergence_history
Flag indicating if the convergence history is to be documented.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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 build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
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...
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1050
Abstract base class for matrices of doubles – adds abstract interfaces for solving, LU decomposition and multiplication by vectors.
Definition: matrices.h:275
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
void build(const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)
build method: vector of values, vector of column indices, vector of row starts and number of rows and...
Definition: matrices.cc:1692
Explicit template specialisation of the Gauss Seidel method for compressed row format matrices...