MFEM  v4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
strumpack.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_STRUMPACK
15 #ifdef MFEM_USE_MPI
16 
17 #include "strumpack.hpp"
18 
19 using namespace std;
20 using namespace strumpack;
21 
22 namespace mfem
23 {
24 
25 STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(MPI_Comm comm,
26  int num_loc_rows, int first_loc_row,
27  int glob_nrows, int glob_ncols,
28  int *I, int *J, double *data)
29  : comm_(comm), A_(NULL)
30 {
31  // Set mfem::Operator member data
32  height = num_loc_rows;
33  width = num_loc_rows;
34 
35  // Allocate STRUMPACK's CSRMatrixMPI
36  int nprocs, rank;
37  MPI_Comm_rank(comm_, &rank);
38  MPI_Comm_size(comm_, &nprocs);
39  int * dist = new int[nprocs + 1];
40  dist[rank + 1] = first_loc_row + num_loc_rows;
41  dist[0] = 0;
42  MPI_Allgather(MPI_IN_PLACE, 0, MPI_INT, dist + 1, 1, MPI_INT, comm_);
43  A_ = new CSRMatrixMPI<double,int>(num_loc_rows, I, J, data, dist, comm_, false);
44  delete[] dist;
45 }
46 
48  : comm_(hypParMat.GetComm()),
49  A_(NULL)
50 {
51  // First cast the parameter to a hypre_ParCSRMatrix
52  hypre_ParCSRMatrix * parcsr_op =
53  (hypre_ParCSRMatrix *)const_cast<HypreParMatrix&>(hypParMat);
54 
55  MFEM_ASSERT(parcsr_op != NULL,"STRUMPACK: const_cast failed in SetOperator");
56 
57  // Create the CSRMatrixMPI A_ by borrowing the internal data from a
58  // hypre_CSRMatrix.
59  hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
60  hypre_CSRMatrixSetDataOwner(csr_op,0);
61 
62  height = csr_op->num_rows;
63  width = csr_op->num_rows;
64 
65  int nprocs, rank;
66  MPI_Comm_rank(comm_, &rank);
67  MPI_Comm_size(comm_, &nprocs);
68  int * dist = new int[nprocs + 1];
69  dist[rank + 1] = parcsr_op->first_row_index + csr_op->num_rows;
70  dist[0] = 0;
71  MPI_Allgather(MPI_IN_PLACE, 0, MPI_INT, dist + 1, 1, MPI_INT, comm_);
72  A_ = new CSRMatrixMPI<double,int>(csr_op->num_rows, csr_op->i, csr_op->j,
73  csr_op->data, dist, comm_, false);
74  delete[] dist;
75 
76  // Everything has been copied or abducted so delete the structure
77  hypre_CSRMatrixDestroy(csr_op);
78 }
79 
81 {
82  // Delete the struct
83  if ( A_ != NULL ) { delete A_; }
84 }
85 
86 STRUMPACKSolver::STRUMPACKSolver( int argc, char* argv[], MPI_Comm comm )
87  : comm_(comm),
88  APtr_(NULL),
89  solver_(NULL)
90 {
91  this->Init(argc, argv);
92 }
93 
95  : comm_(A.GetComm()),
96  APtr_(&A),
97  solver_(NULL)
98 {
99  height = A.Height();
100  width = A.Width();
101 
102  this->Init(0, NULL);
103 }
104 
106 {
107  if ( solver_ != NULL ) { delete solver_; }
108 }
109 
110 void STRUMPACKSolver::Init( int argc, char* argv[] )
111 {
112  MPI_Comm_size(comm_, &numProcs_);
113  MPI_Comm_rank(comm_, &myid_);
114 
115  factor_verbose_ = false;
116  solve_verbose_ = false;
117 
118  solver_ = new StrumpackSparseSolverMPIDist<double,int>(comm_, argc, argv,
119  false);
120 }
121 
123 {
124  solver_->options().set_from_command_line( );
125 }
126 
128 {
129  factor_verbose_ = print_stat;
130 }
131 
133 {
134  solve_verbose_ = print_stat;
135 }
136 
137 void STRUMPACKSolver::SetKrylovSolver( strumpack::KrylovSolver method )
138 {
139  solver_->options().set_Krylov_solver( method );
140 }
141 
142 void STRUMPACKSolver::SetReorderingStrategy( strumpack::ReorderingStrategy
143  method )
144 {
145  solver_->options().set_reordering_method( method );
146 }
147 
149 {
150 #if STRUMPACK_VERSION_MAJOR >= 3
151  solver_->options().set_matching( strumpack::MatchingJob::NONE );
152 #else
153  solver_->options().set_mc64job( strumpack::MC64Job::NONE );
154 #endif
155 }
156 
158 {
159 #if STRUMPACK_VERSION_MAJOR >= 3
160  solver_->options().set_matching
161  ( strumpack::MatchingJob::MAX_DIAGONAL_PRODUCT_SCALING );
162 #else
163  solver_->options().set_mc64job
164  ( strumpack::MC64Job::MAX_DIAGONAL_PRODUCT_SCALING );
165 #endif
166 }
167 
168 #if STRUMPACK_VERSION_MAJOR >= 3
170 {
171  solver_->options().set_matching
172  ( strumpack::MatchingJob::COMBBLAS );
173 }
174 #endif
175 
176 void STRUMPACKSolver::SetRelTol( double rtol )
177 {
178  solver_->options().set_rel_tol( rtol );
179 }
180 
181 void STRUMPACKSolver::SetAbsTol( double atol )
182 {
183  solver_->options().set_abs_tol( atol );
184 }
185 
186 
187 void STRUMPACKSolver::Mult( const Vector & x, Vector & y ) const
188 {
189  MFEM_ASSERT(APtr_ != NULL,
190  "STRUMPACK Error: The operator must be set before"
191  " the system can be solved.");
192  MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
193  << ", expected size = " << Width());
194  MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
195  << ", expected size = " << Height());
196 
197  double* yPtr = (double*)y;
198  double* xPtr = (double*)(const_cast<Vector&>(x));
199 
200  solver_->options().set_verbose( factor_verbose_ );
201  ReturnCode ret = solver_->factor();
202  switch (ret)
203  {
204  case ReturnCode::SUCCESS: break;
205  case ReturnCode::MATRIX_NOT_SET:
206  {
207  MFEM_ABORT("STRUMPACK: Matrix was not set!");
208  }
209  break;
210  case ReturnCode::REORDERING_ERROR:
211  {
212  MFEM_ABORT("STRUMPACK: Matrix reordering failed!");
213  }
214  break;
215  }
216  solver_->options().set_verbose( solve_verbose_ );
217  solver_->solve(xPtr, yPtr);
218 
219 }
220 
222 {
223  // Verify that we have a compatible operator
224  APtr_ = dynamic_cast<const STRUMPACKRowLocMatrix*>(&op);
225  if ( APtr_ == NULL )
226  {
227  mfem_error("STRUMPACKSolver::SetOperator : not STRUMPACKRowLocMatrix!");
228  }
229 
230  solver_->set_matrix( *(APtr_->getA()) );
231 
232  // Set mfem::Operator member data
233  height = op.Height();
234  width = op.Width();
235 
236 }
237 
238 } // mfem namespace
239 
240 #endif // MFEM_USE_MPI
241 #endif // MFEM_USE_STRUMPACK
const STRUMPACKRowLocMatrix * APtr_
Definition: strumpack.hpp:161
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
STRUMPACKSolver(int argc, char *argv[], MPI_Comm comm)
Definition: strumpack.cpp:86
int Size() const
Returns the size of the vector.
Definition: vector.hpp:150
strumpack::StrumpackSparseSolverMPIDist< double, int > * solver_
Definition: strumpack.hpp:162
strumpack::CSRMatrixMPI< double, int > * getA() const
Definition: strumpack.hpp:55
void SetPrintFactorStatistics(bool print_stat)
Definition: strumpack.cpp:127
STRUMPACKRowLocMatrix(MPI_Comm comm, int num_loc_rows, int first_loc_row, int glob_nrows, int glob_ncols, int *I, int *J, double *data)
Definition: strumpack.cpp:25
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:146
void SetPrintSolveStatistics(bool print_stat)
Definition: strumpack.cpp:132
void SetReorderingStrategy(strumpack::ReorderingStrategy method)
Definition: strumpack.cpp:142
void SetRelTol(double rtol)
Definition: strumpack.cpp:176
A class to initialize the size of a Tensor.
Definition: dtensor.hpp:55
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: strumpack.cpp:187
void SetKrylovSolver(strumpack::KrylovSolver method)
Definition: strumpack.cpp:137
void SetAbsTol(double atol)
Definition: strumpack.cpp:181
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: strumpack.cpp:221
Vector data type.
Definition: vector.hpp:48
Abstract operator.
Definition: operator.hpp:21
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:187
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25