24 #ifndef __PASO_MUMPS_H__
25 #define __PASO_MUMPS_H__
31 #ifdef ESYS_HAVE_MUMPS
33 #if defined(MPI_COMM_WORLD)
34 #undef MPI_COMM_WORLD // breaks mumps_mpi.h, defined in escriptcore/src/EsysMPI.h
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] // macro s.t. indices match documentation
51 #endif // ESYS_HAVE_MUMPS
58 #ifdef ESYS_HAVE_MUMPS
60 #ifdef _WIN32 // workaround for d/zmumps dll clash
61 HINSTANCE h_mumps_c_dll;
63 #endif // ESYS_HAVE_MUMPS
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;
88 #ifdef _WIN32 // workaround for d/zmumps dll clash
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;
96 #endif // ESYS_HAVE_MUMPS
102 #ifdef ESYS_HAVE_MUMPS
104 typedef ZMUMPS_COMPLEX mumps_t;
105 #ifdef _WIN32 // workaround for dmumps/zdmumps dll clash
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;
113 #endif // ESYS_HAVE_MUMPS
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 = "
237 #else // ESYS_HAVE_MUMPS
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;
265 #endif // __PASO_MUMPS_H__