241class AlephMatrixHypre
249 : IAlephMatrix(tm, kernel, index)
250 , m_hypre_ijmatrix(0)
252 debug() <<
"[AlephMatrixHypre] new AlephMatrixHypre";
257 debug() <<
"[~AlephMatrixHypre]";
258 if (m_hypre_ijmatrix)
259 HYPRE_IJMatrixDestroy(m_hypre_ijmatrix);
274 void AlephMatrixCreate(
void)
276 debug() <<
"[AlephMatrixHypre::AlephMatrixCreate] HYPRE MatrixCreate idx:" << m_index;
280 for (
int iCpu = 0; iCpu < m_kernel->size(); ++iCpu) {
281 if (m_kernel->rank() != m_kernel->solverRanks(m_index)[iCpu])
284 ilower = m_kernel->topology()->gathered_nb_row(iCpu);
285 iupper = m_kernel->topology()->gathered_nb_row(iCpu + 1) - 1;
287 debug() <<
"[AlephMatrixHypre::AlephMatrixCreate] ilower=" << ilower <<
", iupper=" << iupper;
291 debug() <<
"[AlephMatrixHypre::AlephMatrixCreate] jlower=" << jlower <<
", jupper=" << jupper;
293 hypreCheck(
"HYPRE_IJMatrixCreate",
294 HYPRE_IJMatrixCreate(MPI_COMM_SUB,
299 debug() <<
"[AlephMatrixHypre::AlephMatrixCreate] HYPRE IJMatrixSetObjectType";
300 HYPRE_IJMatrixSetObjectType(m_hypre_ijmatrix, HYPRE_PARCSR);
301 debug() <<
"[AlephMatrixHypre::AlephMatrixCreate] HYPRE IJMatrixSetRowSizes";
302 HYPRE_IJMatrixSetRowSizes(m_hypre_ijmatrix, m_kernel->topology()->gathered_nb_row_elements().data());
303 debug() <<
"[AlephMatrixHypre::AlephMatrixCreate] HYPRE IJMatrixInitialize";
304 HYPRE_IJMatrixInitialize(m_hypre_ijmatrix);
305 HYPRE_IJMatrixGetObject(m_hypre_ijmatrix, &
object);
306 m_hypre_parmatrix = (HYPRE_ParCSRMatrix)
object;
311 void AlephMatrixSetFilled(
bool) {}
315 int AlephMatrixAssemble(
void)
317 debug() <<
"[AlephMatrixHypre::AlephMatrixAssemble]";
318 hypreCheck(
"HYPRE_IJMatrixAssemble",
319 HYPRE_IJMatrixAssemble(m_hypre_ijmatrix));
325 void AlephMatrixFill(
int size, HYPRE_Int* rows, HYPRE_Int* cols,
double* values)
327 debug() <<
"[AlephMatrixHypre::AlephMatrixFill] size=" << size;
329 HYPRE_Int col[1] = { 1 };
330 for (
int i = 0; i < size; i++) {
331 rtn += HYPRE_IJMatrixSetValues(m_hypre_ijmatrix, 1, col, &rows[i], &cols[i], &values[i]);
333 hypreCheck(
"HYPRE_IJMatrixSetValues", rtn);
335 debug() <<
"[AlephMatrixHypre::AlephMatrixFill] done";
347 HYPRE_ClearAllErrors();
348 const bool convergence_analyse = params->convergenceAnalyse();
351 const Real res0 = b->norm_max();
353 if (convergence_analyse)
354 info() <<
"analyse convergence : norme max du second membre res0 : " << res0;
356 const Real considered_as_null = params->minRHSNorm();
357 if (res0 < considered_as_null) {
358 HYPRE_ParVectorSetConstantValues(x->m_hypre_parvector, 0.0);
359 residual_norm[0] = res0;
360 if (convergence_analyse)
361 info() <<
"analyse convergence : le second membre du système linéaire est inférieur à : " << considered_as_null;
365 if (params->xoUser()) {
372 HYPRE_ParCSRMatrixMatvec(1.0, m_hypre_parmatrix, x->m_hypre_parvector, 0., tmp->m_hypre_parvector);
373 HYPRE_ParVectorAxpy(-1.0, b->m_hypre_parvector, tmp->m_hypre_parvector);
375 const Real residu = tmp->norm_max();
378 if (residu < considered_as_null) {
379 if (convergence_analyse) {
380 info() <<
"analyse convergence : |Ax0-b| est inférieur à " << considered_as_null;
381 info() <<
"analyse convergence : x0 est déjà solution du système.";
383 residual_norm[0] = residu;
387 const Real relative_error = residu / res0;
388 if (convergence_analyse)
389 info() <<
"analyse convergence : résidu initial : " << residu
390 <<
" --- residu relatif initial (residu/res0) : " << residu / res0;
392 if (relative_error < (params->epsilon())) {
393 if (convergence_analyse)
394 info() <<
"analyse convergence : X est déjà solution du système";
395 residual_norm[0] = residu;
411 solver_param->setAmgCoarseningMethod(TypesSolver::AMG_COARSENING_AUTO);
412 const String func_name(
"SolverMatrixHypre::solve");
419 HYPRE_IJVector solution = ximpl->m_hypre_ijvector;
420 HYPRE_IJVector RHS = bimpl->m_hypre_ijvector;
423 HYPRE_IJMatrixGetObject(m_hypre_ijmatrix, &
object);
424 HYPRE_ParCSRMatrix M = (HYPRE_ParCSRMatrix)
object;
425 HYPRE_IJVectorGetObject(solution, &
object);
426 HYPRE_ParVector X = (HYPRE_ParVector)
object;
427 HYPRE_IJVectorGetObject(RHS, &
object);
428 HYPRE_ParVector B = (HYPRE_ParVector)
object;
435 if (isAlreadySolved(ximpl2, bimpl2, timpl2, residual_norm, solver_param)) {
436 ItacRegion(isAlreadySolved, AlephMatrixHypre);
437 debug() <<
"[AlephMatrixHypre::AlephMatrixSolve] isAlreadySolved !";
442 TypesSolver::ePreconditionerMethod preconditioner_method = solver_param->precond();
444 TypesSolver::eSolverMethod solver_method = solver_param->method();
447 HYPRE_Solver solver = 0;
449 switch (solver_method) {
450 case TypesSolver::PCG:
451 initSolverPCG(solver_param, solver);
453 case TypesSolver::BiCGStab:
454 initSolverBiCGStab(solver_param, solver);
456 case TypesSolver::GMRES:
457 initSolverGMRES(solver_param, solver);
464 HYPRE_Solver precond = 0;
466 switch (preconditioner_method) {
467 case TypesSolver::NONE:
469 case TypesSolver::DIAGONAL:
470 setDiagonalPreconditioner(solver_method, solver, precond);
472 case TypesSolver::ILU:
473 setILUPreconditioner(solver_method, solver, precond);
475 case TypesSolver::SPAIstat:
476 setSpaiStatPreconditioner(solver_method, solver, solver_param, precond);
478 case TypesSolver::AMG:
479 setAMGPreconditioner(solver_method, solver, solver_param, precond);
481 case TypesSolver::AINV:
483 case TypesSolver::SPAIdyn:
485 case TypesSolver::ILUp:
487 case TypesSolver::IC:
489 case TypesSolver::POLY:
496 HYPRE_Int iteration = 0;
497 double residue = 0.0;
499 switch (solver_method) {
500 case TypesSolver::PCG:
501 ierr = solvePCG(solver_param, solver, M, B, X, iteration, residue);
503 case TypesSolver::BiCGStab:
504 ierr = solveBiCGStab(solver, M, B, X, iteration, residue);
506 case TypesSolver::GMRES:
507 ierr = solveGMRES(solver, M, B, X, iteration, residue);
513 nb_iteration =
static_cast<Integer>(iteration);
514 residual_norm[0] =
static_cast<Real>(residue);
529 switch (preconditioner_method) {
530 case TypesSolver::NONE:
532 case TypesSolver::DIAGONAL:
534 case TypesSolver::ILU:
535 HYPRE_ParCSRPilutDestroy(precond);
537 case TypesSolver::SPAIstat:
538 HYPRE_ParCSRParaSailsDestroy(precond);
540 case TypesSolver::AMG:
541 HYPRE_BoomerAMGDestroy(precond);
547 if (iteration == solver_param->maxIter() && solver_param->stopErrorStrategy()) {
548 info() <<
"\n============================================================";
549 info() <<
"\nCette erreur est retournée après " << iteration <<
"\n";
550 info() <<
"\nOn a atteind le nombre max d'itérations du solveur.";
551 info() <<
"\nIl est possible de demander au code de ne pas tenir compte de cette erreur.";
552 info() <<
"\nVoir la documentation du jeu de données concernant le service solveur.";
553 info() <<
"\n======================================================";
554 throw Exception(
"AlephMatrixHypre::Solve",
"On a atteind le nombre max d'itérations du solveur");
561 void writeToFile(
const String filename)
563 HYPRE_IJMatrixPrint(m_hypre_ijmatrix, filename.
localstr());
568 void initSolverPCG(
const AlephParams* solver_param, HYPRE_Solver& solver)
570 const String func_name =
"SolverMatrixHypre::initSolverPCG";
571 double epsilon = solver_param->epsilon();
572 int max_it = solver_param->maxIter();
573 int output_level = solver_param->getOutputLevel();
575 HYPRE_ParCSRPCGCreate(MPI_COMM_SUB, &solver);
576 HYPRE_ParCSRPCGSetMaxIter(solver, max_it);
577 HYPRE_ParCSRPCGSetTol(solver, epsilon);
578 HYPRE_ParCSRPCGSetTwoNorm(solver, 1);
579 HYPRE_ParCSRPCGSetPrintLevel(solver, output_level);
584 void initSolverBiCGStab(
const AlephParams* solver_param, HYPRE_Solver& solver)
586 const String func_name =
"SolverMatrixHypre::initSolverBiCGStab";
587 double epsilon = solver_param->epsilon();
588 int max_it = solver_param->maxIter();
589 int output_level = solver_param->getOutputLevel();
591 HYPRE_ParCSRBiCGSTABCreate(MPI_COMM_SUB, &solver);
592 HYPRE_ParCSRBiCGSTABSetMaxIter(solver, max_it);
593 HYPRE_ParCSRBiCGSTABSetTol(solver, epsilon);
594 HYPRE_ParCSRBiCGSTABSetPrintLevel(solver, output_level);
599 void initSolverGMRES(
const AlephParams* solver_param, HYPRE_Solver& solver)
601 const String func_name =
"SolverMatrixHypre::initSolverGMRES";
602 double epsilon = solver_param->epsilon();
603 int max_it = solver_param->maxIter();
604 int output_level = solver_param->getOutputLevel();
606 HYPRE_ParCSRGMRESCreate(MPI_COMM_SUB, &solver);
607 const int krylov_dim = 20;
608 HYPRE_ParCSRGMRESSetKDim(solver, krylov_dim);
609 HYPRE_ParCSRGMRESSetMaxIter(solver, max_it);
610 HYPRE_ParCSRGMRESSetTol(solver, epsilon);
611 HYPRE_ParCSRGMRESSetPrintLevel(solver, output_level);
616 void setDiagonalPreconditioner(
const TypesSolver::eSolverMethod solver_method,
617 const HYPRE_Solver& solver,
618 HYPRE_Solver& precond)
620 const String func_name =
"SolverMatrixHypre::setDiagonalPreconditioner";
621 switch (solver_method) {
622 case TypesSolver::PCG:
623 HYPRE_ParCSRPCGSetPrecond(solver,
624 HYPRE_ParCSRDiagScale,
625 HYPRE_ParCSRDiagScaleSetup,
628 case TypesSolver::BiCGStab:
629 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
630 HYPRE_ParCSRDiagScale,
631 HYPRE_ParCSRDiagScaleSetup,
634 case TypesSolver::GMRES:
635 HYPRE_ParCSRGMRESSetPrecond(solver,
636 HYPRE_ParCSRDiagScale,
637 HYPRE_ParCSRDiagScaleSetup,
641 throw ArgumentException(func_name,
"solveur inconnu pour preconditionnement 'Diagonal'");
647 void setILUPreconditioner(
const TypesSolver::eSolverMethod solver_method,
648 const HYPRE_Solver& solver,
649 HYPRE_Solver& precond)
651 const String func_name =
"SolverMatrixHypre::setILUPreconditioner";
652 switch (solver_method) {
653 case TypesSolver::PCG:
654 throw ArgumentException(func_name,
"solveur PCG indisponible avec le preconditionnement 'ILU'");
656 case TypesSolver::BiCGStab:
657 HYPRE_ParCSRPilutCreate(MPI_COMM_SUB, &precond);
658 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
659 HYPRE_ParCSRPilutSolve,
660 HYPRE_ParCSRPilutSetup,
663 case TypesSolver::GMRES:
664 HYPRE_ParCSRPilutCreate(MPI_COMM_SUB,
666 HYPRE_ParCSRGMRESSetPrecond(solver,
667 HYPRE_ParCSRPilutSolve,
668 HYPRE_ParCSRPilutSetup,
672 throw ArgumentException(func_name,
"solveur inconnu pour preconditionnement ILU\n");
678 void setSpaiStatPreconditioner(
const TypesSolver::eSolverMethod solver_method,
679 const HYPRE_Solver& solver,
681 HYPRE_Solver& precond)
683 HYPRE_ParCSRParaSailsCreate(MPI_COMM_SUB, &precond);
684 double alpha = solver_param->alpha();
685 int gamma = solver_param->gamma();
690 HYPRE_ParCSRParaSailsSetParams(precond, alpha, gamma);
691 switch (solver_method) {
692 case TypesSolver::PCG:
693 HYPRE_ParCSRPCGSetPrecond(solver, HYPRE_ParCSRParaSailsSolve, HYPRE_ParCSRParaSailsSetup, precond);
695 case TypesSolver::BiCGStab:
696 throw ArgumentException(
"AlephMatrixHypre::setSpaiStatPreconditioner",
"solveur 'BiCGStab' invalide pour preconditionnement 'SPAIstat'");
698 case TypesSolver::GMRES:
700 HYPRE_ParCSRParaSailsSetSym(precond, 0);
701 HYPRE_ParCSRGMRESSetPrecond(solver, HYPRE_ParaSailsSolve, HYPRE_ParaSailsSetup, precond);
704 throw ArgumentException(
"AlephMatrixHypre::setSpaiStatPreconditioner",
"solveur inconnu pour preconditionnement 'SPAIstat'\n");
711 void setAMGPreconditioner(
const TypesSolver::eSolverMethod solver_method,
712 const HYPRE_Solver& solver,
714 HYPRE_Solver& precond)
718 double trunc_factor = 0.1;
719 int cycle_type = solver_param->getAmgCycle();
720 int coarsen_type = solver_param->amgCoarseningMethod();
722 int relax_default = 3;
728 int measure_type = 1;
729 double max_row_sum = 1.0;
732 const int gamma = solver_param->gamma();
736 double strong_threshold = 0.1;
737 const double alpha = solver_param->alpha();
739 strong_threshold = alpha;
741 Integer output_level = solver_param->getOutputLevel();
743 HYPRE_Int* num_grid_sweeps = _allocHypre<HYPRE_Int>(4);
744 HYPRE_Int* grid_relax_type = _allocHypre<HYPRE_Int>(4);
745 HYPRE_Int** grid_relax_points = _allocHypre<HYPRE_Int*>(4);
746 double* relax_weight = _allocHypre<double>(max_levels);
748 for (
int i = 0; i < max_levels; i++)
749 relax_weight[i] = 1.0;
751 if (coarsen_type == 5) {
753 num_grid_sweeps[0] = 3;
754 grid_relax_type[0] = relax_default;
755 grid_relax_points[0] = _allocHypre<HYPRE_Int>(3);
756 grid_relax_points[0][0] = -2;
757 grid_relax_points[0][1] = -1;
758 grid_relax_points[0][2] = 1;
761 num_grid_sweeps[1] = 4;
762 grid_relax_type[1] = relax_default;
763 grid_relax_points[1] = _callocHypre<HYPRE_Int>(4);
764 grid_relax_points[1][0] = -1;
765 grid_relax_points[1][1] = 1;
766 grid_relax_points[1][2] = -2;
767 grid_relax_points[1][3] = -2;
770 num_grid_sweeps[2] = 4;
771 grid_relax_type[2] = relax_default;
772 grid_relax_points[2] = _allocHypre<HYPRE_Int>(4);
773 grid_relax_points[2][0] = -2;
774 grid_relax_points[2][1] = -2;
775 grid_relax_points[2][2] = 1;
776 grid_relax_points[2][3] = -1;
780 num_grid_sweeps[0] = 2 * num_sweep;
781 grid_relax_type[0] = relax_default;
782 grid_relax_points[0] = _allocHypre<HYPRE_Int>(2 * num_sweep);
783 for (
int i = 0; i < 2 * num_sweep; i += 2) {
784 grid_relax_points[0][i] = -1;
785 grid_relax_points[0][i + 1] = 1;
789 num_grid_sweeps[1] = 2 * num_sweep;
790 grid_relax_type[1] = relax_default;
791 grid_relax_points[1] = _allocHypre<HYPRE_Int>(2 * num_sweep);
792 for (
int i = 0; i < 2 * num_sweep; i += 2) {
793 grid_relax_points[1][i] = -1;
794 grid_relax_points[1][i + 1] = 1;
798 num_grid_sweeps[2] = 2 * num_sweep;
799 grid_relax_type[2] = relax_default;
800 grid_relax_points[2] = _allocHypre<HYPRE_Int>(2 * num_sweep);
801 for (
int i = 0; i < 2 * num_sweep; i += 2) {
802 grid_relax_points[2][i] = -1;
803 grid_relax_points[2][i + 1] = 1;
808 num_grid_sweeps[3] = 1;
809 grid_relax_type[3] = 9;
810 grid_relax_points[3] = _allocHypre<HYPRE_Int>(1);
811 grid_relax_points[3][0] = 0;
815 HYPRE_BoomerAMGCreate(&precond);
816 HYPRE_BoomerAMGSetPrintLevel(precond, output_level);
817 HYPRE_BoomerAMGSetCoarsenType(precond, (hybrid * coarsen_type));
818 HYPRE_BoomerAMGSetMeasureType(precond, measure_type);
819 HYPRE_BoomerAMGSetStrongThreshold(precond, strong_threshold);
820 HYPRE_BoomerAMGSetTruncFactor(precond, trunc_factor);
821 HYPRE_BoomerAMGSetMaxIter(precond, 1);
822 HYPRE_BoomerAMGSetCycleType(precond, cycle_type);
823 HYPRE_BoomerAMGSetNumGridSweeps(precond, num_grid_sweeps);
824 HYPRE_BoomerAMGSetGridRelaxType(precond, grid_relax_type);
825 HYPRE_BoomerAMGSetRelaxWeight(precond, relax_weight);
826 HYPRE_BoomerAMGSetGridRelaxPoints(precond, grid_relax_points);
827 HYPRE_BoomerAMGSetTol(precond, 0.0);
828 HYPRE_BoomerAMGSetMaxLevels(precond, max_levels);
829 HYPRE_BoomerAMGSetMaxRowSum(precond, max_row_sum);
831 switch (solver_method) {
832 case TypesSolver::PCG:
833 HYPRE_ParCSRPCGSetPrecond(solver,
834 HYPRE_BoomerAMGSolve,
835 HYPRE_BoomerAMGSetup,
838 case TypesSolver::BiCGStab:
839 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
840 HYPRE_BoomerAMGSolve,
841 HYPRE_BoomerAMGSetup,
844 case TypesSolver::GMRES:
845 HYPRE_ParCSRGMRESSetPrecond(solver,
846 HYPRE_BoomerAMGSolve,
847 HYPRE_BoomerAMGSetup,
851 throw ArgumentException(
"AlephMatrixHypre::setAMGPreconditioner",
"solveur inconnu pour preconditionnement 'AMG'\n");
858 HYPRE_Solver& solver,
859 HYPRE_ParCSRMatrix& M,
862 HYPRE_Int& iteration,
865 const String func_name =
"SolverMatrixHypre::solvePCG";
866 const bool xo = solver_param->xoUser();
870 HYPRE_ParVectorSetConstantValues(X, 0.0);
871 HYPRE_ParCSRPCGSetup(solver, M, B, X);
872 HYPRE_ParCSRPCGSolve(solver, M, B, X);
873 HYPRE_ParCSRPCGGetNumIterations(solver, &iteration);
874 HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(solver, &residue);
876 HYPRE_Int converged = 0;
877 HYPRE_PCGGetConverged(solver, &converged);
878 error |= (!converged);
880 HYPRE_ParCSRPCGDestroy(solver);
887 bool solveBiCGStab(HYPRE_Solver& solver,
888 HYPRE_ParCSRMatrix& M,
891 HYPRE_Int& iteration,
894 const String func_name =
"SolverMatrixHypre::solveBiCGStab";
896 HYPRE_ParVectorSetRandomValues(X, 775);
897 HYPRE_ParCSRBiCGSTABSetup(solver, M, B, X);
898 HYPRE_ParCSRBiCGSTABSolve(solver, M, B, X);
899 HYPRE_ParCSRBiCGSTABGetNumIterations(solver, &iteration);
900 HYPRE_ParCSRBiCGSTABGetFinalRelativeResidualNorm(solver, &residue);
902 HYPRE_Int converged = 0;
903 hypre_BiCGSTABGetConverged(solver, &converged);
904 error |= (!converged);
906 HYPRE_ParCSRBiCGSTABDestroy(solver);
913 bool solveGMRES(HYPRE_Solver& solver,
914 HYPRE_ParCSRMatrix& M,
917 HYPRE_Int& iteration,
920 const String func_name =
"SolverMatrixHypre::solveGMRES";
922 HYPRE_ParCSRGMRESSetup(solver, M, B, X);
923 HYPRE_ParCSRGMRESSolve(solver, M, B, X);
924 HYPRE_ParCSRGMRESGetNumIterations(solver, &iteration);
925 HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(solver, &residue);
927 HYPRE_Int converged = 0;
928 HYPRE_GMRESGetConverged(solver, &converged);
929 error |= (!converged);
931 HYPRE_ParCSRGMRESDestroy(solver);
937 HYPRE_IJMatrix m_hypre_ijmatrix =
nullptr;
938 HYPRE_ParCSRMatrix m_hypre_parmatrix =
nullptr;