24 #ifndef __PASO_MUMPS_H__
25 #define __PASO_MUMPS_H__
31 #ifdef ESYS_HAVE_MUMPS
33 #if defined(MPI_COMM_WORLD)
36 #include <mumps_mpi.h>
41 #define MUMPS_JOB_INIT -1
42 #define MUMPS_JOB_END -2
43 #define MUMPS_USE_COMM_WORLD -987654
44 #define ICNTL(I) icntl[(I)-1]
58 #ifdef ESYS_HAVE_MUMPS
61 HINSTANCE h_mumps_c_dll;
78 void MUMPS_print_list(
const char* name,
const T* vals,
const int n,
const int max_n=100);
85 #ifdef ESYS_HAVE_MUMPS
87 typedef double mumps_t;
89 typedef HRESULT (CALLBACK* MUMPS_C_FUNC_PTR)(DMUMPS_STRUC_C*);
90 MUMPS_C_FUNC_PTR mumps_c;
91 const char* mumps_lib =
"libdmumps";
92 const char* mumps_proc =
"dmumps_c";
94 void (*mumps_c)(DMUMPS_STRUC_C*) = &dmumps_c;
102 #ifdef ESYS_HAVE_MUMPS
104 typedef ZMUMPS_COMPLEX mumps_t;
106 typedef HRESULT (CALLBACK* MUMPS_C_FUNC_PTR)(ZMUMPS_STRUC_C*);
107 MUMPS_C_FUNC_PTR mumps_c;
108 const char* mumps_lib =
"libzmumps";
109 const char* mumps_proc =
"zmumps_c";
111 void (*mumps_c)(ZMUMPS_STRUC_C*) = &zmumps_c;
117 template <
typename T>
121 #ifdef ESYS_HAVE_MUMPS
125 pt->id.job = MUMPS_JOB_END;
126 pt->mumps_c(&pt->id);
128 FreeLibrary(pt->h_mumps_c_dll);
131 std::string
message = pt->ssExceptMsg.str();
137 MUMPS_INT ierr = MPI_Finalize();
139 std::cout <<
"MUMPS: instance terminated." << std::endl;
148 template <
typename T>
151 #ifdef ESYS_HAVE_MUMPS
153 throw PasoException(
"Paso: MUMPS requires CSR format with index offset 1 and block size 1.");
160 pt->h_mumps_c_dll = LoadLibrary(pt->mumps_lib);
161 if (pt->h_mumps_c_dll == NULL) {
162 std::stringstream ss;
163 ss <<
"Paso: MUMPS LoadLibrary failed - \"" << pt->mumps_lib <<
"\".";
167 if (pt->mumps_c == NULL) {
168 std::stringstream ss;
169 ss <<
"Paso: MUMPS GetProcAddress failed - \"" << pt->mumps_proc <<
"\".";
173 A->solver_p = (
void*) pt;
177 A->pattern->csrToHB();
178 MUMPS_INT n = A->numRows;
179 MUMPS_INT8 nnz = A->pattern->len;
180 MUMPS_INT* irn =
reinterpret_cast<MUMPS_INT*
>(A->pattern->hb_row);
181 MUMPS_INT* jcn =
reinterpret_cast<MUMPS_INT*
>(A->pattern->hb_col);
182 pt->verbose = verbose;
184 std::cout <<
"MUMPS in ===>" << std::endl;
185 std::cout <<
"n = " << n << std::endl;
186 std::cout <<
"nnz = " << nnz << std::endl;
195 std::memcpy(pt->rhs, in, n*
sizeof(T));
197 ierr = MPI_Init(NULL, NULL);
201 pt->id.comm_fortran = MUMPS_USE_COMM_WORLD;
202 pt->id.par = 1; pt->id.sym = 0;
203 pt->id.job = MUMPS_JOB_INIT;
204 pt->mumps_c(&pt->id);
207 pt->id.n = n; pt->id.nnz = nnz;
208 pt->id.irn = irn; pt->id.jcn = jcn;
214 pt->id.ICNTL(1)=-1; pt->id.ICNTL(2)=-1; pt->id.ICNTL(3)=-1; pt->id.ICNTL(4)=0;
219 pt->mumps_c(&pt->id);
220 if (pt->id.infog[0] < 0) {
221 pt->
ssExceptMsg <<
"(PROC " << pt->myid <<
") MUMPS ERROR: INFOG(1)=" << pt->id.infog[0]
222 <<
", INFOG(2)=" << pt->id.infog[1];
224 std::memcpy(out,
reinterpret_cast<T*
>(pt->rhs), n*
sizeof(T));
225 if (pt->id.infog[0] > 0) {
226 std::cout <<
"(PROC " << pt->myid <<
") MUMPS WARNING: INFOG(1)=" << pt->id.infog[0]
227 <<
", INFOG(2)=" << pt->id.infog[1];
230 std::cout <<
"MUMPS out ===>" << std::endl;
232 std::cout <<
"MUMPS: factorization and solve completed (time = "
244 template <
typename T>
247 std::cout << name <<
" = [ ";
248 for (
int i=0; i<n; i++) {
252 std::cout << vals[i];
255 std::cout <<
", ...";
260 std::cout <<
" ]" << std::endl;
#define MPI_COMM_WORLD
Definition: EsysMPI.h:50
#define PASO_MUMPS
Definition: Options.h:57
#define MATRIX_FORMAT_BLK1
Definition: Paso.h:63
#define MATRIX_FORMAT_OFFSET1
Definition: Paso.h:64
PasoException exception class.
Definition: PasoException.h:34
std::complex< real_t > cplx_t
complex data type
Definition: DataTypes.h:55
index_t dim_t
Definition: DataTypes.h:66
double gettime()
returns the current ticks for timing
Definition: EsysMPI.h:192
Definition: BiCGStab.cpp:25
void MUMPS_print_list(const char *name, const T *vals, const int n, const int max_n=100)
Definition: MUMPS.h:245
boost::shared_ptr< SparseMatrix< T > > SparseMatrix_ptr
Definition: SparseMatrix.h:37
std::ostream & operator<<(std::ostream &os, const cplx_t &c)
Definition: MUMPS.cpp:34
void MUMPS_solve(SparseMatrix_ptr< T > A, T *out, T *in, dim_t numRefinements, bool verbose)
calls the solver
Definition: MUMPS.h:149
void MUMPS_free(SparseMatrix< T > *A)
frees any MUMPS related data from the matrix
Definition: MUMPS.h:118
Definition: blocktools.h:70
cplx_t * rhs
Definition: MUMPS.h:101
double * rhs
Definition: MUMPS.h:84
std::stringstream ssExceptMsg
Definition: MUMPS.h:57
bool verbose
Definition: MUMPS.h:56
T * rhs
Definition: MUMPS.h:68
Definition: SparseMatrix.h:45
void * solver_p
pointer to data needed by a solver
Definition: SparseMatrix.h:177