Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Amesos_Mumps.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: Direct Sparse Solver Package
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25//
26// ***********************************************************************
27// @HEADER
28
29#include "Amesos_Mumps.h"
30#include "Epetra_Map.h"
31#include "Epetra_Import.h"
32#include "Epetra_RowMatrix.h"
33#include "Epetra_CrsMatrix.h"
34#include "Epetra_Vector.h"
35#include "Epetra_Util.h"
36#include "Epetra_Time.h"
37#include "Epetra_IntSerialDenseVector.h"
38#include "Epetra_SerialDenseVector.h"
39#include "Epetra_IntVector.h"
40#include "Epetra_Vector.h"
41#include "Epetra_SerialDenseMatrix.h"
42#include "Epetra_Util.h"
43
44#define ICNTL(I) icntl[(I)-1]
45#define CNTL(I) cntl[(I)-1]
46#define INFOG(I) infog[(I)-1]
47#define INFO(I) info[(I)-1]
48#define RINFOG(I) rinfog[(I)-1]
49
50//=============================================================================
51Amesos_Mumps::Amesos_Mumps(const Epetra_LinearProblem &prob ) :
52 IsComputeSchurComplementOK_(false),
53 NoDestroy_(false),
54 MaxProcs_(-1),
55 UseTranspose_(false),
56 MtxConvTime_(-1),
57 MtxRedistTime_(-1),
58 VecRedistTime_(-1),
59 SymFactTime_(-1),
60 NumFactTime_(-1),
61 SolveTime_(-1),
62 RowSca_(0),
63 ColSca_(0),
64 PermIn_(0),
65 NumSchurComplementRows_(-1),
66#ifdef HAVE_MPI
67 MUMPSComm_(0),
68#endif
69 Problem_(&prob)
70{
71 // -777 is for me. It means: I never called MUMPS, so
72 // SymbolicFactorization should not call Destroy() and ask MUMPS to
73 // free its space.
74 MDS.job = -777;
75
76 // load up my default parameters
77 ICNTL[1] = -1; // Turn off error messages 6=on, -1 =off
78 ICNTL[2] = -1; // Turn off diagnostic printing 6=on, -1=off
79 ICNTL[3] = -1; // Turn off global information messages 6=on, -1=off
80 ICNTL[4] = -1; // 3 = most msgs; -1= none
81
82#if defined(MUMPS_4_9) || defined(MUMPS_5_0)
83
84 ICNTL[5] = 0; // Matrix is given in assembled (i.e. triplet) from
85 ICNTL[6] = 7; // Choose column permutation automatically
86 ICNTL[7] = 7; // Choose ordering method automatically
87 ICNTL[8] = 77; // Choose scaling automatically
88 ICNTL[9] = 1; // Compute A x = b
89 ICNTL[10] = 0; // Maximum steps of iterative refinement
90 ICNTL[11] = 0; // Do not collect statistics
91 ICNTL[12] = 0; // Automatic choice of ordering strategy
92 ICNTL[13] = 0; // Use ScaLAPACK for root node
93 ICNTL[14] = 20; // Increase memory allocation 20% at a time
94
95 // 15, 16 and 17 are not used
96 // 18 is set after we know NumProc
97 // 19 is set after we know Schur status
98
99 ICNTL[20] = 0; // RHS is given in dense form
100 ICNTL[21] = 0; // Solution is assembled at end, not left distributed
101 ICNTL[22] = 0; // Do all computations in-core
102 ICNTL[23] = 0; // We don't supply maximum MB of working memory
103 ICNTL[24] = 0; // Do not perform null pivot detection
104 ICNTL[25] = 0; // No null space basis computation, return 1 possible solution
105 ICNTL[26] = 0; // Do not condense/reduce Schur RHS
106 ICNTL[27] = -8; // Blocking factor for multiple RHSs
107 ICNTL[28] = 0; // Automatic choice of sequential/parallel analysis phase
108 ICNTL[29] = 0; // Parallel analysis uses PT-SCOTCH or ParMetis? (none)
109
110 // 30 - 40 are not used
111
112#else
113 ICNTL[5] = 0; // Matrix is given in elemental (i.e. triplet) from
114 ICNTL[6] = 7; // Choose column permutation automatically
115 ICNTL[7] = 7; // Choose ordering method automatically
116 ICNTL[8] = 7; // Choose scaling automatically
117 ICNTL[8] = 7; // Choose scaling automatically
118 ICNTL[9] = 1; // Compute A x = b
119 ICNTL[10] = 0; // Maximum steps of iterative refinement
120 ICNTL[11] = 0; // Do not collect statistics
121 ICNTL[12] = 0; // Use Node level parallelism
122 ICNTL[13] = 0; // Use ScaLAPACK for root node
123 ICNTL[14] = 20; // Increase memory allocation 20% at a time
124 ICNTL[15] = 0; // Minimize memory use (not flop count)
125 ICNTL[16] = 0; // Do not perform null space detection
126 ICNTL[17] = 0; // Unused (null space dimension)
127#endif
128
129 Teuchos::ParameterList ParamList;
130 SetParameters(ParamList);
131}
132
133//=============================================================================
135{
136 if (!NoDestroy_)
137 {
138 // destroy instance of the package
139 MDS.job = -2;
140
141 if (Comm().MyPID() < MaxProcs_) dmumps_c(&(MDS));
142
143#if 0
144 if (IsComputeSchurComplementOK_ && (Comm().MyPID() == 0)
145 && MDS.schur) {
146 delete [] MDS.schur;
147 MDS.schur = 0;
148 }
149#endif
150
151#ifdef HAVE_MPI
152 if (MUMPSComm_)
153 {
154 MPI_Comm_free( &MUMPSComm_ );
155 MUMPSComm_ = 0;
156 }
157#endif
158
159 if( (verbose_ && PrintTiming_) || verbose_ == 2) PrintTiming();
160 if( (verbose_ && PrintStatus_) || verbose_ == 2) PrintStatus();
161 }
162}
163
164//=============================================================================
169
170//=============================================================================
171int Amesos_Mumps::ConvertToTriplet(const bool OnlyValues)
172{
173
174 Epetra_RowMatrix* ptr;
175 if (Comm().NumProc() == 1)
176 ptr = &Matrix();
177 else {
178 ptr = &RedistrMatrix(true);
179 }
180
181 ResetTimer();
182
183#ifdef EXTRA_DEBUG_INFO
184 Epetra_CrsMatrix* Eptr = dynamic_cast<Epetra_CrsMatrix*>( ptr );
185 if ( ptr->NumGlobalNonzeros() < 300 ) SetICNTL(4,3 ); // Enable more debug info for small matrices
186 if ( ptr->NumGlobalNonzeros() < 42 && Eptr ) {
187 std::cout << " Matrix = " << std::endl ;
188 Eptr->Print( std::cout ) ;
189 } else {
190 assert( Eptr );
191 }
192#endif
193
194 Row.resize(ptr->NumMyNonzeros());
195 Col.resize(ptr->NumMyNonzeros());
196 Val.resize(ptr->NumMyNonzeros());
197
198 int MaxNumEntries = ptr->MaxNumEntries();
199 std::vector<int> Indices;
200 std::vector<double> Values;
201 Indices.resize(MaxNumEntries);
202 Values.resize(MaxNumEntries);
203
204 int count = 0;
205
206 for (int i = 0; i < ptr->NumMyRows() ; ++i) {
207
208 int GlobalRow = ptr->RowMatrixRowMap().GID(i);
209
210 int NumEntries = 0;
211 int ierr;
212 ierr = ptr->ExtractMyRowCopy(i, MaxNumEntries,
213 NumEntries, &Values[0],
214 &Indices[0]);
215 AMESOS_CHK_ERR(ierr);
216
217 for (int j = 0 ; j < NumEntries ; ++j) {
218 if (OnlyValues == false) {
219 Row[count] = GlobalRow + 1;
220 Col[count] = ptr->RowMatrixColMap().GID(Indices[j]) + 1;
221 }
222
223 // MS // Added on 15-Mar-05.
224 if (AddToDiag_ && Indices[j] == i)
225 Values[j] += AddToDiag_;
226
227 Val[count] = Values[j];
228 count++;
229 }
230 }
231
232 MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
233
234 assert (count <= ptr->NumMyNonzeros());
235
236 return(0);
237}
238
239//=============================================================================
240// NOTE: suppose first position is 1 (as in FORTRAN)
241void Amesos_Mumps::SetICNTL(int pos, int value)
242{
243 ICNTL[pos] = value;
244}
245
246//=============================================================================
247// NOTE: suppose first position is 1 (as in FORTRAN)
248void Amesos_Mumps::SetCNTL(int pos, double value)
249{
250 CNTL[pos] = value;
251}
252
253//=============================================================================
255{
256 std::map<int,int>::iterator i_iter;
257 for (i_iter = ICNTL.begin() ; i_iter != ICNTL.end() ; ++i_iter)
258 {
259 int pos = i_iter->first;
260 int val = i_iter->second;
261 if (pos < 1 || pos > 40) continue;
262 MDS.ICNTL(pos) = val;
263 }
264
265 std::map<int,double>::iterator d_iter;
266 for (d_iter = CNTL.begin() ; d_iter != CNTL.end() ; ++d_iter)
267 {
268 int pos = d_iter->first;
269 double val = d_iter->second;
270 if (pos < 1 || pos > 5) continue;
271 MDS.CNTL(pos) = val;
272 }
273
274 // fix some options
275
276 if (Comm().NumProc() != 1) MDS.ICNTL(18)= 3;
277 else MDS.ICNTL(18)= 0;
278
279 if (UseTranspose()) MDS.ICNTL(9) = 0;
280 else MDS.ICNTL(9) = 1;
281
282 MDS.ICNTL(5) = 0;
283
284#if 0
285 if (IsComputeSchurComplementOK_) MDS.ICNTL(19) = 1;
286 else MDS.ICNTL(19) = 0;
287
288 // something to do if the Schur complement is required.
289 if (IsComputeSchurComplementOK_ && Comm().MyPID() == 0)
290 {
291 MDS.size_schur = NumSchurComplementRows_;
292 MDS.listvar_schur = SchurComplementRows_;
294 }
295#endif
296}
297
298//=============================================================================
299int Amesos_Mumps::SetParameters( Teuchos::ParameterList & ParameterList)
300{
301 // ========================================= //
302 // retrive MUMPS' parameters from list. //
303 // default values defined in the constructor //
304 // ========================================= //
305
306 // retrive general parameters
307
308 SetStatusParameters(ParameterList);
309
310 SetControlParameters(ParameterList);
311
312 if (ParameterList.isParameter("NoDestroy"))
313 NoDestroy_ = ParameterList.get<bool>("NoDestroy");
314
315 // can be use using mumps sublist, via "ICNTL(2)"
316 // if (debug_)
317 // MDS.ICNTL(2) = 6; // Turn on Mumps verbose output
318
319 // retrive MUMPS' specific parameters
320
321 if (ParameterList.isSublist("mumps"))
322 {
323 Teuchos::ParameterList MumpsParams = ParameterList.sublist("mumps") ;
324
325 // ICNTL
326 for (int i = 1 ; i <= 40 ; ++i)
327 {
328 char what[80];
329 sprintf(what, "ICNTL(%d)", i);
330 if (MumpsParams.isParameter(what))
331 SetICNTL(i, MumpsParams.get<int>(what));
332 }
333
334 // CNTL
335 for (int i = 1 ; i <= 5 ; ++i)
336 {
337 char what[80];
338 sprintf(what, "CNTL(%d)", i);
339 if (MumpsParams.isParameter(what))
340 SetCNTL(i, MumpsParams.get<double>(what));
341 }
342
343 // ordering
344 if (MumpsParams.isParameter("PermIn"))
345 {
346 int* PermIn = 0;
347 PermIn = MumpsParams.get<int*>("PermIn");
348 if (PermIn) SetOrdering(PermIn);
349 }
350 // Col scaling
351 if (MumpsParams.isParameter("ColScaling"))
352 {
353 double * ColSca = 0;
354 ColSca = MumpsParams.get<double *>("ColScaling");
355 if (ColSca) SetColScaling(ColSca);
356 }
357 // Row scaling
358 if (MumpsParams.isParameter("RowScaling"))
359 {
360 double * RowSca = 0;
361 RowSca = MumpsParams.get<double *>("RowScaling");
362 if (RowSca) SetRowScaling(RowSca);
363 }
364 // that's all folks
365 }
366
367 return(0);
368}
369
370//=============================================================================
371
373{
374#ifndef HAVE_AMESOS_MPI_C2F
375 MaxProcs_ = -3;
376#endif
377
378 // check parameters and fix values of MaxProcs_
379
380 int NumGlobalNonzeros, NumRows;
381
382 NumGlobalNonzeros = Matrix().NumGlobalNonzeros();
383 NumRows = Matrix().NumGlobalRows();
384
385 // optimal value for MaxProcs == -1
386
387 int OptNumProcs1 = 1 + EPETRA_MAX(NumRows/10000, NumGlobalNonzeros/100000);
388 OptNumProcs1 = EPETRA_MIN(Comm().NumProc(),OptNumProcs1);
389
390 // optimal value for MaxProcs == -2
391
392 int OptNumProcs2 = (int)sqrt(1.0 * Comm().NumProc());
393 if (OptNumProcs2 < 1) OptNumProcs2 = 1;
394
395 // fix the value of MaxProcs
396
397 switch (MaxProcs_) {
398 case -1:
399 MaxProcs_ = OptNumProcs1;
400 break;
401 case -2:
402 MaxProcs_ = OptNumProcs2;
403 break;
404 case -3:
405 MaxProcs_ = Comm().NumProc();
406 break;
407 }
408
409 // few checks
410 if (MaxProcs_ > Comm().NumProc()) MaxProcs_ = Comm().NumProc();
411// if ( MaxProcs_ > 1 ) MaxProcs_ = Comm().NumProc(); // Bug - bogus kludge here - didn't work anyway
412}
413
414//=============================================================================
416{
417
418 // erase data if present.
419 if (IsSymbolicFactorizationOK_ && MDS.job != -777)
420 Destroy();
421
424
425 CreateTimer(Comm());
426
429
430#if defined(HAVE_MPI) && defined(HAVE_AMESOS_MPI_C2F)
431 if (MaxProcs_ != Comm().NumProc())
432 {
433 if(MUMPSComm_)
434 MPI_Comm_free(&MUMPSComm_);
435
436 std::vector<int> ProcsInGroup(MaxProcs_);
437 for (int i = 0 ; i < MaxProcs_ ; ++i)
438 ProcsInGroup[i] = i;
439
440 MPI_Group OrigGroup, MumpsGroup;
441 MPI_Comm_group(MPI_COMM_WORLD, &OrigGroup);
442 MPI_Group_incl(OrigGroup, MaxProcs_, &ProcsInGroup[0], &MumpsGroup);
443 MPI_Comm_create(MPI_COMM_WORLD, MumpsGroup, &MUMPSComm_);
444
445#ifdef MUMPS_4_9
446 MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f( MUMPSComm_);
447#else
448
449#ifndef HAVE_AMESOS_OLD_MUMPS
450 MDS.comm_fortran = (DMUMPS_INT) MPI_Comm_c2f( MUMPSComm_);
451#else
452 MDS.comm_fortran = (F_INT) MPI_Comm_c2f( MUMPSComm_);
453#endif
454
455#endif
456
457 }
458 else
459 {
460 const Epetra_MpiComm* MpiComm = dynamic_cast<const Epetra_MpiComm*>(&Comm());
461 assert (MpiComm != 0);
462#ifdef MUMPS_4_9
463 MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f(MpiComm->GetMpiComm());
464#else
465
466#ifndef HAVE_AMESOS_OLD_MUMPS
467 MDS.comm_fortran = (DMUMPS_INT) MPI_Comm_c2f(MpiComm->GetMpiComm());
468#else
469 MDS.comm_fortran = (F_INT) MPI_Comm_c2f(MpiComm->GetMpiComm());
470#endif
471
472#endif
473 }
474#else
475 // This next three lines of code were required to make Amesos_Mumps work
476 // with Ifpack_SubdomainFilter. They is usefull in all cases
477 // when using MUMPS on a subdomain.
478 const Epetra_MpiComm* MpiComm = dynamic_cast<const Epetra_MpiComm*>(&Comm());
479 assert (MpiComm != 0);
480 MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f(MpiComm->GetMpiComm());
481 // only thing I can do, use MPI_COMM_WORLD. This will work in serial as well
482 // Previously, the next line was uncommented, but we don't want MUMPS to work
483 // on the global MPI comm, but on the comm associated with the matrix
484 // MDS.comm_fortran = -987654;
485#endif
486
487 MDS.job = -1 ; // Initialization
488 MDS.par = 1 ; // Host IS involved in computations
489// MDS.sym = MatrixProperty_;
490 MDS.sym = 0; // MatrixProperty_ is unititalized. Furthermore MUMPS
491 // expects only half of the matrix to be provided for
492 // symmetric matrices. Hence setting MDS.sym to be non-zero
493 // indicating that the matrix is symmetric will only work
494 // if we change ConvertToTriplet to pass only half of the
495 // matrix. Bug #2331 and Bug #2332 - low priority
496
497
498 RedistrMatrix(true);
499
500 if (Comm().MyPID() < MaxProcs_)
501 {
502 dmumps_c(&(MDS)); // Initialize MUMPS
503 static_cast<void>( CheckError( ) );
504 }
505
506 MDS.n = Matrix().NumGlobalRows();
507
508 // fix pointers for nonzero pattern of A. Numerical values
509 // will be entered in PerformNumericalFactorization()
510 if (Comm().NumProc() != 1)
511 {
512 MDS.nz_loc = RedistrMatrix().NumMyNonzeros();
513
514 if (Comm().MyPID() < MaxProcs_)
515 {
516 MDS.irn_loc = &Row[0];
517 MDS.jcn_loc = &Col[0];
518 }
519 }
520 else
521 {
522 if (Comm().MyPID() == 0)
523 {
524 MDS.nz = Matrix().NumMyNonzeros();
525 MDS.irn = &Row[0];
526 MDS.jcn = &Col[0];
527 }
528 }
529
530 // scaling if provided by the user
531 if (RowSca_ != 0)
532 {
533 MDS.rowsca = RowSca_;
534 MDS.colsca = ColSca_;
535 }
536
537 // given ordering if provided by the user
538 if (PermIn_ != 0)
539 {
540 MDS.perm_in = PermIn_;
541 }
542
543 MDS.job = 1; // Request symbolic factorization
544
546
547 // Perform symbolic factorization
548
549 ResetTimer();
550
551 if (Comm().MyPID() < MaxProcs_)
552 dmumps_c(&(MDS));
553
554 SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
555
556 int IntWrong = CheckError()?1:0 ;
557 int AnyWrong;
558 Comm().SumAll( &IntWrong, &AnyWrong, 1 ) ;
559 bool Wrong = AnyWrong > 0 ;
560
561
562 if ( Wrong ) {
564 }
565
568
569 return 0;
570}
571
572//=============================================================================
573
575{
577
578 if (IsSymbolicFactorizationOK_ == false)
580
581 RedistrMatrix(true);
583
584 if (Comm().NumProc() != 1)
585 {
586 if (Comm().MyPID() < MaxProcs_)
587 MDS.a_loc = &Val[0];
588 }
589 else
590 MDS.a = &Val[0];
591
592 // Request numeric factorization
593 MDS.job = 2;
594 // Perform numeric factorization
595 ResetTimer();
596
597 if (Comm().MyPID() < MaxProcs_) {
598 dmumps_c(&(MDS));
599 }
600
601 NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
602
603 int IntWrong = CheckError()?1:0 ;
604 int AnyWrong;
605 Comm().SumAll( &IntWrong, &AnyWrong, 1 ) ;
606 bool Wrong = AnyWrong > 0 ;
607
608
609 if ( Wrong ) {
611 }
612
615 return(0);
616}
617
618//=============================================================================
619
621{
622 if (IsNumericFactorizationOK_ == false)
624
625 Epetra_MultiVector* vecX = Problem_->GetLHS();
626 Epetra_MultiVector* vecB = Problem_->GetRHS();
627 int NumVectors = vecX->NumVectors();
628
629 if ((vecX == 0) || (vecB == 0))
630 AMESOS_CHK_ERR(-1);
631
632 if (NumVectors != vecB->NumVectors())
633 AMESOS_CHK_ERR(-1);
634
635 if (Comm().NumProc() == 1)
636 {
637 // do not import any data
638 for (int j = 0 ; j < NumVectors; j++)
639 {
640 ResetTimer();
641
642 MDS.job = 3; // Request solve
643
644 for (int i = 0 ; i < Matrix().NumMyRows() ; ++i)
645 (*vecX)[j][i] = (*vecB)[j][i];
646 MDS.rhs = (*vecX)[j];
647
648 dmumps_c(&(MDS)) ; // Perform solve
649 static_cast<void>( CheckError( ) ); // Can hang
650 SolveTime_ = AddTime("Total solve time", SolveTime_);
651 }
652 }
653 else
654 {
655 Epetra_MultiVector SerialVector(SerialMap(),NumVectors);
656
657 ResetTimer();
658 AMESOS_CHK_ERR(SerialVector.Import(*vecB,SerialImporter(),Insert));
659 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
660
661 for (int j = 0 ; j < NumVectors; j++)
662 {
663 if (Comm().MyPID() == 0)
664 MDS.rhs = SerialVector[j];
665
666 // solve the linear system and take time
667 MDS.job = 3;
668 ResetTimer();
669 if (Comm().MyPID() < MaxProcs_)
670 dmumps_c(&(MDS)) ; // Perform solve
671 static_cast<void>( CheckError( ) ); // Can hang
672
673 SolveTime_ = AddTime("Total solve time", SolveTime_);
674 }
675
676 // ship solution back and take timing
677 ResetTimer();
678 AMESOS_CHK_ERR(vecX->Export(SerialVector,SerialImporter(),Insert));
679 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
680 }
681
683 ComputeTrueResidual(Matrix(), *vecX, *vecB, UseTranspose(), "Amesos_Mumps");
684
686 ComputeVectorNorms(*vecX, *vecB, "Amesos_Mumps");
687
688 NumSolve_++;
689
690 return(0) ;
691}
692
693#if 0
694//=============================================================================
695Epetra_CrsMatrix * Amesos_Mumps::GetCrsSchurComplement()
696{
697
699
700 if( Comm().MyPID() != 0 ) NumSchurComplementRows_ = 0;
701
702 Epetra_Map SCMap(-1,NumSchurComplementRows_, 0, Comm());
703
704 CrsSchurComplement_ = new Epetra_CrsMatrix(Copy,SCMap,NumSchurComplementRows_);
705
706 if( Comm().MyPID() == 0 )
707 for( int i=0 ; i<NumSchurComplementRows_ ; ++i ) {
708 for( int j=0 ; j<NumSchurComplementRows_ ; ++j ) {
709 int pos = i+ j *NumSchurComplementRows_;
710 CrsSchurComplement_->InsertGlobalValues(i,1,&(MDS.schur[pos]),&j);
711 }
712 }
713
714 CrsSchurComplement_->FillComplete();
715
716 return CrsSchurComplement_;
717 }
718
719 return 0;
720
721}
722
723//=============================================================================
724Epetra_SerialDenseMatrix * Amesos_Mumps::GetDenseSchurComplement()
725{
726
728
729 if( Comm().MyPID() != 0 ) return 0;
730
731 DenseSchurComplement_ = new Epetra_SerialDenseMatrix(NumSchurComplementRows_,
733
734 for( int i=0 ; i<NumSchurComplementRows_ ; ++i ) {
735 for( int j=0 ; j<NumSchurComplementRows_ ; ++j ) {
736 int pos = i+ j *NumSchurComplementRows_;
737 (*DenseSchurComplement_)(i,j) = MDS.schur[pos];
738 }
739 }
740
742
743 }
744
745 return(0);
746}
747
748//=============================================================================
749int Amesos_Mumps::ComputeSchurComplement(bool flag, int NumSchurComplementRows,
750 int * SchurComplementRows)
751{
752 NumSchurComplementRows_ = NumSchurComplementRows;
753
754 SchurComplementRows_ = SchurComplementRows;
755
756 // modify because MUMPS is fortran-driven
757 if( Comm().MyPID() == 0 )
758 for( int i=0 ; i<NumSchurComplementRows ; ++i ) SchurComplementRows_[i]++;
759
761
762 return 0;
763}
764#endif
765
766//=============================================================================
768{
769
770 if (Comm().MyPID() != 0 ) return;
771
772 // The following lines are commented out to deal with bug #1887 - kss
773#ifndef IRIX64
774 PrintLine();
775 std::cout << "Amesos_Mumps : Matrix has " << Matrix().NumGlobalRows() << " rows"
776 << " and " << Matrix().NumGlobalNonzeros() << " nonzeros" << std::endl;
777 std::cout << "Amesos_Mumps : Nonzero elements per row = "
778 << 1.0*Matrix().NumGlobalNonzeros()/Matrix().NumGlobalRows() << std::endl;
779 std::cout << "Amesos_Mumps : Percentage of nonzero elements = "
780 << 100.0*Matrix().NumGlobalNonzeros()/(pow(Matrix().NumGlobalRows(),2.0)) << std::endl;
781 std::cout << "Amesos_Mumps : Use transpose = " << UseTranspose_ << std::endl;
782// MatrixProperty_ is unused - see bug #2331 and bug #2332 in this file and bugzilla
783 if (MatrixProperty_ == 0) std::cout << "Amesos_Mumps : Matrix is general unsymmetric" << std::endl;
784 if (MatrixProperty_ == 2) std::cout << "Amesos_Mumps : Matrix is general symmetric" << std::endl;
785 if (MatrixProperty_ == 1) std::cout << "Amesos_Mumps : Matrix is SPD" << std::endl;
786 std::cout << "Amesos_Mumps : Available process(es) = " << Comm().NumProc() << std::endl;
787 std::cout << "Amesos_Mumps : Using " << MaxProcs_ << " process(es)" << std::endl;
788
789 std::cout << "Amesos_Mumps : Estimated FLOPS for elimination = "
790 << MDS.RINFOG(1) << std::endl;
791 std::cout << "Amesos_Mumps : Total FLOPS for assembly = "
792 << MDS.RINFOG(2) << std::endl;
793 std::cout << "Amesos_Mumps : Total FLOPS for elimination = "
794 << MDS.RINFOG(3) << std::endl;
795
796 std::cout << "Amesos_Mumps : Total real space to store the LU factors = "
797 << MDS.INFOG(9) << std::endl;
798 std::cout << "Amesos_Mumps : Total integer space to store the LU factors = "
799 << MDS.INFOG(10) << std::endl;
800 std::cout << "Amesos_Mumps : Total number of iterative steps refinement = "
801 << MDS.INFOG(15) << std::endl;
802 std::cout << "Amesos_Mumps : Estimated size of MUMPS internal data\n"
803 << "Amesos_Mumps : for running factorization = "
804 << MDS.INFOG(16) << " Mbytes" << std::endl;
805 std::cout << "Amesos_Mumps : for running factorization = "
806 << MDS.INFOG(17) << " Mbytes" << std::endl;
807 std::cout << "Amesos_Mumps : Allocated during factorization = "
808 << MDS.INFOG(19) << " Mbytes" << std::endl;
809 PrintLine();
810#endif
811}
812
813//=============================================================================
815{
816 bool Wrong = ((MDS.INFOG(1) != 0) || (MDS.INFO(1) != 0))
817 && (Comm().MyPID() < MaxProcs_);
818
819 // an error occurred in MUMPS. Print out information and quit.
820
821 if (Comm().MyPID() == 0 && Wrong)
822 {
823 std::cerr << "Amesos_Mumps : ERROR" << std::endl;
824 std::cerr << "Amesos_Mumps : INFOG(1) = " << MDS.INFOG(1) << std::endl;
825 std::cerr << "Amesos_Mumps : INFOG(2) = " << MDS.INFOG(2) << std::endl;
826 }
827
828 if (MDS.INFO(1) != 0 && Wrong)
829 {
830 std::cerr << "Amesos_Mumps : On process " << Comm().MyPID()
831 << ", INFO(1) = " << MDS.INFO(1) << std::endl;
832 std::cerr << "Amesos_Mumps : On process " << Comm().MyPID()
833 << ", INFO(2) = " << MDS.INFO(2) << std::endl;
834 }
835
836 if (Wrong)
837 return 1 ;
838 else
839 return 0 ;
840}
841
842// ======================================================================
844{
845 if (!Problem_->GetMatrix() || Comm().MyPID() != 0)
846 return;
847
848 double ConTime = GetTime(MtxConvTime_);
849 double MatTime = GetTime(MtxRedistTime_);
850 double VecTime = GetTime(VecRedistTime_);
851 double SymTime = GetTime(SymFactTime_);
852 double NumTime = GetTime(SymFactTime_);
853 double SolTime = GetTime(SolveTime_);
854
855 if (NumSymbolicFact_) SymTime /= NumSymbolicFact_;
856 if (NumNumericFact_) NumTime /= NumNumericFact_;
857 if (NumSolve_) SolTime /= NumSolve_;
858
859 std::string p = "Amesos_Mumps : ";
860 PrintLine();
861
862 std::cout << p << "Time to convert matrix to MUMPS format = "
863 << ConTime << " (s)" << std::endl;
864 std::cout << p << "Time to redistribute matrix = "
865 << MatTime << " (s)" << std::endl;
866 std::cout << p << "Time to redistribute vectors = "
867 << VecTime << " (s)" << std::endl;
868 std::cout << p << "Number of symbolic factorizations = "
869 << NumSymbolicFact_ << std::endl;
870 std::cout << p << "Time for sym fact = "
871 << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
872 std::cout << p << "Number of numeric factorizations = "
873 << NumNumericFact_ << std::endl;
874 std::cout << p << "Time for num fact = "
875 << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
876 std::cout << p << "Number of solve phases = "
877 << NumSolve_ << std::endl;
878 std::cout << p << "Time for solve = "
879 << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
880
881 PrintLine();
882}
883
884// ===========================================================================
885Epetra_RowMatrix& Amesos_Mumps::Matrix()
886{
887 Epetra_RowMatrix* Matrix = Problem_->GetMatrix();
888 assert (Matrix != 0);
889 return(*Matrix);
890}
891
892// ===========================================================================
893const Epetra_RowMatrix& Amesos_Mumps::Matrix() const
894{
895 Epetra_RowMatrix* Matrix = Problem_->GetMatrix();
896 assert (Matrix != 0);
897 return(*Matrix);
898}
899
900// ===========================================================================
902{
903 assert (Comm().NumProc() != 1);
904
905 if (RedistrMap_ == Teuchos::null) {
906 int i = Matrix().NumGlobalRows() / MaxProcs_;
907 if (Comm().MyPID() == 0)
908 i += Matrix().NumGlobalRows() % MaxProcs_;
909 else if (Comm().MyPID() >= MaxProcs_)
910 i = 0;
911
912 RedistrMap_ = rcp(new Epetra_Map(Matrix().NumGlobalRows(),i,0,Comm()));
913 assert (RedistrMap_.get() != 0);
914 }
915 return(*RedistrMap_);
916}
917
918// ===========================================================================
920{
921 assert (Comm().NumProc() != 1);
922
923 if (RedistrImporter_ == null) {
924 RedistrImporter_ = rcp(new Epetra_Import(RedistrMap(),Matrix().RowMatrixRowMap()));
925 assert (RedistrImporter_.get() != 0);
926 }
927 return(*RedistrImporter_);
928}
929
930// ===========================================================================
931Epetra_RowMatrix& Amesos_Mumps::RedistrMatrix(const bool ImportMatrix)
932{
933 if (Comm().NumProc() == 1) return(Matrix());
934
935 if (ImportMatrix) RedistrMatrix_ = null;
936
937 if (RedistrMatrix_ == null)
938 {
939 RedistrMatrix_ = rcp(new Epetra_CrsMatrix(Copy,RedistrMap(),0));
940 if (ImportMatrix) {
941 int ierr = RedistrMatrix_->Import(Matrix(),RedistrImporter(),Insert);
942 assert (ierr == 0);
943 ierr = RedistrMatrix_->FillComplete();
944 assert (ierr == 0);
945 }
946 }
947
948 return(*RedistrMatrix_);
949}
950
951// ===========================================================================
953{
954 if (SerialMap_ == null)
955 {
956 int i = Matrix().NumGlobalRows();
957 if (Comm().MyPID()) i = 0;
958 SerialMap_ = rcp(new Epetra_Map(-1,i,0,Comm()));
959 assert (SerialMap_ != null);
960 }
961 return(*SerialMap_);
962}
963
964// ===========================================================================
966{
967 if (SerialImporter_ == null) {
968 SerialImporter_ = rcp(new Epetra_Import(SerialMap(),Matrix().OperatorDomainMap()));
969 assert (SerialImporter_ != null);
970 }
971 return(*SerialImporter_);
972}
973
974// ===========================================================================
976{
977 return ( MDS.rinfo);
978}
979
980// ===========================================================================
982{
983 return (MDS.info);
984}
985
986// ===========================================================================
988{
989 return (MDS.rinfog);
990}
991
992// ===========================================================================
994{
995 return (MDS.infog);
996}
997
const int NumericallySingularMatrixError
const int StructurallySingularMatrixError
#define AMESOS_CHK_ERR(a)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
double AddToDiag_
Add this value to the diagonal.
int MatrixProperty_
Set the matrix property.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int Solve()
Solves A X = B (or AT x = B)
RCP< Epetra_Map > SerialMap_
Map with all elements on process 0 (for solution and rhs).
double * GetRINFO()
Gets the pointer to the RINFO array (defined on all processes).
bool UseTranspose_
If true, solve the problem with AT.
Epetra_RowMatrix & Matrix()
Returns a reference to the linear system matrix.
void SetCNTL(int pos, double value)
Set CNTL[pos] to value. pos is expressed in FORTRAN style (starting from 1).
DMUMPS_STRUC_C MDS
void SetICNTL(int pos, int value)
Set ICNTL[pos] to value. pos is expressed in FORTRAN style (starting from 1).
Epetra_Map & RedistrMap()
Returns a reference to the map for redistributed matrix.
double * RowSca_
Row and column scaling.
int MtxConvTime_
Quick access pointers to the internal timers.
void PrintStatus() const
Prints information about the factorization and solution phases.
std::map< int, double > CNTL
Epetra_Import & RedistrImporter()
Returns a reference for the redistributed importer.
void SetICNTLandCNTL()
std::map< int, int > ICNTL
int NumericFactorization()
Performs NumericFactorization on the matrix A.
Epetra_RowMatrix & RedistrMatrix(const bool ImportMatrix=false)
Returns a reference to the redistributed matrix, imports it is ImportMatrix is true.
int * GetINFO()
Gets the pointer to the INFO array (defined on all processes).
int * SchurComplementRows_
Rows for the Schur complement (if required)
bool IsComputeSchurComplementOK_
true if the Schur complement has been computed (need to free memory)
RCP< Epetra_CrsMatrix > CrsSchurComplement_
Pointer to the Schur complement, as CrsMatrix.
int SetColScaling(double *ColSca)
Set column scaling.
int NumSchurComplementRows_
Number of rows in the Schur complement (if required)
std::vector< double > Val
values of nonzero elements
int SetOrdering(int *PermIn)
Sets ordering.
Amesos_Mumps(const Epetra_LinearProblem &LinearProblem)
Amesos_Mumps Constructor.
int MaxProcs_
Maximum number of processors in the MUMPS' communicator.
double * GetRINFOG()
Gets the pointer to the RINFOG array (defined on host only).
RCP< Epetra_Map > RedistrMap_
Redistributed matrix.
RCP< Epetra_CrsMatrix > RedistrMatrix_
Redistributed matrix (only if MaxProcs_ > 1).
double * ColSca_
RCP< Epetra_Import > RedistrImporter_
Redistributed importer (from Matrix().RowMatrixRowMap() to RedistrMatrix().RowMatrixRowMap()).
std::vector< int > Col
column indices of nonzero elements
int * GetINFOG()
Get the pointer to the INFOG array (defined on host only).
std::vector< int > Row
row indices of nonzero elements
void PrintTiming() const
Prints timing information.
bool UseTranspose() const
Returns the current UseTranspose setting.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int * PermIn_
PermIn for MUMPS.
Epetra_Import & SerialImporter()
Returns a reference to the importer for SerialMap().
int CheckError()
Checks for MUMPS error, prints them if any. See MUMPS' manual.
int ConvertToTriplet(const bool OnlyValues)
Converts to MUMPS format (COO format).
const Epetra_LinearProblem * Problem_
Pointer to the linear problem to be solved.
void CheckParameters()
Check input parameters.
RCP< Epetra_Import > SerialImporter_
Importer from Matrix.OperatorDomainMap() to SerialMap_.
RCP< Epetra_SerialDenseMatrix > DenseSchurComplement_
Pointer to the Schur complement,as DenseMatrix.
int SetRowScaling(double *RowSca)
Set row scaling.
Epetra_Map & SerialMap()
Returns a reference to the map with all elements on process 0.
~Amesos_Mumps(void)
Amesos_Mumps Destructor.
void Destroy()
Destroys all data associated with \sl this object.
bool PrintTiming_
If true, prints timing information in the destructor.
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
int verbose_
Toggles the output level.
int NumSymbolicFact_
Number of symbolic factorization phases.
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
int NumSolve_
Number of solves.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
bool PrintStatus_
If true, print additional information in the destructor.
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
int NumNumericFact_
Number of numeric factorization phases.
double GetTime(const std::string what) const
Gets the cumulative time using the string.
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
Definition Amesos_Time.h:74
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
Definition Amesos_Time.h:80
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
Definition Amesos_Time.h:64
void ComputeTrueResidual(const Epetra_RowMatrix &Matrix, const Epetra_MultiVector &X, const Epetra_MultiVector &B, const bool UseTranspose, const std::string prefix) const
Computes the true residual, B - Matrix * X, and prints the results.
void ComputeVectorNorms(const Epetra_MultiVector &X, const Epetra_MultiVector &B, const std::string prefix) const
Computes the norms of X and B and print the results.
void PrintLine() const
Prints line on std::cout.