33 extern void dmumps_c(DMUMPS_STRUC_C *mumps_param_ptr);
34 extern void zmumps_c(ZMUMPS_STRUC_C *mumps_param_ptr);
37 #define USE_COMM_WORLD -987654
45 static int find_position(
int *Ai,
int Alen,
int idx)
49 register int lo = 0, hi = Alen - 1, mid;
55 if(idx < Ai[mid]) hi = mid - 1;
56 else if(idx > Ai[mid]) lo = mid + 1;
70 template<
typename Scalar>
71 MumpsMatrix<Scalar>::MumpsMatrix()
82 template<
typename Scalar>
83 MumpsMatrix<Scalar>::~MumpsMatrix()
88 template<
typename Scalar>
89 void MumpsMatrix<Scalar>::alloc()
91 assert(this->pages != NULL);
94 Ap =
new unsigned int[this->size + 1];
95 int aisize = this->get_num_indices();
99 unsigned int i, pos = 0;
100 for (i = 0; i < this->size; i++)
103 pos += sort_and_store_indices(this->pages[i], Ai + pos, Ai + aisize);
107 delete [] this->pages;
110 nnz = Ap[this->size];
112 Ax =
new typename mumps_type<Scalar>::mumps_Scalar[nnz];
113 memset(Ax, 0,
sizeof(Scalar) * nnz);
117 for (
unsigned int i = 0; i < nnz; i++)
124 template<
typename Scalar>
125 void MumpsMatrix<Scalar>::free()
128 delete [] Ap; Ap = NULL;
129 delete [] Ai; Ai = NULL;
130 delete [] Ax; Ax = NULL;
131 delete [] irn; irn = NULL;
132 delete [] jcn; jcn = NULL;
135 inline double mumps_to_Scalar(
double x)
140 inline std::complex<double> mumps_to_Scalar(ZMUMPS_COMPLEX x)
142 return std::complex<double>(x.r, x.i);
145 template<
typename Scalar>
146 Scalar MumpsMatrix<Scalar>::get(
unsigned int m,
unsigned int n)
149 int mid = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
151 if(mid < 0)
return 0.0;
153 if(mid >= 0) mid += Ap[n];
154 return mumps_to_Scalar(Ax[mid]);
157 template<
typename Scalar>
158 void MumpsMatrix<Scalar>::zero()
160 memset(Ax, 0,
sizeof(Scalar) * Ap[this->size]);
163 inline ZMUMPS_COMPLEX& operator +=(ZMUMPS_COMPLEX &a, std::complex<double> b)
170 template<
typename Scalar>
171 void MumpsMatrix<Scalar>::add(
unsigned int m,
unsigned int n, Scalar v)
176 int pos = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
182 #pragma omp critical (MumpsMatrix_add)
188 template<
typename Scalar>
189 void MumpsMatrix<Scalar>::add(
unsigned int m,
unsigned int n, Scalar **mat,
int *rows,
int *cols)
191 for (
unsigned int i = 0; i < m; i++)
192 for (
unsigned int j = 0; j < n; j++)
193 if(rows[i] >= 0 && cols[j] >= 0)
194 add(rows[i], cols[j], mat[i][j]);
199 template<
typename Scalar>
200 void MumpsMatrix<Scalar>::add_to_diagonal(Scalar v)
202 for (
unsigned int i = 0; i < this->size; i++)
210 template<
typename Scalar>
211 bool MumpsMatrix<Scalar>::dump(FILE *file,
const char *var_name,
EMatrixDumpFormat fmt,
char* number_format)
217 fprintf(file,
"%d\n", this->size);
218 fprintf(file,
"%d\n", nnz);
219 for (
unsigned int i = 0; i < nnz; i++)
221 fprintf(file,
"%d %d ", irn[i], jcn[i]);
222 Hermes::Helpers::fprint_num(file, mumps_to_Scalar(Ax[i]));
228 fprintf(file,
"%% Size: %dx%d\n%% Nonzeros: %d\ntemp = zeros(%d, 3);\ntemp =[\n", this->size, this->size, Ap[this->size], Ap[this->size]);
229 for (
unsigned int j = 0; j < this->size; j++)
230 for (
unsigned int i = Ap[j]; i < Ap[j + 1]; i++)
232 fprintf(file,
"%d %d ", Ai[i] + 1, j + 1);
233 Hermes::Helpers::fprint_num(file, mumps_to_Scalar(Ax[i]));
236 fprintf(file,
"];\n%s = spconvert(temp);\n", var_name);
242 this->hermes_fwrite(
"HERMESX\001", 1, 8, file);
243 int ssize =
sizeof(Scalar);
244 this->hermes_fwrite(&ssize,
sizeof(
int), 1, file);
245 this->hermes_fwrite(&this->size,
sizeof(
int), 1, file);
246 this->hermes_fwrite(&nnz,
sizeof(
int), 1, file);
247 this->hermes_fwrite(Ap,
sizeof(
int), this->size + 1, file);
248 this->hermes_fwrite(Ai,
sizeof(
int), nnz, file);
249 this->hermes_fwrite(Ax,
sizeof(Scalar), nnz, file);
258 template<
typename Scalar>
259 unsigned int MumpsMatrix<Scalar>::get_matrix_size()
const
264 template<
typename Scalar>
265 unsigned int MumpsMatrix<Scalar>::get_nnz()
const
270 template<
typename Scalar>
271 double MumpsMatrix<Scalar>::get_fill_in()
const
273 return Ap[this->size] / (double) (this->size * this->size);
276 template<
typename Scalar>
277 void MumpsMatrix<Scalar>::add_matrix(MumpsMatrix<Scalar>* mat)
279 add_as_block(0, 0, mat);
282 template<
typename Scalar>
283 void MumpsMatrix<Scalar>::add_to_diagonal_blocks(
int num_stages, MumpsMatrix<Scalar>* mat)
285 int ndof = mat->get_size();
286 if(this->get_size() != (
unsigned int) num_stages * ndof)
289 for (
int i = 0; i < num_stages; i++)
291 this->add_as_block(ndof*i, ndof*i, mat);
295 template<
typename Scalar>
296 void MumpsMatrix<Scalar>::add_sparse_to_diagonal_blocks(
int num_stages, SparseMatrix<Scalar>* mat)
298 add_to_diagonal_blocks(num_stages, static_cast<MumpsMatrix*>(mat));
301 inline ZMUMPS_COMPLEX& operator +=(ZMUMPS_COMPLEX &a, ZMUMPS_COMPLEX b)
308 template<
typename Scalar>
309 void MumpsMatrix<Scalar>::add_as_block(
unsigned int i,
unsigned int j, MumpsMatrix<Scalar>* mat)
312 for (
unsigned int col = 0;col<mat->get_size();col++)
314 for (
unsigned int n = mat->Ap[col];n<mat->Ap[col + 1];n++)
316 idx = find_position(Ai + Ap[col + j], Ap[col + 1 + j] - Ap[col + j], mat->Ai[n] + i);
320 Ax[idx] +=mat->Ax[n];
326 template<
typename Scalar>
327 void MumpsMatrix<Scalar>::multiply_with_vector(Scalar* vector_in, Scalar* vector_out)
329 for(
unsigned int i = 0;i<this->size;i++)
334 for (
unsigned int i = 0;i<nnz;i++)
336 a = mumps_to_Scalar(Ax[i]);
337 vector_out[jcn[i]-1] +=vector_in[irn[i]-1]*a;
342 void MumpsMatrix<double>::multiply_with_Scalar(
double value)
345 for(
int i = 0;i<n;i++)
352 void MumpsMatrix<std::complex<double> >::multiply_with_Scalar(std::complex<double> value)
355 std::complex<double> a;
356 for(
int i = 0;i<n;i++)
358 a = std::complex<double>(Ax[i].r, Ax[i].i);
365 inline void mumps_assign_Scalar(ZMUMPS_COMPLEX & a, std::complex<double> b)
371 inline void mumps_assign_Scalar(
double & a,
double b)
377 template<
typename Scalar>
378 void MumpsMatrix<Scalar>::create(
unsigned int size,
unsigned int nnz,
int* ap,
int* ai, Scalar* ax)
382 this->Ap =
new unsigned int[this->size + 1]; assert(this->Ap != NULL);
383 this->Ai =
new int[nnz]; assert(this->Ai != NULL);
384 this->Ax =
new typename mumps_type<Scalar>::mumps_Scalar[nnz]; assert(this->Ax != NULL);
385 irn =
new int[nnz]; assert(this->irn !=NULL);
386 jcn =
new int[nnz]; assert(this->jcn !=NULL);
388 for (
unsigned int i = 0; i < this->size; i++)
391 for (
int j = ap[i];j<ap[i + 1];j++) jcn[j] = i;
393 this->Ap[this->size] = ap[this->size];
394 for (
unsigned int i = 0; i < nnz; i++)
396 mumps_assign_Scalar(this->Ax[i], ax[i]);
402 template<
typename Scalar>
403 MumpsMatrix<Scalar>* MumpsMatrix<Scalar>::duplicate()
405 MumpsMatrix<Scalar> * nmat =
new MumpsMatrix<Scalar>();
408 nmat->size = this->size;
409 nmat->Ap =
new unsigned int[this->size + 1]; assert(nmat->Ap != NULL);
410 nmat->Ai =
new int[nnz]; assert(nmat->Ai != NULL);
411 nmat->Ax =
new typename mumps_type<Scalar>::mumps_Scalar[nnz]; assert(nmat->Ax != NULL);
412 nmat->irn =
new int[nnz]; assert(nmat->irn !=NULL);
413 nmat->jcn =
new int[nnz]; assert(nmat->jcn !=NULL);
414 for (
unsigned int i = 0;i<nnz;i++)
418 nmat->irn[i] = irn[i];
419 nmat->jcn[i] = jcn[i];
421 for (
unsigned int i = 0;i<this->size + 1;i++)
430 template<
typename Scalar>
431 MumpsVector<Scalar>::MumpsVector()
437 template<
typename Scalar>
438 MumpsVector<Scalar>::~MumpsVector()
443 template<
typename Scalar>
444 void MumpsVector<Scalar>::alloc(
unsigned int n)
452 template<
typename Scalar>
453 void MumpsVector<Scalar>::change_sign()
455 for (
unsigned int i = 0; i < this->size; i++) v[i] *= -1.;
458 template<
typename Scalar>
459 void MumpsVector<Scalar>::zero()
461 memset(v, 0, this->size *
sizeof(Scalar));
464 template<
typename Scalar>
465 void MumpsVector<Scalar>::free()
472 template<
typename Scalar>
473 void MumpsVector<Scalar>::set(
unsigned int idx, Scalar y)
478 template<
typename Scalar>
479 void MumpsVector<Scalar>::add(
unsigned int idx, Scalar y)
481 #pragma omp critical (MumpsVector_add)
485 template<
typename Scalar>
486 void MumpsVector<Scalar>::add(
unsigned int n,
unsigned int *idx, Scalar *y)
488 for (
unsigned int i = 0; i < n; i++)
494 template<
typename Scalar>
495 Scalar MumpsVector<Scalar>::get(
unsigned int idx)
500 template<
typename Scalar>
501 void MumpsVector<Scalar>::extract(Scalar *v)
const
503 memcpy(v, this->v, this->size *
sizeof(Scalar));
506 template<
typename Scalar>
507 void MumpsVector<Scalar>::add_vector(Vector<Scalar>* vec)
509 assert(this->length() == vec->length());
510 for (
unsigned int i = 0; i < this->length(); i++) this->add(i, vec->get(i));
513 template<
typename Scalar>
514 void MumpsVector<Scalar>::add_vector(Scalar* vec)
516 for (
unsigned int i = 0; i < this->length(); i++) this->add(i, vec[i]);
519 template<
typename Scalar>
520 bool MumpsVector<Scalar>::dump(FILE *file,
const char *var_name,
EMatrixDumpFormat fmt,
char* number_format)
525 for (
unsigned int i = 0; i < this->size; i++)
527 Hermes::Helpers::fprint_num(file, v[i], number_format);
534 fprintf(file,
"%% Size: %dx1\n%s =[\n", this->size, var_name);
535 for (
unsigned int i = 0; i < this->size; i++)
537 Hermes::Helpers::fprint_num(file, v[i], number_format);
540 fprintf(file,
" ];\n");
545 this->hermes_fwrite(
"HERMESR\001", 1, 8, file);
546 int ssize =
sizeof(Scalar);
547 this->hermes_fwrite(&ssize,
sizeof(
int), 1, file);
548 this->hermes_fwrite(&this->size,
sizeof(
int), 1, file);
549 this->hermes_fwrite(v,
sizeof(Scalar), this->size, file);
558 template class HERMES_API MumpsMatrix<double>;
559 template class HERMES_API MumpsMatrix<std::complex<double> >;
560 template class HERMES_API MumpsVector<double>;
561 template class HERMES_API MumpsVector<std::complex<double> >;
566 #define ICNTL(I) icntl[(I)-1]
567 #define MUMPS_INFO(param, I) (param).infog[(I)-1]
568 #define INFOG(I) infog[(I)-1]
573 #define JOB_ANALYZE_FACTORIZE_SOLVE 6
574 #define JOB_FACTORIZE_SOLVE 5
578 void MumpsSolver<double>::mumps_c(mumps_type<double>::mumps_struct * param)
584 void MumpsSolver<std::complex<double> >::mumps_c(mumps_type<std::complex<double> >::mumps_struct * param)
589 template<
typename Scalar>
590 bool MumpsSolver<Scalar>::check_status()
592 switch (param.INFOG(1))
595 case -1: this->warn(
"Error occured on processor %d", MUMPS_INFO(param, 2));
break;
597 default: this->warn(
"INFOG(1) = %d", param.INFOG(1));
break;
602 template<
typename Scalar>
603 bool MumpsSolver<Scalar>::reinit()
613 param.job = JOB_INIT;
616 param.comm_fortran = USE_COMM_WORLD;
619 inited = check_status();
643 template<
typename Scalar>
644 MumpsSolver<Scalar>::MumpsSolver(MumpsMatrix<Scalar> *m, MumpsVector<Scalar> *rhs) :
645 DirectSolver<Scalar>(), m(m), rhs(rhs)
652 param.INFOG(33) = -999;
656 template<
typename Scalar>
657 MumpsSolver<Scalar>::~MumpsSolver()
666 if(param.rhs != NULL)
delete [] param.rhs;
669 template<
typename Scalar>
670 int MumpsSolver<Scalar>::get_matrix_size()
675 template<
typename Scalar>
676 bool MumpsSolver<Scalar>::solve()
687 if( !setup_factorization() )
693 param.rhs =
new typename mumps_type<Scalar>::mumps_Scalar[m->size];
694 memcpy(param.rhs, rhs->v, m->size *
sizeof(Scalar));
699 ret = check_status();
704 this->sln =
new Scalar[m->size];
705 for (
unsigned int i = 0; i < rhs->size; i++)
706 this->sln[i] = mumps_to_Scalar(param.rhs[i]);
710 this->time = this->accumulated();
718 template<
typename Scalar>
719 bool MumpsSolver<Scalar>::setup_factorization()
723 int eff_fact_scheme = this->factorization_scheme;
725 if( this->factorization_scheme == HERMES_REUSE_MATRIX_REORDERING ||
726 this->factorization_scheme == HERMES_REUSE_FACTORIZATION_COMPLETELY )
729 switch (eff_fact_scheme)
738 param.job = JOB_ANALYZE_FACTORIZE_SOLVE;
747 param.job = JOB_FACTORIZE_SOLVE;
754 if(param.INFOG(33) != -2)
758 param.job = JOB_ANALYZE_FACTORIZE_SOLVE;
763 param.job = JOB_FACTORIZE_SOLVE;
767 param.job = JOB_SOLVE;
774 template class HERMES_API MumpsSolver<double>;
775 template class HERMES_API MumpsSolver<std::complex<double> >;