Amesos Package Browser (Single Doxygen Collection)  Development
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 //=============================================================================
51 Amesos_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 //=============================================================================
166 {
167  Destroy();
168 }
169 
170 //=============================================================================
171 int 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)
241 void 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)
248 void 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 //=============================================================================
300 {
301  // ========================================= //
302  // retrive MUMPS' parameters from list. //
303  // default values defined in the constructor //
304  // ========================================= //
305 
306  // retrive general parameters
307 
309 
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 
427  CheckParameters();
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 
545  SetICNTLandCNTL();
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 
567  NumSymbolicFact_++;
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 
614  NumNumericFact_++;
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 //=============================================================================
695 Epetra_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 //=============================================================================
724 Epetra_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 
741  return DenseSchurComplement_;
742 
743  }
744 
745  return(0);
746 }
747 
748 //=============================================================================
749 int 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_->GetOperator() == 0 || 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 // ===========================================================================
885 Epetra_RowMatrix& Amesos_Mumps::Matrix()
886 {
887  Epetra_RowMatrix* Matrix = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
888  assert (Matrix != 0);
889  return(*Matrix);
890 }
891 
892 // ===========================================================================
893 const Epetra_RowMatrix& Amesos_Mumps::Matrix() const
894 {
895  Epetra_RowMatrix* Matrix = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
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 // ===========================================================================
931 Epetra_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 
RCP< Epetra_Import > SerialImporter_
Importer from Matrix.OperatorDomainMap() to SerialMap_.
Definition: Amesos_Mumps.h:415
int NumSymbolicFact_
Number of symbolic factorization phases.
Definition: Amesos_Status.h:67
~Amesos_Mumps(void)
Amesos_Mumps Destructor.
double GetTime(const std::string what) const
Gets the cumulative time using the string.
Definition: Amesos_Time.h:102
int Solve()
Solves A X = B (or AT x = B)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
std::map< int, int > ICNTL
Definition: Amesos_Mumps.h:419
RCP< Epetra_CrsMatrix > CrsSchurComplement_
Pointer to the Schur complement, as CrsMatrix.
Definition: Amesos_Mumps.h:394
RCP< Epetra_Map > SerialMap_
Map with all elements on process 0 (for solution and rhs).
Definition: Amesos_Mumps.h:413
int * GetINFO()
Gets the pointer to the INFO array (defined on all processes).
int NumSchurComplementRows_
Number of rows in the Schur complement (if required)
Definition: Amesos_Mumps.h:389
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
Definition: Amesos_Status.h:48
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
T & get(ParameterList &l, const std::string &name)
std::map< int, double > CNTL
Definition: Amesos_Mumps.h:420
int * PermIn_
PermIn for MUMPS.
Definition: Amesos_Mumps.h:386
double * GetRINFO()
Gets the pointer to the RINFO array (defined on all processes).
int MtxConvTime_
Quick access pointers to the internal timers.
Definition: Amesos_Mumps.h:379
bool isSublist(const std::string &name) const
void SetCNTL(int pos, double value)
Set CNTL[pos] to value. pos is expressed in FORTRAN style (starting from 1).
RCP< ParameterList > sublist(const RCP< ParameterList > &paramList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
int * SchurComplementRows_
Rows for the Schur complement (if required)
Definition: Amesos_Mumps.h:391
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
Definition: Amesos_Time.h:64
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
RCP< Epetra_Map > RedistrMap_
Redistributed matrix.
Definition: Amesos_Mumps.h:407
const Epetra_LinearProblem * Problem_
Pointer to the linear problem to be solved.
Definition: Amesos_Mumps.h:404
std::vector< int > Col
column indices of nonzero elements
Definition: Amesos_Mumps.h:368
int NumNumericFact_
Number of numeric factorization phases.
Definition: Amesos_Status.h:69
T * get() const
Epetra_Import & RedistrImporter()
Returns a reference for the redistributed importer.
int NumSolve_
Number of solves.
Definition: Amesos_Status.h:71
Epetra_Import & SerialImporter()
Returns a reference to the importer for SerialMap().
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
Definition: Amesos_Status.h:56
int MatrixProperty_
Set the matrix property.
int SetColScaling(double *ColSca)
Set column scaling.
Definition: Amesos_Mumps.h:258
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
#define AMESOS_CHK_ERR(a)
Epetra_RowMatrix & RedistrMatrix(const bool ImportMatrix=false)
Returns a reference to the redistributed matrix, imports it is ImportMatrix is true.
int ConvertToTriplet(const bool OnlyValues)
Converts to MUMPS format (COO format).
void PrintStatus() const
Prints information about the factorization and solution phases.
int MaxProcs_
Maximum number of processors in the MUMPS&#39; communicator.
Definition: Amesos_Mumps.h:373
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.
Definition: Amesos_Utils.h:29
bool PrintTiming_
If true, prints timing information in the destructor.
Definition: Amesos_Status.h:52
bool PrintStatus_
If true, print additional information in the destructor.
Definition: Amesos_Status.h:54
void SetICNTL(int pos, int value)
Set ICNTL[pos] to value. pos is expressed in FORTRAN style (starting from 1).
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
DMUMPS_STRUC_C MDS
Definition: Amesos_Mumps.h:417
RCP< Epetra_CrsMatrix > RedistrMatrix_
Redistributed matrix (only if MaxProcs_ > 1).
Definition: Amesos_Mumps.h:411
double * ColSca_
Definition: Amesos_Mumps.h:383
bool UseTranspose() const
Returns the current UseTranspose setting.
Definition: Amesos_Mumps.h:153
int SetOrdering(int *PermIn)
Sets ordering.
Definition: Amesos_Mumps.h:269
Amesos_Mumps(const Epetra_LinearProblem &LinearProblem)
Amesos_Mumps Constructor.
Epetra_RowMatrix & Matrix()
Returns a reference to the linear system matrix.
std::vector< int > Row
row indices of nonzero elements
Definition: Amesos_Mumps.h:366
int CheckError()
Checks for MUMPS error, prints them if any. See MUMPS&#39; manual.
Epetra_Map & SerialMap()
Returns a reference to the map with all elements on process 0.
void PrintTiming() const
Prints timing information.
const int NumericallySingularMatrixError
RCP< Epetra_SerialDenseMatrix > DenseSchurComplement_
Pointer to the Schur complement,as DenseMatrix.
Definition: Amesos_Mumps.h:396
Copy
void SetICNTLandCNTL()
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.
Definition: Amesos_Utils.h:52
const int StructurallySingularMatrixError
bool isParameter(const std::string &name) const
std::vector< double > Val
values of nonzero elements
Definition: Amesos_Mumps.h:370
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
Definition: Amesos_Time.h:74
bool IsComputeSchurComplementOK_
true if the Schur complement has been computed (need to free memory)
Definition: Amesos_Mumps.h:360
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
Definition: Amesos_Mumps.h:319
void CheckParameters()
Check input parameters.
int * GetINFOG()
Get the pointer to the INFOG array (defined on host only).
int SetRowScaling(double *RowSca)
Set row scaling.
Definition: Amesos_Mumps.h:246
int verbose_
Toggles the output level.
Definition: Amesos_Status.h:61
RCP< Epetra_Import > RedistrImporter_
Redistributed importer (from Matrix().RowMatrixRowMap() to RedistrMatrix().RowMatrixRowMap()).
Definition: Amesos_Mumps.h:409
double * RowSca_
Row and column scaling.
Definition: Amesos_Mumps.h:383
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
Definition: Amesos_Status.h:50
bool UseTranspose_
If true, solve the problem with AT.
Definition: Amesos_Mumps.h:376
void PrintLine() const
Prints line on std::cout.
Definition: Amesos_Utils.h:71
void Destroy()
Destroys all data associated with this object.
double * GetRINFOG()
Gets the pointer to the RINFOG array (defined on host only).
Epetra_Map & RedistrMap()
Returns a reference to the map for redistributed matrix.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
Definition: Amesos_Status.h:58
double AddToDiag_
Add this value to the diagonal.