Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_MUMPS_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Amesos2: Templated Direct Sparse Solver Package
4//
5// Copyright 2011 NTESS and the Amesos2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
18#ifndef AMESOS2_MUMPS_DEF_HPP
19#define AMESOS2_MUMPS_DEF_HPP
20
21#include <Teuchos_Tuple.hpp>
22#include <Teuchos_ParameterList.hpp>
23#include <Teuchos_StandardParameterEntryValidators.hpp>
24
25#ifdef HAVE_MPI
26#include <Teuchos_DefaultMpiComm.hpp>
27#endif
28
29#include <limits>
30
31#include "Amesos2_SolverCore_def.hpp"
32#include "Amesos2_MUMPS_decl.hpp"
33
34namespace Amesos2
35{
36
37
38 template <class Matrix, class Vector>
39 MUMPS<Matrix,Vector>::MUMPS(
40 Teuchos::RCP<const Matrix> A,
41 Teuchos::RCP<Vector> X,
42 Teuchos::RCP<const Vector> B )
43 : SolverCore<Amesos2::MUMPS,Matrix,Vector>(A, X, B)
44 , schur_size(0)
45 , schur_out_ptr(nullptr)
46 , only_forward_solve(false)
47 , only_backward_solve(false)
48 , is_contiguous_(true)
49 {
50
51 typedef FunctionMap<MUMPS,scalar_type> function_map;
52
53 MUMPS_MATRIX_LOAD = false;
54 MUMPS_STRUCT = false;
55 MUMPS_MATRIX_LOAD_PREORDERING = false;
56
57 #ifdef HAVE_MPI
58 using Teuchos::Comm;
59 using Teuchos::MpiComm;
60 using Teuchos::RCP;
61 using Teuchos::rcp;
62 using Teuchos::rcp_dynamic_cast;
63
64 //Comm
65 mumps_par.comm_fortran = -987654;
66 RCP<const Comm<int> > matComm = this->matrixA_->getComm();
67 //Add Exception Checking
68 TEUCHOS_TEST_FOR_EXCEPTION(
69 matComm.is_null(), std::logic_error, "Amesos2::MUMPS:: matComm is null");
70 RCP<const MpiComm<int> > matMpiComm =
71 rcp_dynamic_cast<const MpiComm<int> >(matComm);
72 //Add Exception Checking
73 TEUCHOS_TEST_FOR_EXCEPTION(
74 matMpiComm.is_null(),
75 std::logic_error, "Amesos2::MUMPS:: matMPIComm is null");
76 TEUCHOS_TEST_FOR_EXCEPTION(
77 matMpiComm->getRawMpiComm().is_null(),
78 std::logic_error, "Amesos2::MUMPS:: rawMPIComm is null");
79 MPI_Comm rawMpiComm = (* (matMpiComm->getRawMpiComm()) )();
80 mumps_par.comm_fortran = (int) MPI_Comm_c2f(rawMpiComm);
81 #endif
82
83 mumps_par.job = -1;
84 mumps_par.par = 1;
85 mumps_par.sym = 0;
86
87 function_map::mumps_c(&(mumps_par)); //init
88 MUMPS_ERROR();
89
90 mumps_par.n = this->globalNumCols_;
91
92 /*Default parameters*/
93 mumps_par.icntl[0] = -1; // Turn off error messages
94 mumps_par.icntl[1] = -1; // Turn off diagnositic printing
95 mumps_par.icntl[2] = -1; // Turn off global information messages
96 mumps_par.icntl[3] = 1; // No messages
97 mumps_par.icntl[4] = 0; // Matrix given in assembled (Triplet) form
98 mumps_par.icntl[5] = 7; // Choose column permuation automatically
99 mumps_par.icntl[6] = 7; // Choose ordering methods automatically
100 mumps_par.icntl[7] = 7; // Choose scaling automatically
101 mumps_par.icntl[8] = 1; // Compuate Ax = b;
102 mumps_par.icntl[9] = 0; // iterative refinement
103 mumps_par.icntl[10] = 0; //Do not collect statistics
104 mumps_par.icntl[11] = 0; // Automatic choice of ordering strategy
105 mumps_par.icntl[12] = 0; //Use ScaLAPACK for root node
106 mumps_par.icntl[13] = 20; // Increase memory allocation 20% at a time
107 mumps_par.icntl[17] = 0;
108 mumps_par.icntl[18] = 0; // do not provide back the Schur complement
109 mumps_par.icntl[19] = 0; // RHS is given in dense form
110 mumps_par.icntl[20] = 0; //RHS is in dense form
111 mumps_par.icntl[21] = 0; // Do all compuations in-core
112 mumps_par.icntl[22] = 0; // No max MB for work
113 mumps_par.icntl[23] = 0; // Do not perform null pivot detection
114 mumps_par.icntl[24] = 0; // No null space basis compuation
115 mumps_par.icntl[25] = 0; // Do not condense/reduce Schur RHS
116 mumps_par.icntl[27] = 1; // sequential analysis
117 mumps_par.icntl[28] = 0; //
118 mumps_par.icntl[29] = 0; //
119 mumps_par.icntl[30] = 0; //
120 mumps_par.icntl[31] = 0;
121 mumps_par.icntl[32] = 0;
122 }
123
124 template <class Matrix, class Vector>
125 MUMPS<Matrix,Vector>::~MUMPS( )
126 {
127 /* Clean up the struc*/
128 typedef FunctionMap<MUMPS,scalar_type> function_map;
129
130 if(MUMPS_STRUCT == true)
131 {
132 free(mumps_par.a);
133 free(mumps_par.jcn);
134 free(mumps_par.irn);
135 }
136 mumps_par.job = -2;
137 if (this->rank_ < this->nprocs_) {
138 function_map::mumps_c(&(mumps_par));
139 }
140
141 }
142
143 template<class Matrix, class Vector>
144 int
145 MUMPS<Matrix,Vector>::preOrdering_impl()
146 {
147 /* TODO: Define what it means for MUMPS
148 */
149 #ifdef HAVE_AMESOS2_TIMERS
150 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
151 #endif
152
153 return(0);
154 }//end preOrdering_impl()
155
156 template <class Matrix, class Vector>
157 int
158 MUMPS<Matrix,Vector>::symbolicFactorization_impl()
159 {
160 #ifdef HAVE_AMESOS2_TIMERS
161 Teuchos::TimeMonitor symFactTime( this->timers_.symFactTime_ );
162 #endif
163
164 typedef FunctionMap<MUMPS,scalar_type> function_map;
165 if ( this->globalNumRows_ > 0 ) {
166 mumps_par.par = 1;
167 mumps_par.job = 1; // sym factor
168 function_map::mumps_c(&(mumps_par));
169 MUMPS_ERROR();
170 }
171 return(0);
172 }//end symblicFactortion_impl()
173
174
175 template <class Matrix, class Vector>
176 int
177 MUMPS<Matrix,Vector>::numericFactorization_impl()
178 {
179 #ifdef HAVE_AMESOS2_TIMERS
180 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
181 #endif
182
183 typedef FunctionMap<MUMPS,scalar_type> function_map;
184 if ( this->globalNumRows_ > 0 ) {
185 mumps_par.job = 2;
186 function_map::mumps_c(&(mumps_par));
187 MUMPS_ERROR();
188
189 if (schur_size > 0 && schur_out_ptr != nullptr) {
190 // copy schur out if the output pointer is valid (assuming enough space has been allocated)
191 // MUMPS returns the Schur complement in row-major
192 // So, we transpose it to store in column-major
193 for (size_t i = 0; i < schur_size; i++) {
194 for (size_t j = 0; j < schur_size; j++) {
195 schur_out_ptr[j+i*schur_size] = schur_out(i+j*schur_size);
196 }
197 }
198 }
199 }
200 return(0);
201 }//end numericFactorization_impl()
202
203
204 template <class Matrix, class Vector>
205 int
206 MUMPS<Matrix,Vector>::solve_impl(
207 const Teuchos::Ptr< MultiVecAdapter<Vector> > X,
208 const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
209 {
210 typedef FunctionMap<MUMPS,scalar_type> function_map;
211
212 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
213 const size_t nrhs = X->getGlobalNumVectors();
214 const size_t val_store_size = Teuchos::as<size_t>(ld_rhs * nrhs);
215
216 const bool initialize_data = true;
217 const bool do_not_initialize_data = false;
218 {
219 #ifdef HAVE_AMESOS2_TIMERS
220 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
221 #endif
222
223 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
224 host_mv_view>::do_get(initialize_data, B, bvals_,
225 Teuchos::as<size_t>(ld_rhs),
226 (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
227 this->rowIndexBase_);
228 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
229 host_mv_view>::do_get(do_not_initialize_data, X, xvals_,
230 Teuchos::as<size_t>(ld_rhs),
231 (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
232 this->rowIndexBase_);
233
234 // copy b to x to give x to MUMPS
235 Kokkos::deep_copy(xvals_, bvals_);
236 }
237 int ierr = 0; // returned error code
238 if ( this->globalNumRows_ > 0 ) {
239 #ifdef HAVE_AMESOS2_TIMERS
240 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
241 #endif
242
243 mumps_par.nrhs = nrhs;
244 mumps_par.lrhs = mumps_par.n;
245 mumps_par.job = 3;
246 if ( this->root_ )
247 {
248 // MUMPS is called with all MPI ranks
249 // but only root rank has the actual data.
250 // (b has been copied into x)
251 mumps_par.rhs = xvals_.data();
252 }
253
254 if (only_forward_solve || only_backward_solve) {
255 // do only forward or backward solve
256 if (only_forward_solve) {
257 mumps_par.icntl[25] = 1;
258 } else {
259 mumps_par.icntl[25] = 2;
260 }
261
262 // copy rhs to reduced-rhs
263 Kokkos::resize(schur_rhs, schur_size*nrhs);
264 size_t n = this->globalNumRows_;
265 size_t n2 = schur_size;
266 size_t n1 = n - n2;
267 for (size_t i = 0; i < n2; i++) {
268 for (size_t j = 0; j < nrhs; j++) schur_rhs[i+j*schur_size] = xvals_(n1+i, j);
269 }
270 mumps_par.redrhs = schur_rhs.data();
271 } else {
272 // forward solve, followed by backward solve
273 mumps_par.icntl[25] = 0;
274 }
275
276 function_map::mumps_c(&(mumps_par));
277 MUMPS_ERROR();
278 }
279 {
280 #ifdef HAVE_AMESOS2_TIMERS
281 Teuchos::TimeMonitor redistTimer2(this->timers_.vecRedistTime_);
282 #endif
283 if (only_forward_solve) {
284 // copy reduced-sol to sol
285 size_t n = this->globalNumRows_;
286 size_t n2 = schur_size;
287 size_t n1 = n - n2;
288 for (size_t i = 0; i < n2; i++) {
289 for (size_t j = 0; j < nrhs; j++) xvals_(n1+i, j) = schur_rhs[i+j*schur_size];
290 }
291 }
292
293 Util::put_1d_data_helper_kokkos_view<MultiVecAdapter<Vector>,
294 host_mv_view>::do_put(X, xvals_,
295 Teuchos::as<size_t>(ld_rhs),
296 (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
297 this->rowIndexBase_);
298 }
299 // ch: see function loadA_impl()
300 MUMPS_MATRIX_LOAD_PREORDERING = false;
301 return(ierr);
302 }//end solve()
303
304
305 template <class Matrix, class Vector>
306 bool
307 MUMPS<Matrix,Vector>::matrixShapeOK_impl() const
308 {
309 // The MUMPS can only handle square for right now
310 return( this->globalNumRows_ == this->globalNumCols_ );
311 }
312
313
314 template <class Matrix, class Vector>
315 void
316 MUMPS<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
317 {
318 using Teuchos::RCP;
319 using Teuchos::getIntegralValue;
320 using Teuchos::ParameterEntryValidator;
321
322 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
323 /*To Do --- add support for parameters */
324 if(parameterList->isParameter("ICNTL(1)")){
325 mumps_par.icntl[0] = parameterList->get<int>("ICNTL(1)", -1);
326 }
327 if(parameterList->isParameter("ICNTL(2)")){
328 mumps_par.icntl[1] = parameterList->get<int>("ICNTL(2)", -1);
329 }
330 if(parameterList->isParameter("ICNTL(3)")){
331 mumps_par.icntl[2] = parameterList->get<int>("ICNTL(3)", -1);
332 }
333 if(parameterList->isParameter("ICNTL(4)")){
334 mumps_par.icntl[3] = parameterList->get<int>("ICNTL(4)", 1);
335 }
336 if(parameterList->isParameter("ICNTL(6)")){
337 mumps_par.icntl[5] = parameterList->get<int>("ICNTL(6)", 0);
338 }
339 if(parameterList->isParameter("ICNTL(7)")){
340 mumps_par.icntl[6] = parameterList->get<int>("ICNTL(7)", 7);
341 }
342 if(parameterList->isParameter("ICNTL(9)")){
343 mumps_par.icntl[8] = parameterList->get<int>("ICNTL(9)", 1);
344 }
345 if(parameterList->isParameter("ICNTL(11)")){
346 mumps_par.icntl[10] = parameterList->get<int>("ICNTL(11)", 0);
347 }
348 if(parameterList->isParameter("ICNTL(14)")){
349 mumps_par.icntl[13] = parameterList->get<int>("ICNTL(14)", 20);
350 }
351 /* For Partial factorization */
352 if(parameterList->isParameter("PartialFacto")) {
353 int partial_facto = parameterList->get<int>("PartialFacto");
354 if (partial_facto == 1 || partial_facto == 2) {
355 mumps_par.icntl[18] = 1; // centralized schur complement
356 }
357 }
358 if(parameterList->isParameter("SchurPart")) {
359 // copy schur-part to the internal view (in case the user-pointer go out of scope?)
360 auto schur_part_ptr = parameterList->get<const local_ordinal_type*>("SchurPart");
361 schur_size = 0;
362 for (global_size_type i=0; i<this->globalNumCols_; i++) {
363 if (schur_part_ptr[i] == 1) schur_size ++;
364 }
365 // allocate internal storage to store the schur complement (user may not want it and may not provide a valid pointer?)
366 Kokkos::resize(schur_part, schur_size);
367 Kokkos::resize(schur_out, schur_size*schur_size);
368 schur_size = 0;
369 for (global_size_type i=0; i<this->globalNumCols_; i++) {
370 if (schur_part_ptr[i] == 1) {
371 schur_part(schur_size) = i+1; // TODO: base-1 ??
372 schur_size ++;
373 }
374 }
375 mumps_par.size_schur = schur_size;
376 mumps_par.listvar_schur = schur_part.data();
377 mumps_par.schur = schur_out.data();
378 }
379 if(parameterList->isParameter("SchurOut")) {
380 // store schur-part to the internal view (if user wants the output, then the pointer should stay, so no need for internal view?)
381 schur_out_ptr = parameterList->get<scalar_type*>("SchurOut");
382 }
383 if(parameterList->isParameter("OnlyForwardSolve")) {
384 only_forward_solve = parameterList->get<bool>("OnlyForwardSolve");
385 }
386 if(parameterList->isParameter("OnlyBackwardSolve")) {
387 only_backward_solve = parameterList->get<bool>("OnlyBackwardSolve");
388 }
389
390 if(parameterList->isParameter("verbose")){
391 bool verbose = parameterList->get<bool>("verbose");
392 if (verbose) mumps_par.icntl[3] = 2;
393 }
394 if( parameterList->isParameter("IsContiguous") ){
395 is_contiguous_ = parameterList->get<bool>("IsContiguous");
396 }
397 }//end set parameters()
398
399
400 template <class Matrix, class Vector>
401 Teuchos::RCP<const Teuchos::ParameterList>
402 MUMPS<Matrix,Vector>::getValidParameters_impl() const
403 {
404 using Teuchos::ParameterList;
405
406 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
407
408 if( is_null(valid_params) ){
409 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
410
411 pl->set("ICNTL(1)", -1, "Output stream for error messages." );
412 pl->set("ICNTL(2)", -1, "Output stream for diagnostic." );
413 pl->set("ICNTL(3)", -1, "Output stream for global information." );
414 pl->set("ICNTL(4)", 1, "Level of printing." );
415 pl->set("ICNTL(6)", 0, "Permutes the matrix to a zero-free diagonal" );
416 pl->set("ICNTL(9)", 1, "Transpose solve, if not 1" );
417 pl->set("ICNTL(11)", 0, "Computes statistics for error analysis" );
418 pl->set("ICNTL(14)", 20, "Percentage increase in the estimated working space" );
419
420 // For partial factorization
421 scalar_type *dummy_scalar_ptr;
422 const local_ordinal_type *dummy_ordinal_ptr;
423 pl->set("PartialFacto", 0,
424 "Perform partial factorization to extract dense Schur complement (0: no, 1: form + factor Schur, 2: ony form");
425 pl->set("SchurPart", dummy_ordinal_ptr,
426 "Specify rows/columns belonging to Schur complement for partial factorization");
427 pl->set("SchurOut", dummy_scalar_ptr,
428 "Store output Schur complement from partial factorization");
429 pl->set("OnlyForwardSolve", false,
430 "Perform only the forward substitution");
431 pl->set("OnlyBackwardSolve", false,
432 "Perform only the backward substitution");
433
434 pl->set("verbose", false, "Whether verbose");
435 pl->set("IsContiguous", true, "Whether GIDs contiguous");
436 valid_params = pl;
437 }
438
439 return valid_params;
440 }//end getValidParmaeters_impl()
441
442
443 template <class Matrix, class Vector>
444 bool
445 MUMPS<Matrix,Vector>::loadA_impl(EPhase current_phase)
446 {
447 #ifdef HAVE_AMESOS2_TIMERS
448 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
449 #endif
450 if(MUMPS_MATRIX_LOAD == false || current_phase==NUMFACT)
451 {
452 // Only the root image needs storage allocated
453 if( !MUMPS_MATRIX_LOAD && this->root_ ) {
454 Kokkos::resize(host_nzvals_view_, this->globalNumNonZeros_);
455 Kokkos::resize(host_rows_view_, this->globalNumNonZeros_);
456 Kokkos::resize(host_col_ptr_view_, this->globalNumRows_ + 1);
457 }
458
459 #ifdef HAVE_AMESOS2_TIMERS
460 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
461 #endif
462
463 local_ordinal_type nnz_ret = 0;
464 Util::get_ccs_helper_kokkos_view<
465 MatrixAdapter<Matrix>,host_value_type_view,host_ordinal_type_view,host_ordinal_type_view>
466 ::do_get(this->matrixA_.ptr(), host_nzvals_view_, host_rows_view_, host_col_ptr_view_, nnz_ret,
467 (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
468 ARBITRARY,
469 this->rowIndexBase_);
470
471 Teuchos::broadcast(*(this->getComm()), 0, &nnz_ret); // We only care about root's value
472 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != Teuchos::as<local_ordinal_type>(this->globalNumNonZeros_),
473 std::runtime_error, "Did not get the expected number of non-zero vals");
474
475 if( this->root_ ){
476 ConvertToTriplet();
477 }
478 /* ch: In general, the matrix is loaded during the preordering phase.
479 However, if the matrix pattern has not changed during consecutive calls of numeric factorizations
480 we can reuse the previous symbolic factorization. In this case, the matrix is not loaded in the preordering phase,
481 because it is not called. Therefore, we need to load the matrix in the numeric factorization phase.
482 */
483 if (current_phase==PREORDERING){
484 MUMPS_MATRIX_LOAD_PREORDERING = true;
485 }
486 }
487
488 MUMPS_MATRIX_LOAD = true;
489 return (true);
490 }//end loadA_impl()
491
492 template <class Matrix, class Vector>
493 void
494 MUMPS<Matrix,Vector>::describe_impl(Teuchos::FancyOStream &out,
495 const Teuchos::EVerbosityLevel verbLevel) const
496 {
497 out << " MUMPS current parameters:" << std::endl;
498 out << " > ICNTL(1) = " << mumps_par.icntl[0] << std::endl;
499 out << " > ICNTL(2) = " << mumps_par.icntl[1] << std::endl;
500 out << " > ICNTL(3) = " << mumps_par.icntl[2] << std::endl;
501 out << " > ICNTL(4) = " << mumps_par.icntl[3] << std::endl;
502 out << " > ICNTL(6) = " << mumps_par.icntl[5] << std::endl;
503 out << " > ICNTL(7) = " << mumps_par.icntl[6] << std::endl;
504 out << " > ICNTL(9) = " << mumps_par.icntl[8] << std::endl;
505 out << " > ICNTL(11) = " << mumps_par.icntl[10] << std::endl;
506 out << " > ICNTL(14) = " << mumps_par.icntl[13] << std::endl;
507 out << " > IsContiguous = " << (is_contiguous_ ? "YES" : "NO") << std::endl;
508 out << std::endl;
509 }
510
511 template <class Matrix, class Vector>
512 int
513 MUMPS<Matrix,Vector>::ConvertToTriplet()
514 {
515 if ( !MUMPS_STRUCT ) {
516 MUMPS_STRUCT = true;
517 mumps_par.n = this->globalNumCols_;
518 mumps_par.nz = this->globalNumNonZeros_;
519 mumps_par.a = (mumps_type*)malloc(mumps_par.nz * sizeof(mumps_type));
520 mumps_par.irn = (MUMPS_INT*)malloc(mumps_par.nz *sizeof(MUMPS_INT));
521 mumps_par.jcn = (MUMPS_INT*)malloc(mumps_par.nz * sizeof(MUMPS_INT));
522 }
523 if((mumps_par.a == NULL) || (mumps_par.irn == NULL)
524 || (mumps_par.jcn == NULL)) {
525 return -1;
526 }
527 /* Going from full CSC to full Triplet */
528 /* Will have to add support for symmetric case*/
529 local_ordinal_type tri_count = 0;
530 local_ordinal_type i,j;
531 local_ordinal_type max_local_ordinal = 0;
532
533 for(i = 0; i < (local_ordinal_type)this->globalNumCols_; i++) {
534 for( j = host_col_ptr_view_(i); j < host_col_ptr_view_(i+1)-1; j++) {
535 mumps_par.jcn[tri_count] = (MUMPS_INT)i+1; //Fortran index
536 mumps_par.irn[tri_count] = (MUMPS_INT)host_rows_view_(j)+1; //Fortran index
537 mumps_par.a[tri_count] = host_nzvals_view_(j);
538
539 tri_count++;
540 }
541
542 j = host_col_ptr_view_(i+1)-1;
543 mumps_par.jcn[tri_count] = (MUMPS_INT)i+1; //Fortran index
544 mumps_par.irn[tri_count] = (MUMPS_INT)host_rows_view_(j)+1; //Fortran index
545 mumps_par.a[tri_count] = host_nzvals_view_(j);
546
547 tri_count++;
548
549 if(host_rows_view_(j) > max_local_ordinal) {
550 max_local_ordinal = host_rows_view_(j);
551 }
552 }
553 TEUCHOS_TEST_FOR_EXCEPTION(std::numeric_limits<MUMPS_INT>::max() <= max_local_ordinal,
554 std::runtime_error,
555 "Matrix index larger than MUMPS_INT");
556
557 return 0;
558 }//end Convert to Trip()
559
560 template<class Matrix, class Vector>
561 void
562 MUMPS<Matrix,Vector>::MUMPS_ERROR()const
563 {
564 using Teuchos::Comm;
565 using Teuchos::RCP;
566 bool Wrong = ((mumps_par.info[0] != 0) || (mumps_par.infog[0] != 0)) && (this->rank_ < this->nprocs_);
567 if(Wrong){
568 if (this->rank_==0) {
569 std::cerr << "Amesos2_Mumps : ERROR" << std::endl;
570 if ( this->status_.getNumSolve() > 0) {
571 std::cerr << " Last Phase : SOLVE" << std::endl;
572 } else if( this->status_.numericFactorizationDone() ){
573 std::cerr << " Last Phase : NUMFACT" << std::endl;
574 } else if( this->status_.symbolicFactorizationDone() ){
575 std::cerr << " Last Phase : SYMBFACT" << std::endl;
576 } else if( this->status_.preOrderingDone() ){
577 std::cerr << " Last Phase : PREORDERING" << std::endl;
578 } else {
579 std::cerr << " Last Phase : CLEAN" << std::endl;
580 }
581 std::cerr << "Amesos2_Mumps : INFOG(1) = " << mumps_par.infog[0] << std::endl;
582 std::cerr << "Amesos2_Mumps : INFOG(2) = " << mumps_par.infog[1] << std::endl;
583 }
584 if (mumps_par.info[0] != 0 && Wrong) {
585 std::cerr << "Amesos2_Mumps : On process " << this->matrixA_->getComm()->getRank()
586 << ", INFO(1) = " << mumps_par.info[0] << std::endl;
587 std::cerr << "Amesos2_Mumps : On process " << this->matrixA_->getComm()->getRank()
588 << ", INFO(2) = " << mumps_par.info[1] << std::endl;
589 }
590
591
592 }
593 // Throw on all ranks
594 int WrongInt = Wrong;
595 RCP<const Comm<int> > matComm = this->matrixA_->getComm();
596 Teuchos::broadcast<int,int>(*matComm,0,1,&WrongInt);
597 TEUCHOS_TEST_FOR_EXCEPTION(WrongInt>0,
598 std::runtime_error,
599 "MUMPS error");
600
601 }//end MUMPS_ERROR()
602
603
604 template<class Matrix, class Vector>
605 const char* MUMPS<Matrix,Vector>::name = "MUMPS";
606
607} // end namespace Amesos2
608
609#endif // AMESOS2_MUMPS_DEF_HPP
Amesos2_MUMPS_decl.hpp
Amesos2 MUMPS declarations.
Amesos2_SolverCore_def.hpp
Amesos2::ROOTED
@ ROOTED
Definition Amesos2_TypeDecl.hpp:93
Amesos2::CONTIGUOUS_AND_ROOTED
@ CONTIGUOUS_AND_ROOTED
Definition Amesos2_TypeDecl.hpp:94
Amesos2::ARBITRARY
@ ARBITRARY
Definition Amesos2_TypeDecl.hpp:109
Amesos2::MUMPS
Amesos2 interface to the MUMPS package.
Definition Amesos2_MUMPS_decl.hpp:51
Amesos2::MatrixAdapter
A Matrix adapter interface for Amesos2.
Definition Amesos2_MatrixAdapter_decl.hpp:42
Amesos2::EPhase
EPhase
Used to indicate a phase in the direct solution.
Definition Amesos2_TypeDecl.hpp:31
Amesos2::FunctionMap
Passes functions to TPL functions based on type.
Definition Amesos2_FunctionMap.hpp:43
Amesos2::MultiVecAdapter
A templated MultiVector class adapter for Amesos2.
Definition Amesos2_MultiVecAdapter_decl.hpp:142
Amesos2::Util::get_ccs_helper_kokkos_view
A generic helper class for getting a CCS representation of a Matrix.
Definition Amesos2_Util.hpp:589
Generated by doxygen 1.9.8