Loading...
Searching...
No Matches
solver_mumps.cc
Go to the documentation of this file.
// direct solver MUMPS, mpi implementations
//
// Note : why no seq implementation based on mumps ?
//
// 1. mumps has a pure seq library implementation.
// Both seq and mpi libraries cannot be linked together because they define
// both the same dmumps() function.
// 2. but, when using mpi, the mpi implementation is unable to solve efficiently
// in // some independant sparse linear systems (eg. sparse mass matrix local
// on an element) because the right-hand side may be merged on the proc=0
// for an obscure raison.
//
// Thus rheolef uses mumps for mpi-only implementation of the direct solver
// and use another library (eg umfpack) for the seq direct solver.
//
// The present 'seq' implementation has been tested and is invalid :
// the code is maintained (but no more instanciated) for two reasons:
// 1. it can be used when mpi is not present, as a seq solver with the pure seq impl of mumps
// 2. mumps will propose in the future to provide a distributed RHS:
// https://listes.ens-lyon.fr/sympa/arc/mumps-users
// From: Alfredo Buttari <alfredo.buttari@enseeiht.fr>
// Date: Mon, 31 Jan 2011 19:07:03 +0100
// Subject: Re: [mumps-users] Scalability of MUMPS solution phase
// Jack,
// yes, adding support for distributed RHS is on our TODO list but I
// can't say at this moment when it will be available.
//
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_MUMPS
#include "solver_mumps.h"
#undef DEBUG_MUMPS_SCOTCH
#ifdef DEBUG_MUMPS_SCOTCH
#undef _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH // has failed on some linear systems: prefer seq-scotch when available
#endif // DEBUG_MUMPS
namespace rheolef {
using namespace std;
// =========================================================================
// 1. local functions
// =========================================================================
// count non-zero entries of the upper part of A
// 1) diagonal block case
template<class T, class M>
static
typename csr<T,M>::size_type
nnz_dia_upper (const csr<T,M>& a)
{
typedef typename csr<T,M>::size_type size_type;
size_type dia_nnz = 0;
typename csr<T,M>::const_iterator dia_ia = a.begin();
for (size_type i = 0, n = a.nrow(); i < n; i++) {
for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
size_type j = (*p).first;
if (j >= i) dia_nnz++;
}
}
return dia_nnz;
}
// 2) extra-diagonal block case, for distributed csr
template<class T, class M>
static
typename csr<T,M>::size_type
nnz_ext_upper (const csr<T,M>& a)
{
return 0;
}
#ifdef _RHEOLEF_HAVE_MPI
template<class T>
static
typename csr<T,distributed>::size_type
nnz_ext_upper (const csr<T,distributed>& a)
{
typedef typename csr<T,distributed>::size_type size_type;
size_type first_i = a.row_ownership().first_index();
size_type ext_nnz = 0;
typename csr<T,distributed>::const_iterator ext_ia = a.ext_begin();
for (size_type i = 0, n = a.nrow(); i < n; i++) {
for (typename csr<T,distributed>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
size_type dis_i = first_i + i;
size_type dis_j = a.jext2dis_j ((*p).first);
if (dis_j >= dis_i) ext_nnz++;
}
}
return ext_nnz;
}
#endif // _RHEOLEF_HAVE_MPI
template<class T, class M>
static
void
csr2mumps_struct_seq (const csr<T,M>& a, vector<int>& row, vector<int>& col, bool use_symmetry)
{
typedef typename csr<T,M>::size_type size_type;
typename csr<T,M>::const_iterator dia_ia = a.begin();
for (size_type i = 0, n = a.nrow(), q = 0; i < n; i++) {
for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
size_type j = (*p).first;
if (!use_symmetry || j >= i) {
row[q] = i + 1;
col[q] = j + 1;
q++;
}
}
}
}
template<class T, class M>
static
void
csr2mumps_values_seq (const csr<T,M>& a, vector<T>& val, bool use_symmetry)
{
typedef typename csr<T,M>::size_type size_type;
typename csr<T,M>::const_iterator dia_ia = a.begin();
for (size_type i = 0, n = a.nrow(), q = 0; i < n; i++) {
for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
size_type j = (*p).first;
if (!use_symmetry || j >= i) {
val[q] = (*p).second;
q++;
}
}
}
}
#ifdef _RHEOLEF_HAVE_MPI
template<class T>
static
void
csr2mumps_struct_dis (const csr<T,sequential>& a, vector<int>& row, vector<int>& col, bool use_symmetry, bool)
{
// dummy case: for the compiler to be happy
}
template<class T>
static
void
csr2mumps_struct_dis (const csr<T,distributed>& a, vector<int>& row, vector<int>& col, bool use_symmetry, bool drop_ext_nnz)
{
typedef typename csr<T,distributed>::size_type size_type;
typename csr<T,distributed>::const_iterator dia_ia = a.begin();
typename csr<T,distributed>::const_iterator ext_ia = a.ext_begin();
size_type dis_nr = a.dis_nrow();
size_type dis_nc = a.dis_ncol();
size_type nc = a.ncol();
size_type first_i = a.row_ownership().first_index();
for (size_type i = 0, n = a.nrow(), q = 0; i < n; i++) {
for (typename csr<T,distributed>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
size_type j = (*p).first;
size_type dis_i = first_i + i;
size_type dis_j = first_i + j;
assert_macro (j < nc, "invalid matrix");
assert_macro (dis_i < dis_nr && dis_j < dis_nc, "invalid matrix");
if (!use_symmetry || dis_j >= dis_i) {
row[q] = dis_i + 1;
col[q] = dis_j + 1;
q++;
}
}
if (drop_ext_nnz) continue;
for (typename csr<T,distributed>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
size_type dis_i = first_i + i;
size_type dis_j = a.jext2dis_j ((*p).first);
assert_macro (dis_i < dis_nr && dis_j < dis_nc, "invalid matrix");
if (!use_symmetry || dis_j >= dis_i) {
row[q] = dis_i + 1;
col[q] = dis_j + 1;
q++;
}
}
}
}
template<class T>
static
void
csr2mumps_values_dis (const csr<T,sequential>& a, vector<T>& val, bool use_symmetry, bool)
{
// dummy case: for the compiler to be happy
}
template<class T>
static
void
csr2mumps_values_dis (const csr<T,distributed>& a, vector<T>& val, bool use_symmetry, bool drop_ext_nnz)
{
#ifdef TO_CLEAN
std:fill (val.begin(), val.end(), std::numeric_limits<T>::max());
#endif // TO_CLEAN
typedef typename csr<T,distributed>::size_type size_type;
size_type first_i = a.row_ownership().first_index();
typename csr<T,distributed>::const_iterator dia_ia = a.begin();
typename csr<T,distributed>::const_iterator ext_ia = a.ext_begin();
for (size_type i = 0, n = a.nrow(), q = 0; i < n; i++) {
for (typename csr<T,distributed>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
size_type dis_i = first_i + i;
size_type dis_j = first_i + (*p).first;
if (!use_symmetry || dis_j >= dis_i) {
#ifdef TO_CLEAN
check_macro(q < val.size(), "invalid index");
#endif // TO_CLEAN
val[q] = (*p).second;
q++;
}
}
if (drop_ext_nnz) continue;
for (typename csr<T,distributed>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
size_type dis_i = first_i + i;
size_type dis_j = a.jext2dis_j ((*p).first);
if (!use_symmetry || dis_j >= dis_i) {
#ifdef TO_CLEAN
check_macro(q < val.size(), "invalid index");
#endif // TO_CLEAN
val[q] = (*p).second;
q++;
}
}
}
#ifdef TO_CLEAN
T val_min = std::numeric_limits<T>::max(), val_max = 0;
for (size_t q = 0; q < val.size(); q++) {
val_min = std::min (val_min, abs(val[q]));
val_max = std::max (val_max, abs(val[q]));
}
warning_macro ("val_min="<<val_min<<", val_max="<<val_max);
#endif // TO_CLEAN
}
#endif // _RHEOLEF_HAVE_MPI
// =========================================================================
// 2. the class interface
// =========================================================================
template<class T, class M>
solver_mumps_rep<T,M>::solver_mumps_rep (const csr<T,M>& a, const solver_option& opt)
: solver_abstract_rep<T,M>(opt),
_has_mumps_instance(false),
_drop_ext_nnz(false),
_mumps_par(),
_row(),
_col(),
_val(),
_a00(0),
_det()
{
_drop_ext_nnz = is_distributed<M>::value && opt.force_seq;
update_values (a);
}
template<class T, class M>
void
solver_mumps_rep<T,M>::update_values (const csr<T,M>& a)
{
size_type nproc = communicator().size();
if (a.dis_nrow() <= 1) {
// skip 0 or 1 dis_nrow, where mumps core dump
if (a.nrow() == 1) {
_a00 = (*(a.begin()[0])).second;
}
_has_mumps_instance = false;
return;
}
_has_mumps_instance = true;
// ----------------------------------
// step 0 : init a mumps instance
// ----------------------------------
bool use_symmetry = a.is_symmetric();
bool use_definite_positive = a.is_definite_positive();
_mumps_par.par = 1; // use parallel
#ifdef _RHEOLEF_HAVE_MPI
_mumps_par.comm_fortran = MPI_Comm_c2f (a.row_ownership().comm());
#endif // _RHEOLEF_HAVE_MPI
if (use_symmetry && !use_definite_positive) {
// sym invertible but could be undef: could have either < 0 or > 0 eigenvalues
// => leads to similar cpu with s.d.p. but without the asumption of positive or negative :
_mumps_par.sym = 2;
} else if (use_symmetry && use_definite_positive) {
#if 0
// sym def pos
// a) standard approach
// Note: could fail in scalapack when sym def negative
_mumps_par.sym = 1;
// b) second approach: from AmeButGue-2011-mumps, page 15:
// Another approach to suppress numerical pivoting which works with ScaLAPACK
// for both positive definite and negative definite matrices consists in setting SYM=2 and
// CNTL(1)=0.0D0 (recommended strategy).
_mumps_par.sym = 2;
_mumps_par.cntl[1-1] = 0.0;
#endif // 0
// c) third approach: do not use the knowledge of positive or negative, since
// this knowledge is contained in bilinear forms defined by the user
// use only the general symetry:
// tests show CPU is very similar than two previous approaches
// and the user do not need know about positive_definite or negative_definite
_mumps_par.sym = 2;
} else {
// unsym invertible
_mumps_par.sym = 0;
}
_mumps_par.job = -1;
dmumps_c(&_mumps_par);
check_macro (_mumps_par.infog[1-1] == 0, "mumps: an error has occurred");
// ----------------------------------
// step 1 : symbolic factorization
// ----------------------------------
if (base::option().verbose_level != 0) {
_mumps_par.icntl[1-1] = 1; // error output
_mumps_par.icntl[2-1] = 1; // verbose output
_mumps_par.icntl[3-1] = 1; // global infos
_mumps_par.icntl[4-1] = base::option().verbose_level; // verbosity level
// strcpy (_mumps_par.write_problem, "dump_mumps"); // dump sparse structure
} else {
_mumps_par.icntl[1-1] = -1;
_mumps_par.icntl[2-1] = -1;
_mumps_par.icntl[3-1] = -1;
_mumps_par.icntl[4-1] = 0;
}
if (base::option().compute_determinant) {
_mumps_par.icntl[33-1] = 1;
}
_mumps_par.icntl[5-1] = 0; // matrix is in assembled format
#if defined(_RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_SCOTCH)
_mumps_par.icntl[7-1] = 3; // sequential ordering: 3==scotch
#elif defined(_RHEOLEF_HAVE_MUMPS_WITH_PARMETIS) || defined(_RHEOLEF_HAVE_MUMPS_WITH_METIS)
_mumps_par.icntl[7-1] = 5; // sequential ordering: 5==metis
#else // _RHEOLEF_HAVE_MUMPS_WITH_PARMETIS
_mumps_par.icntl[7-1] = 7; // ordering: 7=auto
#endif // _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
_mumps_par.icntl[14-1] = 40; // default 20 seems to be insufficient in some cases
_mumps_par.icntl[29-1] = 0; // 0: auto; 1: ptscotch; 2: parmetis
_mumps_par.icntl[22-1] = 0; // 0: in-core ; 1: out-of-core TODO
if (is_distributed<M>::value) {
_mumps_par.icntl[18-1] = 3; // distributed
#if defined(_RHEOLEF_HAVE_MUMPS_WITH_SCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_METIS)
_mumps_par.icntl[28-1] = 1; // sequential ordering (preferred, more robust)
#elif defined(_RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH) || defined(_RHEOLEF_HAVE_MUMPS_WITH_PARMETIS)
_mumps_par.icntl[28-1] = 2; // parallel ordering (less robust, but faster)
#endif // _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH || _RHEOLEF_HAVE_MUMPS_WITH_PARMETIS
#ifdef _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
_mumps_par.icntl[29-1] = 1; // ptscotch=1
#elif _RHEOLEF_HAVE_MUMPS_WITH_PARMETIS
_mumps_par.icntl[29-1] = 2; // parmetis=2
#else
_mumps_par.icntl[29-1] = 0; // automatic choice
#endif // _RHEOLEF_HAVE_MUMPS_WITH_PTSCOTCH
} else {
_mumps_par.icntl[18-1] = 0; // sequential
}
// load the matrix in mumps
_mumps_par.n = a.dis_nrow();
size_type dia_nnz = ((use_symmetry) ? nnz_dia_upper(a) : a.nnz());
size_type ext_nnz = ((use_symmetry) ? nnz_ext_upper(a) : a.ext_nnz());
size_type dis_nnz = dia_nnz + (_drop_ext_nnz ? 0 : ext_nnz);
if (is_distributed<M>::value) {
#ifdef _RHEOLEF_HAVE_MPI
_row.resize (dis_nnz);
_col.resize (dis_nnz);
csr2mumps_struct_dis (a, _row, _col, use_symmetry, _drop_ext_nnz);
_mumps_par.nz_loc = dis_nnz;
_mumps_par.irn_loc = _row.begin().operator->();
_mumps_par.jcn_loc = _col.begin().operator->();
#endif // _RHEOLEF_HAVE_MPI
} else { // sequential
_row.resize (dia_nnz);
_col.resize (dia_nnz);
csr2mumps_struct_seq (a, _row, _col, use_symmetry);
_mumps_par.nz = dia_nnz;
#if (_RHEOLEF_MUMPS_VERSION_MAJOR >= 5)
_mumps_par.nnz = dia_nnz;
#endif
_mumps_par.irn = _row.begin().operator->();
_mumps_par.jcn = _col.begin().operator->();
}
_mumps_par.job = 1;
dmumps_c(&_mumps_par);
trace_macro ("symbolic: ordering effectively used = " << _mumps_par.infog[7-1]); // 3=scoth, etc
check_macro (_mumps_par.infog[1-1] == 0, "mumps: error infog(1)="<<_mumps_par.infog[1-1]<<" has occurred (see mumps/infog(1) documentation; infog(2)="<<_mumps_par.infog[2-1]<<")");
if (a.dis_nrow() <= 1) {
// skip 0 or 1 dis_nrow, where mumps core dump
if (a.nrow() == 1) {
_a00 = (*(a.begin()[0])).second;
}
return;
}
// ----------------------------------
// step 2 : numeric factorization
// ----------------------------------
if (is_distributed<M>::value) {
#ifdef _RHEOLEF_HAVE_MPI
_val.resize (dis_nnz);
csr2mumps_values_dis (a, _val, use_symmetry, _drop_ext_nnz);
_mumps_par.a_loc = _val.begin().operator->();
#endif // _RHEOLEF_HAVE_MPI
} else { // sequential
_val.resize (dia_nnz);
csr2mumps_values_seq (a, _val, use_symmetry);
_mumps_par.a = _val.begin().operator->();
}
_mumps_par.job = 2;
bool finished = false;
while (!finished && _mumps_par.icntl[14-1] <= 2000) {
dmumps_c(&_mumps_par);
if ((_mumps_par.infog[1-1] == -8 || _mumps_par.infog[1-1] == -9) &&
_mumps_par.infog[2-1] >= 0) {
// not enought working space:
// default 20 seems to be insufficient in some cases (periodic,surfacic,p-laplacian)
_mumps_par.icntl[14-1] += 10;
if (_mumps_par.icntl[14-1] > 100) {
dis_warning_macro ("numeric factorization: increase working space of "<<_mumps_par.icntl[14-1]<< "% and retry...");
}
} else {
finished = true;
}
}
check_macro (_mumps_par.infog[1-1] == 0,
"solver failed: mumps error infog(1)=" <<_mumps_par.infog[1-1]
<< ", infog(2)=" <<_mumps_par.infog[2-1]
<< ", icntl(14)="<<_mumps_par.icntl[14-1]
<< ", icntl(23)="<<_mumps_par.icntl[23-1]);
if (_mumps_par.icntl[33-1] != 0) {
// get determinant from infos:
_det.mantissa = _mumps_par.rinfog[12-1];
// imag part _mumps_par.rinfog[13-1] is zero as matrix is real here
_det.exponant = _mumps_par.infog[34-1];
_det.base = 2;
}
}
template<class T, class M>
vec<T,M>
solver_mumps_rep<T,M>::solve (const vec<T,M>& b) const
{
if (b.dis_size() <= 1) {
// skip 0 or 1 dis_nrow, where mumps core dump
if (b.size() == 1) {
return vec<T,M>(b.ownership(), b[0]/_a00);
} else {
return vec<T,M>(b.ownership());
}
}
// ----------------------------------
// step 3 : define rhs & solve
// ----------------------------------
vec<T,M> x;
vector<T> sol;
vector<int> perm;
size_type sol_size = 0;
vec<T,M> b_host;
T dummy; // mumps faild when zero sizes and 0 pointers...
int idummy;
if (is_distributed<M>::value) {
// 3.1a. merge the rhs b on proc=0 as required by distributed mumps (why ?)
_mumps_par.icntl[21-1] = 1; // 1: solution is distributed
_mumps_par.icntl[10-1] = 0; // nstep of iterative refinement (not in distributed version?)
communicator comm = b.ownership().comm();
distributor host_ownership (b.dis_size(), comm, comm.rank() == 0 ? b.dis_size() : 0);
b_host.resize (host_ownership);
size_type first_i = b.ownership().first_index();
for (size_type i = 0, n = b.size(); i < n; i++) {
size_type dis_i = first_i + i;
b_host.dis_entry(dis_i) = b[i];
}
b_host.dis_entry_assembly();
if (comm.rank() == 0) {
_mumps_par.nrhs = 1;
_mumps_par.rhs = b_host.begin().operator->();
}
sol_size = _mumps_par.info[23-1];
sol.resize (sol_size);
perm.resize (sol_size);
_mumps_par.lsol_loc = sol_size;
_mumps_par.sol_loc = (sol_size > 0) ? sol.begin().operator->() : &dummy;
_mumps_par.isol_loc = (sol_size > 0) ? perm.begin().operator->() : &idummy;
} else { // sequential
// 3.1b. inplace resolution: copy b in x and put x as rhs
_mumps_par.icntl[21-1] = 0; // 0: sol goes on the proc=0, inplace in rhs
_mumps_par.icntl[10-1] = 0; // nstep of iterative refinement (not in distributed version?)
x = b;
_mumps_par.nrhs = 1;
_mumps_par.rhs = x.begin().operator->();
}
// 3.2. run mumps solver
_mumps_par.icntl [9-1] = 1; // solve A*x=b : ctrl=1 ; otherwise solve A^t*x=b TODO
_mumps_par.job = 3;
dmumps_c(&_mumps_par);
check_macro (_mumps_par.infog[1-1] >= 0,
"solver failed: mumps error infog(1)="<<_mumps_par.infog[1-1]
<< ", infog(2)="<<_mumps_par.infog[2-1]);
if (_mumps_par.infog[1-1] > 0) {
warning_macro ("mumps warning infog(1)="<<_mumps_par.infog[1-1]
<< ", infog(2)="<<_mumps_par.infog[2-1]);
}
if (is_distributed<M>::value) {
// 3.3. redistribute the solution as b.ownership
x.resize (b.ownership());
for (size_t i = 0; i < sol_size; i++) {
size_type dis_i = perm[i] - 1;
assert_macro (dis_i < x.dis_size(), "invalid index");
x.dis_entry(dis_i) = sol[i];
}
x.dis_entry_assembly();
}
return x;
}
template<class T, class M>
vec<T,M>
solver_mumps_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
error_macro ("not yet");
return vec<T,M>();
}
template<class T, class M>
solver_mumps_rep<T,M>::~solver_mumps_rep ()
{
if (!_has_mumps_instance) return;
_has_mumps_instance = false;
// ----------------------------------
// step F : close the mumps instance
// ----------------------------------
_mumps_par.job = -2;
dmumps_c(&_mumps_par);
check_macro (_mumps_par.infog[1-1] == 0, "mumps: an error has occurred");
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
// TODO: code is only valid here for T=double
template class solver_mumps_rep<double,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_mumps_rep<double,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
#endif // MUMPS
size_type
field::size_type size_type
Definition branch.cc:430
communicator
see the communicator page for the full documentation
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::solver_abstract_rep
Definition solver.h:191
rheolef::solver_mumps_rep
Definition solver_mumps.h:39
rheolef::solver_mumps_rep::~solver_mumps_rep
~solver_mumps_rep()
Definition solver_mumps.cc:527
rheolef::solver_mumps_rep::update_values
void update_values(const csr< T, M > &a)
Definition solver_mumps.cc:265
rheolef::solver_mumps_rep::size_type
base::size_type size_type
Definition solver_mumps.h:44
rheolef::solver_mumps_rep::solver_mumps_rep
solver_mumps_rep(const csr< T, M > &a, const solver_option &opt=solver_option())
Definition solver_mumps.cc:249
rheolef::solver_mumps_rep::_drop_ext_nnz
bool _drop_ext_nnz
Definition solver_mumps.h:64
rheolef::solver_mumps_rep::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition solver_mumps.cc:521
rheolef::solver_mumps_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition solver_mumps.cc:445
rheolef::solver_option
see the solver_option page for the full documentation
Definition solver_option.h:155
rheolef::solver_option::force_seq
bool force_seq
Definition solver_option.h:173
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
rheolef::vec::resize
void resize(const distributor &ownership, const T &init_val=std::numeric_limits< T >::max())
Definition vec.h:199
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
assert_macro
#define assert_macro(ok_condition, message)
Definition dis_macros.h:113
error_macro
#define error_macro(message)
Definition dis_macros.h:49
dis_warning_macro
#define dis_warning_macro(message)
Definition dis_macros.h:66
warning_macro
#define warning_macro(message)
Definition dis_macros.h:53
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.n
int n
Definition mkgeo_ball.sh:150
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
std
STL namespace.
solver_mumps.h
p
Definition sphere.icc:25
rheolef::is_distributed
Definition communicator.h:70
M
Expr1::memory_type M
Definition vec_expr_v2.h:416