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