mumps_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: 1110 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-01-19 07:46:59 +0000 (Tue, 19 Jan 2016) $
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 //Interface to MUMPS solver
31 
32 
33 #ifdef OOMPH_HAS_MPI
34 #include "mpi.h"
35 #endif
36 
37 
38 #include<iostream>
39 #include<vector>
40 
41 
42 //oomph-lib headers
43 #include "mumps_solver.h"
44 #include "Vector.h"
45 #include "oomph_utilities.h"
46 #include "problem.h"
47 
48 
49 namespace oomph
50 {
51 
52 
53 
54 //////////////////////////////////////////////////////////////////////////
55 //////////////////////////////////////////////////////////////////////////
56 //////////////////////////////////////////////////////////////////////////
57 
58 
59 
60  //=========================================================================
61  /// \short Default factor for workspace -- static so it can be overwritten
62  /// globally.
63  //=========================================================================
65 
66 
67 //=========================================================================
68 ///Static warning to suppress warnings about incorrect distribution of
69 ///RHS vector. Default is false
70 //=========================================================================
72 =false;
73 
74 //=============================================================================
75 /// Constructor: Call setup
76 //=============================================================================
78  {
79  Doc_stats=false;
80  Suppress_solve=false;
81  Delete_matrix_data=false;
87  }
88 
89 
90 //=============================================================================
91 /// Initialise instance of mumps data structure
92 //=============================================================================
94  {
95 
96  // Make new instance of Mumps data structure
97  Mumps_struc_pt=new DMUMPS_STRUC_C;
98 
99  // Mumps' hack to indicate that mpi_comm_world is used. Conversion of any
100  // other communicators appears to be non-portable, so we simply
101  // issue a warning if we later discover that oomph-lib's communicator
102  // is not MPI_COMM_WORLD
103  Mumps_struc_pt->comm_fortran = -987654;
104 
105  // Root processor participates in solution
106  Mumps_struc_pt->par = 1;
107 
108  // Matrix is unsymmetric
109  Mumps_struc_pt->sym = 0;
110 
111  // First call does the initialise phase
112  Mumps_struc_pt->job = -1;
113 
114  // Call c interface to double precision mumps to initialise
115  dmumps_c(Mumps_struc_pt);
117 
118  //Output stream for global info on host. Negative value suppresses printing
119  Mumps_struc_pt->icntl[2]=-1;
120 
121  //Only show error messages and stats
122  Mumps_struc_pt->icntl[3]=2;// 4 for full doc; creates huge amount of output
123 
124  //Write matrix
125  // sprintf(Mumps_struc_pt->write_problem,"/work/e173/e173/mheil/matrix");
126 
127  // Assembled matrix (rather than element-by_element)
128  Mumps_struc_pt->icntl[4]=0;
129 
130  // Distributed problem with user-specified distribution
131  Mumps_struc_pt->icntl[17]=3;
132 
133  // Dense RHS
134  Mumps_struc_pt->icntl[19]=0;
135 
136  // Non-distributed solution
137  Mumps_struc_pt->icntl[20]=0;
138  }
139 
140 
141 
142 //=============================================================================
143 /// Shutdown mumps
144 //=============================================================================
146  {
148  {
149  // Shut down flag
150  Mumps_struc_pt->job = -2;
151 
152  // Call c interface to double precision mumps to shut down
153  dmumps_c(Mumps_struc_pt);
154 
155  // Cleanup
156  delete Mumps_struc_pt;
157  Mumps_struc_pt=0;
158 
159  Mumps_is_initialised=false;
160 
161  // Cleanup storage
162  Irn_loc.clear();
163  Jcn_loc.clear();
164  A_loc.clear();
165 
166  }
167 
168  }
169 
170 
171 //=============================================================================
172 /// Destructor: Shutdown mumps
173 //=============================================================================
175  {
176  shutdown_mumps();
177  }
178 
179 
180 //=============================================================================
181 /// LU decompose the matrix addressed by matrix_pt using
182 /// mumps. The resulting matrix factors are stored internally.
183 /// Note: if Delete_matrix_data is true the function
184 /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
185 //=============================================================================
187 {
188 
189 
190  // Initialise timer
191  double t_start = TimingHelpers::timer();
192 
193  // set the distribution
194  DistributableLinearAlgebraObject* dist_matrix_pt=
195  dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
196  if (dist_matrix_pt)
197  {
198  // the solver has the same distribution as the matrix if possible
199  this->build_distribution(dist_matrix_pt->distribution_pt());
200  }
201  else
202  {
203  throw OomphLibError("Matrix must be distributable",
204  OOMPH_CURRENT_FUNCTION,
205  OOMPH_EXCEPTION_LOCATION);
206  }
207 
208 
209  //Check that we have a square matrix
210 #ifdef PARANOID
211  int n = matrix_pt->nrow();
212  int m = matrix_pt->ncol();
213  if(n != m)
214  {
215  std::ostringstream error_message_stream;
216  error_message_stream << "Can only solve for square matrices\n"
217  << "N, M " << n << " " << m << std::endl;
218 
219  throw OomphLibError(error_message_stream.str(),
220  OOMPH_CURRENT_FUNCTION,
221  OOMPH_EXCEPTION_LOCATION);
222  }
224  {
225  if (this->distribution_pt()->communicator_pt()->mpi_comm()!=
226  MPI_COMM_WORLD)
227  {
228  std::ostringstream error_message_stream;
229  error_message_stream
230  << "Warning: Mumps wrapper assumes that communicator is MPI_COMM_WORLD\n"
231  << " which is not the case, so mumps may die...\n"
232  << " If it does initialise oomph-lib's mpi without requesting\n"
233  << " the communicator to be a duplicate of MPI_COMM_WORLD\n"
234  << " (done via an optional boolean to MPI_Helpers::init(...)\n\n"
235  << " (You can suppress this warning by recompiling without\n"
236  << " paranoia or calling \n\n"
237  << " MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
238  << " \n";
239  OomphLibWarning(error_message_stream.str(),
240  "MumpsSolver::factorise()",
241  OOMPH_EXCEPTION_LOCATION);
242  }
243  }
244 
245 #endif
246 
247  // Initialise
249  {
250  shutdown_mumps();
251  }
253 
254 
255  // Doc stats?
256  if (Doc_stats)
257  {
258  //Output stream for global info on host. Negative value suppresses printing
259  Mumps_struc_pt->icntl[2]=6;
260  }
261 
262  // number of processors
263  unsigned nproc = 1;
264  if (dist_matrix_pt != 0)
265  {
266  nproc = dist_matrix_pt->distribution_pt()->communicator_pt()->nproc();
267  }
268 
269  // Is it a CRDoubleMatrix?
270  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
271  if (cr_matrix_pt != 0)
272  {
273 #ifdef PARANOID
274  // paranoid check that the matrix has been set up
275  if (!cr_matrix_pt->built())
276  {
277  throw OomphLibError
278  ("To apply MumpsSolver to a CRDoubleMatrix - it must be built",
279  OOMPH_CURRENT_FUNCTION,OOMPH_EXCEPTION_LOCATION);
280  }
281 #endif
282 
283  // if the matrix is distributed then set up solver
284  if ((nproc==1)||(cr_matrix_pt->distributed()))
285  {
286  double t_start_copy = TimingHelpers::timer();
287 
288  // Find the number of rows and non-zero entries in the matrix
289  const int nnz_loc = int(cr_matrix_pt->nnz());
290  const int n = matrix_pt->nrow();
291 
292  // Create mumps storage
293 
294  // Create vector for row numbers
295  Irn_loc.resize(nnz_loc);
296 
297  // Create vector for column numbers
298  Jcn_loc.resize(nnz_loc);
299 
300  // Vector for entries
301  A_loc.resize(nnz_loc);
302 
303  // First row held on this processor
304  int first_row = cr_matrix_pt->first_row();
305 
306  // Copy into coordinate storage scheme using pointer arithmetic
307  double* matrix_value_pt = cr_matrix_pt->value();
308  int* matrix_index_pt = cr_matrix_pt->column_index();
309  int* matrix_start_pt = cr_matrix_pt->row_start();
310  int i_row=0;
311  for(int count=0;count<nnz_loc;count++)
312  {
313  A_loc[count] = matrix_value_pt[count];
314  Jcn_loc[count] = matrix_index_pt[count]+1;
315  if (count<matrix_start_pt[i_row+1])
316  {
317  Irn_loc[count] = first_row+i_row+1;
318  }
319  else
320  {
321  i_row++;
322  Irn_loc[count] = first_row+i_row+1;
323  }
324  }
325 
326  // Now delete the matrix if we are allowed
327  if (Delete_matrix_data==true)
328  {
329  cr_matrix_pt->clear();
330  }
331 
332  if ((Doc_time) &&
333  (this->distribution_pt()->communicator_pt()->my_rank()==0))
334  {
335  double t_end_copy = TimingHelpers::timer();
336  oomph_info
337  << "Time for copying matrix into MumpsSolver data structure [sec] : "
338  << t_end_copy-t_start_copy << std::endl;
339  }
340 
341  // Call mumps factorisation
342  //-------------------------
343 
344  //Specify size of system
345  Mumps_struc_pt->n=n;
346 
347  // Number of nonzeroes
348  Mumps_struc_pt->nz_loc=nnz_loc;
349 
350  //The entries
351  Mumps_struc_pt->irn_loc=&Irn_loc[0];
352  Mumps_struc_pt->jcn_loc=&Jcn_loc[0];
353  Mumps_struc_pt->a_loc=&A_loc[0];
354 
355  double t_start_analyse = TimingHelpers::timer();
356 
357  // Do analysis
358  Mumps_struc_pt->job = 1;
359  dmumps_c(Mumps_struc_pt);
360 
361 
362  if ((Doc_time) &&
363  (this->distribution_pt()->communicator_pt()->my_rank()==0))
364  {
365  double t_end_analyse = TimingHelpers::timer();
366  oomph_info
367  << "Time for mumps analysis stage in MumpsSolver [sec] : "
368  << t_end_analyse-t_start_analyse << std::endl;
369  }
370 
371 
372  int my_rank=this->distribution_pt()->communicator_pt()->my_rank();
373 
374  // Document estimate for size of workspace
375  if (my_rank==0)
376  {
378  {
379  oomph_info << "Estimated max. workspace in MB: "
380  << Mumps_struc_pt->infog[25] << std::endl;
381  }
382  }
383 
384  double t_start_factor = TimingHelpers::timer();
385 
386  // Loop until successfully factorised
387  bool factorised=false;
388  while (!factorised)
389  {
390 
391  // Set workspace to multiple of that -- ought to be "significantly
392  // larger than infog(26)", according to the manual :(
394  Mumps_struc_pt->infog[25];
395 
397  {
398  oomph_info << "Attempting factorisation with workspace estimate: "
399  << Mumps_struc_pt->icntl[22] << " MB\n";
400  oomph_info << "corresponding to Workspace_scaling_factor = "
401  << Workspace_scaling_factor << "\n";
402  }
403 
404  // Do factorisation
405  Mumps_struc_pt->job = 2;
406  dmumps_c(Mumps_struc_pt);
407 
408  // Check for error
409  if (Mumps_struc_pt->infog[0]!=0)
410  {
412  {
413  oomph_info << "Error during mumps factorisation!\n";
414  oomph_info << "Error codes: "
415  << Mumps_struc_pt->info[0] << " "
416  << Mumps_struc_pt->info[1] << std::endl;
417  }
418 
419  // Increase scaling factor for workspace and run again
421 
423  {
424  oomph_info << "Increasing workspace_scaling_factor to "
425  << Workspace_scaling_factor << std::endl;
426  }
427  }
428  else
429  {
430  factorised=true;
431  }
432  }
433 
434 
435  if ((Doc_time) &&
436  (this->distribution_pt()->communicator_pt()->my_rank()==0))
437  {
438  double t_end_factor = TimingHelpers::timer();
439  oomph_info
440  << "Time for actual mumps factorisation in MumpsSolver [sec] : "
441  << t_end_factor-t_start_factor<< std::endl;
442  }
443  }
444  // else the CRDoubleMatrix is not distributed
445  else
446  {
447  std::ostringstream error_message_stream;
448  error_message_stream << "MumpsSolver only works for a "
449  << " distributed CRDoubleMatrix\n";
450  throw OomphLibError(error_message_stream.str(),
451  OOMPH_CURRENT_FUNCTION,
452  OOMPH_EXCEPTION_LOCATION);
453 
454  }
455  }
456  // Otherwise throw an error
457  else
458  {
459  std::ostringstream error_message_stream;
460  error_message_stream << "MumpsSolver implemented only for "
461  << "distributed CRDoubleMatrix. \n";
462  throw OomphLibError(error_message_stream.str(),
463  OOMPH_CURRENT_FUNCTION,
464  OOMPH_EXCEPTION_LOCATION);
465  }
466 
467 
468  if ((Doc_time) && (this->distribution_pt()->communicator_pt()->my_rank()==0))
469  {
470 
471  double t_end = TimingHelpers::timer();
472  oomph_info << "Time for MumpsSolver factorisation [sec] : "
473  << t_end-t_start << std::endl;
474  }
475 
476  // Switch off docing again by setting output stream for global info on
477  // to negative number
478  Mumps_struc_pt->icntl[2]=-1;
479 
480 }
481 
482 //=============================================================================
483 /// Do the backsubstitution for mumps solver.
484 /// This does not make any assumption about the distribution of the
485 /// vectors
486 //=============================================================================
488  DoubleVector &result)
489  {
490 
491  double t_start = TimingHelpers::timer();
492 
493 #ifdef PARANOID
495  {
496  if (this->distribution_pt()->communicator_pt()->mpi_comm()!=MPI_COMM_WORLD)
497  {
498  std::ostringstream error_message_stream;
499  error_message_stream
500  << "Warning: Mumps wrapper assumes that communicator is MPI_COMM_WORLD\n"
501  << " which is not the case, so mumps may die...\n"
502  << " If it does initialise oomph-lib's mpi without requesting\n"
503  << " the communicator to be a duplicate of MPI_COMM_WORLD\n"
504  << " (done via an optional boolean to MPI_Helpers::init(...)\n\n"
505  << " (You can suppress this warning by recompiling without\n"
506  << " paranoia or calling \n\n"
507  << " MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
508  << " \n";
509  OomphLibWarning(error_message_stream.str(),
510  "MumpsSolver::backsub()",
511  OOMPH_EXCEPTION_LOCATION);
512  }
513  }
514 
515  // Initialised?
517  {
518  std::ostringstream error_message_stream;
519  error_message_stream
520  << "Mumps has not been initialised.";
521  throw OomphLibError(error_message_stream.str(),
522  OOMPH_CURRENT_FUNCTION,
523  OOMPH_EXCEPTION_LOCATION);
524  }
525 
526  // check that the rhs vector is setup
527  if (!rhs.distribution_pt()->built())
528  {
529  std::ostringstream error_message_stream;
530  error_message_stream
531  << "The vectors rhs must be setup";
532  throw OomphLibError(error_message_stream.str(),
533  OOMPH_CURRENT_FUNCTION,
534  OOMPH_EXCEPTION_LOCATION);
535  }
536 
537 #endif
538 
539  // Check that the rhs distribution is the same as the distribution as this
540  // solver. If not redistribute and issue a warning
541  LinearAlgebraDistribution rhs_distribution(rhs.distribution_pt());
542  if (!(*rhs.distribution_pt() == *this->distribution_pt()))
543  {
545  {
546  std::ostringstream warning_stream;
547  warning_stream
548  << "The distribution of rhs vector does not match that of the solver.\n";
549  warning_stream
550  << "The rhs may have to be redistributed but we're not doing this because\n"
551  << "I'm no longer convinced it's necessary. Keep an eye on this...\n";
552  warning_stream
553  << "To remove this warning you can either:\n"
554  << " i) Ensure that the rhs vector has the correct distribution\n"
555  << " before calling the resolve() function\n"
556  << "or ii) Set the flag \n"
557  << " MumpsSolver::Suppress_incorrect_rhs_distribution_warning_in_resolve\n"
558  << " to be true\n\n";
559 
560  OomphLibWarning(warning_stream.str(),
561  "MumpsSolver::resolve()",
562  OOMPH_EXCEPTION_LOCATION);
563  }
564 
565  // //Have to cast away const-ness (which tells us that we shouldn't really
566  // //be doing this!)
567  // const_cast<DoubleVector&>(rhs).redistribute(this->distribution_pt());
568  }
569 
570 
571 
572 #ifdef PARANOID
573  // if the result vector is setup then check it has the same distribution
574  // as the rhs
575  if (result.distribution_built())
576  {
577  if (!(*result.distribution_pt() == *rhs.distribution_pt()))
578  {
579  std::ostringstream error_message_stream;
580  error_message_stream
581  << "The result vector distribution has been setup; it must have the "
582  << "same distribution as the rhs vector.";
583  throw OomphLibError(error_message_stream.str(),
584  OOMPH_CURRENT_FUNCTION,
585  OOMPH_EXCEPTION_LOCATION);
586  }
587  }
588 #endif
589 
590 
591  // Doc stats?
592  if (Doc_stats)
593  {
594  //Output stream for global info on host. Negative value suppresses printing
595  Mumps_struc_pt->icntl[2]=6;
596  }
597 
598  // number of DOFs
599  int ndof = this->distribution_pt()->nrow();
600 
601  // Make backup to avoid over-writing
602  DoubleVector tmp_rhs;
603  tmp_rhs=rhs;
604 
605  // Now turn this into a global (non-distributed) vector
606  // because that's what mumps needs
607 
608  // Make a global distribution (i.e. one that isn't distributed)
609  LinearAlgebraDistribution global_distribution(
610  this->distribution_pt()->communicator_pt(),ndof,false);
611 
612  // Redistribute the tmp_rhs vector with this distribution -- it's
613  // now "global", as required for mumps
614  tmp_rhs.redistribute(&global_distribution);
615 
616  // Do the backsubsitution phase -- overwrites the tmp_rhs vector with the
617  // solution
618  Mumps_struc_pt->rhs=&tmp_rhs[0];
619  Mumps_struc_pt->job=3;
620  dmumps_c(Mumps_struc_pt);
621 
622  // Copy back
623  unsigned n=Mumps_struc_pt->n;
624  for (unsigned j=0;j<n;j++)
625  {
626  tmp_rhs[j]=Mumps_struc_pt->rhs[j];
627  }
628 
629  // Broadcast the result which is only held on root
630  MPI_Bcast(&tmp_rhs[0],ndof,MPI_DOUBLE,0,
631  this->distribution_pt()->communicator_pt()->mpi_comm());
632 
633  // If the result vector is distributed, re-distribute the
634  // non-distributed tmp_rhs vector to match
635  if (result.built())
636  {
637  tmp_rhs.redistribute(result.distribution_pt());
638  }
639  else
640  {
641  tmp_rhs.redistribute(this->distribution_pt());
642  }
643 
644  // Now copy the tmp_rhs vector into the (matching) result
645  result = tmp_rhs;
646 
647  if ((Doc_time) && (this->distribution_pt()->communicator_pt()->my_rank()==0))
648  {
649  double t_end = TimingHelpers::timer();
650  oomph_info << "Time for MumpsSolver backsub [sec] : "
651  << t_end-t_start << std::endl;
652  }
653 
654  // Switch off docing again by setting output stream for global info on
655  // to negative number
656  Mumps_struc_pt->icntl[2]=-1;
657 
658 }
659 
660 //=========================================================================
661 /// Clean up the memory allocated for the MumpsSolver solver
662 //=========================================================================
664 {
665  shutdown_mumps();
666 }
667 
668 
669 //=========================================================================
670 /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
671 /// vector and returns the solution of the linear system. Problem pointer
672 /// defaults to NULL and can be omitted. The function returns the global
673 /// result vector. Matrix must be CRDoubleMatrix.
674 /// Note: if Delete_matrix_data is true the function
675 /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
676 //=========================================================================
677 void MumpsSolver::solve(DoubleMatrixBase* const &matrix_pt,
678  const DoubleVector &rhs,
679  DoubleVector &result)
680 {
681 #ifdef PARANOID
683  {
684  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
685  ->distribution_pt()->communicator_pt()->mpi_comm()!=
686  MPI_COMM_WORLD)
687  {
688  std::ostringstream error_message_stream;
689  error_message_stream
690  << "Warning: Mumps wrapper assumes that communicator is MPI_COMM_WORLD\n"
691  << " which is not the case, so mumps may die...\n"
692  << " If it does initialise oomph-lib's mpi without requesting\n"
693  << " the communicator to be a duplicate of MPI_COMM_WORLD\n"
694  << " (done via an optional boolean to MPI_Helpers::init(...)\n\n"
695  << " (You can suppress this warning by recompiling without\n"
696  << " paranoia or calling \n\n"
697  << " MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
698  << " \n";
699  OomphLibWarning(error_message_stream.str(),
700  "MumpsSolver::solve()",
701  OOMPH_EXCEPTION_LOCATION);
702  }
703  }
704 #endif
705 
706  // Initialise timer
707  double t_start = TimingHelpers::timer();
708 
709 #ifdef PARANOID
710  // check that the rhs vector is setup
711  if (!rhs.distribution_pt()->built())
712  {
713  std::ostringstream error_message_stream;
714  error_message_stream
715  << "The vectors rhs must be setup";
716  throw OomphLibError(error_message_stream.str(),
717  OOMPH_CURRENT_FUNCTION,
718  OOMPH_EXCEPTION_LOCATION);
719  }
720 
721  // check that the matrix is square
722  if (matrix_pt->nrow() != matrix_pt->ncol())
723  {
724  std::ostringstream error_message_stream;
725  error_message_stream
726  << "The matrix at matrix_pt must be square.";
727  throw OomphLibError(error_message_stream.str(),
728  OOMPH_CURRENT_FUNCTION,
729  OOMPH_EXCEPTION_LOCATION);
730  }
731 
732  // check that the matrix and the rhs vector have the same nrow()
733  if (matrix_pt->nrow() != rhs.nrow())
734  {
735  std::ostringstream error_message_stream;
736  error_message_stream
737  << "The matrix and the rhs vector must have the same number of rows.";
738  throw OomphLibError(error_message_stream.str(),
739  OOMPH_CURRENT_FUNCTION,
740  OOMPH_EXCEPTION_LOCATION);
741  }
742 
743 
744  // if the matrix is distributable then should have the same distribution
745  // as the rhs vector
746  DistributableLinearAlgebraObject* ddist_matrix_pt =
747  dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
748  if (ddist_matrix_pt != 0)
749  {
750  if (!(*ddist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
751  {
752  std::ostringstream error_message_stream;
753  error_message_stream
754  << "The matrix matrix_pt must have the same distribution as the "
755  << "rhs vector.";
756  throw OomphLibError(error_message_stream.str(),
757  OOMPH_CURRENT_FUNCTION,
758  OOMPH_EXCEPTION_LOCATION);
759  }
760  }
761  // if the matrix is not distributable then it the rhs vector should not be
762  // distributed
763  else
764  {
765  if (rhs.distribution_pt()->distributed())
766  {
767  std::ostringstream error_message_stream;
768  error_message_stream
769  << "The matrix (matrix_pt) is not distributable and therefore the rhs"
770  << " vector must not be distributed";
771  throw OomphLibError(error_message_stream.str(),
772  OOMPH_CURRENT_FUNCTION,
773  OOMPH_EXCEPTION_LOCATION);
774  }
775  }
776  // if the result vector is setup then check it has the same distribution
777  // as the rhs
778  if (result.built())
779  {
780  if (!(*result.distribution_pt() == *rhs.distribution_pt()))
781  {
782  std::ostringstream error_message_stream;
783  error_message_stream
784  << "The result vector distribution has been setup; it must have the "
785  << "same distribution as the rhs vector.";
786  throw OomphLibError(error_message_stream.str(),
787  OOMPH_CURRENT_FUNCTION,
788  OOMPH_EXCEPTION_LOCATION);
789  }
790  }
791 
792 #endif
793 
794 
795 
796  // set the distribution
797  DistributableLinearAlgebraObject* dist_matrix_pt=
798  dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
799  if (dist_matrix_pt)
800  {
801  // the solver has the same distribution as the matrix if possible
802  this->build_distribution(dist_matrix_pt->distribution_pt());
803  }
804  else
805  {
806  // the solver has the same distribution as the RHS
807  this->build_distribution(rhs.distribution_pt());
808  }
809 
810  // Factorise the matrix
811  factorise(matrix_pt);
812 
813  //Now do the back solve
814  backsub(rhs,result);
815 
816  // Doc time for solve
817  double t_end = TimingHelpers::timer();
818  Solution_time = t_end-t_start;
819 
820  if ((Doc_time) && (this->distribution_pt()->communicator_pt()->my_rank()==0))
821  {
822 
823  oomph_info << "Time for MumpsSolver solve [sec] : "
824  << t_end-t_start << std::endl;
825  }
826 
827  // Switch off docing again by setting output stream for global info on
828  // to negative number
829  Mumps_struc_pt->icntl[2]=-1;
830 
831  // If we are not storing the solver data for resolves, delete it
832  if (!Enable_resolve)
833  {
834  clean_up_memory();
835  }
836 
837 
838  }
839 
840 //==================================================================
841 /// Solver: Takes pointer to problem and returns the results Vector
842 /// which contains the solution of the linear system defined by
843 /// the problem's fully assembled Jacobian and residual Vector.
844 //==================================================================
845 void MumpsSolver::solve(Problem* const &problem_pt, DoubleVector &result)
846 {
847 
848 #ifdef PARANOID
850  {
851  if (problem_pt->communicator_pt()->mpi_comm()!=MPI_COMM_WORLD)
852  {
853  std::ostringstream error_message_stream;
854  error_message_stream
855  << "Warning: Mumps wrapper assumes that communicator is MPI_COMM_WORLD\n"
856  << " which is not the case, so mumps may die...\n"
857  << " If it does initialise oomph-lib's mpi without requesting\n"
858  << " the communicator to be a duplicate of MPI_COMM_WORLD\n"
859  << " (done via an optional boolean to MPI_Helpers::init(...)\n\n"
860  << " (You can suppress this warning by recompiling without\n"
861  << " paranoia or calling \n\n"
862  << " MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
863  << " \n";
864  OomphLibWarning(error_message_stream.str(),
865  "MumpsSolver::solve()",
866  OOMPH_EXCEPTION_LOCATION);
867  }
868  }
869 #endif
870 
871  // Initialise timer
872  double t_start = TimingHelpers::timer();
873 
874  // number of dofs
875  unsigned n_dof = problem_pt->ndof();
876 
877  // Set the distribution for the solver.
878  this->distribution_pt()->build(problem_pt->communicator_pt(),n_dof);
879 
880  // Take a copy of Delete_matrix_data
881  bool copy_of_Delete_matrix_data = Delete_matrix_data;
882 
883  // Set Delete_matrix to true
884  Delete_matrix_data = true;
885 
886  // Initialise timer
887  t_start = TimingHelpers::timer();
888 
889  // Storage for the distributed residuals vector
890  DoubleVector residuals(this->distribution_pt(),0.0);
891 
892  // Get the sparse jacobian and residuals of the problem
893  CRDoubleMatrix jacobian(this->distribution_pt());
894  problem_pt->get_jacobian(residuals, jacobian);
895 
896  // Doc time for setup
897  double t_end = TimingHelpers::timer();
898  Jacobian_setup_time = t_end-t_start;
899  int my_rank = this->distribution_pt()->communicator_pt()->my_rank();
900  if ((Doc_time) && (my_rank==0))
901  {
902  oomph_info << "Time to set up CRDoubleMatrix Jacobian [sec] : "
903  << Jacobian_setup_time << std::endl;
904  }
905 
906 
907  //Now call the linear algebra solve, if desired
908  if(!Suppress_solve)
909  {
910  //If the distribution of the result has been build and
911  //does not match that of
912  //the solver then redistribute before the solve and return
913  //to the incoming distribution afterwards.
914  if((result.built()) &&
915  (!(*result.distribution_pt() == *this->distribution_pt())))
916  {
918  temp_global_dist(result.distribution_pt());
919  result.build(this->distribution_pt(),0.0);
920  solve(&jacobian,residuals,result);
921  result.redistribute(&temp_global_dist);
922  }
923  else
924  {
925  solve(&jacobian,residuals,result);
926  }
927  }
928 
929  // Set Delete_matrix back to original value
930  Delete_matrix_data = copy_of_Delete_matrix_data;
931 
932  // Finalise/doc timings
933  if ((Doc_time) && (my_rank==0))
934  {
935  double t_end = TimingHelpers::timer();
936  oomph_info << "Total time for MumpsSolver " << "(np="
937  << this->distribution_pt()->communicator_pt()->nproc()
938  << ",N=" << problem_pt->ndof()
939  <<") [sec] : " << t_end-t_start << std::endl;
940  }
941 
942 }
943 
944 
945 //===============================================================
946 /// Resolve the system defined by the last assembled jacobian
947 /// and the specified rhs vector if resolve has been enabled.
948 /// Note: returns the global result Vector.
949 //===============================================================
951  {
952 
953 #ifdef PARANOID
955  {
956  if (this->distribution_pt()->communicator_pt()->mpi_comm()!=MPI_COMM_WORLD)
957  {
958  std::ostringstream error_message_stream;
959  error_message_stream
960  << "Warning: Mumps wrapper assumes communicator is MPI_COMM_WORLD\n"
961  << " which is not the case, so mumps may die...\n"
962  << " If it does initialise oomph-lib's mpi without requesting\n"
963  << " the communicator to be a duplicate of MPI_COMM_WORLD\n"
964  << " (done via an optional boolean to MPI_Helpers::init(...)\n\n"
965  << " (You can suppress this warning by recompiling without\n"
966  << " paranoia or calling\n\n"
967  << " MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
968  << " \n";
969  OomphLibWarning(error_message_stream.str(),
970  "MumpsSolver::resolve()",
971  OOMPH_EXCEPTION_LOCATION);
972  }
973  }
974 #endif
975 
976  // Store starting time for solve
977  double t_start = TimingHelpers::timer();
978 
979  // Doc stats?
980  if (Doc_stats)
981  {
982  //Output stream for global info on host. Negative value suppresses printing
983  Mumps_struc_pt->icntl[2]=6;
984  }
985 
986  //Now do the back substitution phase
987  backsub(rhs,result);
988 
989  // Doc time for solve
990  double t_end = TimingHelpers::timer();
991  Solution_time = t_end-t_start;
992 
993  // Switch off docing again by setting output stream for global info on
994  // to negative number
995  Mumps_struc_pt->icntl[2]=-1;
996 
997  if ((Doc_time) && (this->distribution_pt()->communicator_pt()->my_rank()==0))
998  {
999  oomph_info << "Time for MumpsSolver solve [sec]: "
1000  << t_end-t_start << std::endl;
1001  }
1002 
1003  }
1004 
1005 
1006 
1007 
1008 }
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1056
bool built() const
double Solution_time
Solution time.
Definition: mumps_solver.h:205
unsigned Workspace_scaling_factor
Definition: mumps_solver.h:223
bool Doc_stats
Set to true for MumpsSolver to output statistics (false by default).
Definition: mumps_solver.h:211
unsigned first_row() const
access function for the first row on this processor
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...
static bool Suppress_incorrect_rhs_distribution_warning_in_resolve
Static flag that determines whether the warning about incorrect distribution of RHSs will be printed ...
Definition: mumps_solver.h:67
static int Default_workspace_scaling_factor
Default factor for workspace – static so it can be overwritten globally.
Definition: mumps_solver.h:191
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
void initialise_mumps()
Initialise instance of mumps data structure.
Definition: mumps_solver.cc:93
Vector< double > A_loc
Definition: mumps_solver.h:240
void build(const OomphCommunicator *const comm_pt, const unsigned &first_row, const unsigned &nrow_local, const unsigned &nrow=0)
Sets the distribution. Takes first_row, nrow_local and nrow as arguments. If nrow is not provided or ...
bool distributed() const
distribution is serial or distributed
Vector< int > Jcn_loc
Definition: mumps_solver.h:237
MumpsSolver()
Constructor: Call setup.
Definition: mumps_solver.cc:77
OomphInfo oomph_info
double Jacobian_setup_time
Jacobian setup time.
Definition: mumps_solver.h:202
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
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
bool Delete_matrix_data
Delete_matrix_data flag. MumpsSolver needs its own copy of the input matrix, therefore a copy must be...
Definition: mumps_solver.h:231
void backsub(const DoubleVector &rhs, DoubleVector &result)
Do the backsubstitution for mumps solver Note: returns the global result Vector.
DMUMPS_STRUC_C * Mumps_struc_pt
Pointer to MUMPS struct that contains the solver data.
Definition: mumps_solver.h:243
Vector< int > Irn_loc
Vector for row numbers.
Definition: mumps_solver.h:234
bool Mumps_is_initialised
Has mumps been initialised?
Definition: mumps_solver.h:220
void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled Jacobian and the specified rhs vector if resolve has...
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1068
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void clear()
clear
Definition: matrices.cc:1677
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...
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:84
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1574
unsigned nrow() const
access function to the number of global rows.
void factorise(DoubleMatrixBase *const &matrix_pt)
Do the factorisation stage Note: if Delete_matrix_data is true the function matrix_pt->clean_up_memor...
bool Suppress_mumps_info_during_solve
Boolean to suppress info being printed to screen during MUMPS solve.
Definition: mumps_solver.h:217
bool Suppress_warning_about_MPI_COMM_WORLD
Boolean to suppress warning about communicator not equal to MPI_COMM_WORLD.
Definition: mumps_solver.h:214
double timer()
returns the time in seconds after some point in past
bool built() const
access function to the Built flag - indicates whether the matrix has been build - i...
Definition: matrices.h:1177
void clean_up_memory()
Clean up the memory allocated by the mumps solver.
~MumpsSolver()
Destructor: Cleanup.
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
Definition: problem.cc:3851
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results Vector which contains the solution of the li...
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
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.
bool Suppress_solve
Suppress solve?
Definition: mumps_solver.h:208
void shutdown_mumps()
Shutdown mumps.