14#include "arcane/aleph/AlephArcane.h"
16#include "arcane/utils/StringList.h"
18#define MPI_COMM_SUB *(MPI_Comm*)(m_kernel->subParallelMng(m_index)->getMPICommunicator())
24#if PETSC_VERSION_GE(3,6,1)
26#elif PETSC_VERSION_(3,3,0)
28#elif PETSC_VERSION_(3,0,0)
31#error "Bad version of PETSc: should be 3.0.0, 3.3.0 or greater than 3.6.1."
45 static_assert(
sizeof(PetscInt)==
sizeof(
AlephInt),
"AlephInt and PetscInt should have the same size");
46 return reinterpret_cast<PetscInt*
>(v);
49const PetscInt* _toPetscInt(
const AlephInt* v)
51 static_assert(
sizeof(PetscInt)==
sizeof(
AlephInt),
"AlephInt and PetscInt should have the same size");
52 return reinterpret_cast<const PetscInt*
>(v);
56petscCheck(PetscErrorCode error_code)
61#if PETSC_VERSION_GE(3,19,0)
62 if (error_code==PETSC_SUCCESS)
64 PetscCallVoid(error_code);
65 ARCANE_FATAL(
"Error in PETSc backend errcode={0}",error_code);
68 ARCANE_FATAL(
"Error in PETSc backend errcode={0}",error_code);
77class AlephTopologyPETSc
78:
public IAlephTopology
83 ~AlephTopologyPETSc()
override
85 debug() <<
"\33[1;5;32m\t\t\t[~AlephTopologyPETSc]\33[0m";
90 void backupAndInitialize()
override {}
91 void restore()
override {}
95 void _checkInitPETSc();
101class AlephVectorPETSc
108 void AlephVectorCreate(
void)
override;
110 int AlephVectorAssemble(
void)
override;
112 void writeToFile(
const String)
override;
113 Real LinftyNorm(
void);
117 Vec m_petsc_vector =
nullptr;
120 PetscInt jLower = -1;
126class AlephMatrixPETSc
133 void AlephMatrixCreate(
void)
override;
134 void AlephMatrixSetFilled(
bool)
override;
135 int AlephMatrixAssemble(
void)
override;
142 void writeToFile(
const String)
override;
146 Mat m_petsc_matrix =
nullptr;
147 KSP m_ksp_solver =
nullptr;
155bool global_need_petsc_finalize =
false;
161class PETScAlephFactoryImpl
169 ~PETScAlephFactoryImpl()
override
171 for (
auto* v : m_IAlephVectors)
173 for (
auto* v : m_IAlephMatrixs)
175 for (
auto* v : m_IAlephTopologys)
177 if (global_need_petsc_finalize){
183 global_need_petsc_finalize =
false;
189 void initialize()
override {}
196 m_IAlephTopologys.add(new_topology);
204 m_IAlephVectors.add(new_vector);
213 m_IAlephMatrixs.add(new_matrix);
235 if (!m_participating_in_solver) {
236 debug() <<
"\33[1;32m\t[AlephTopologyPETSc] Not concerned with this solver, returning\33[0m";
239 debug() <<
"\33[1;32m\t\t[AlephTopologyPETSc] @" <<
this <<
"\33[0m";
240 if (!m_kernel->isParallel()) {
241 PETSC_COMM_WORLD = PETSC_COMM_SELF;
253void AlephTopologyPETSc::
257 PetscBool is_petsc_initialized = PETSC_FALSE;
258 PetscInitialized(&is_petsc_initialized);
259 if (is_petsc_initialized==PETSC_TRUE)
262 AlephKernelSolverInitializeArguments& init_args = m_kernel->solverInitializeArgs();
263 bool has_args = init_args.hasValues();
264 debug() << Trace::Color::cyan() <<
"[AlephTopologyPETSc] Initializing PETSc. UseArg=" << has_args;
267 const CommandLineArguments& args = init_args.commandLineArguments();
268 PetscInitialize(args.commandLineArgc(),args.commandLineArgv(),
nullptr,
nullptr);
271 PetscInitializeNoArguments();
273 global_need_petsc_finalize =
true;
290 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorPETSc] new SolverVectorPETSc";
299 VecDestroy(&m_petsc_vector);
308void AlephVectorPETSc::
309AlephVectorCreate(
void)
311 for(
int iCpu=0;iCpu<m_kernel->size();++iCpu){
312 if (m_kernel->rank()!=m_kernel->solverRanks(m_index)[iCpu])
315 jLower=m_kernel->topology()->gathered_nb_row(iCpu);
316 jUpper=m_kernel->topology()->gathered_nb_row(iCpu+1);
320 VecCreateMPI(MPI_COMM_SUB,
322 m_kernel->topology()->nb_row_size(),
324 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorCreate] PETSC VectorCreate"
325 <<
" of local size=" <<jSize<<
"/"<<m_kernel->topology()->nb_row_size()
332void AlephVectorPETSc::
333AlephVectorSet(
const double *bfr_val,
const AlephInt *bfr_idx,
Integer size)
335 VecSetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val,INSERT_VALUES);
336 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorSet] "<<size<<
" values inserted!";
338 AlephVectorAssemble();
345int AlephVectorPETSc::
346AlephVectorAssemble(
void)
348 VecAssemblyBegin(m_petsc_vector);
349 VecAssemblyEnd(m_petsc_vector);
350 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorAssemble]";
358void AlephVectorPETSc::
361 VecGetValues(m_petsc_vector,size,_toPetscInt(bfr_idx),bfr_val);
362 debug()<<
"\t\t[AlephVectorPETSc::AlephVectorGet] fetched "<< size <<
" values!";
369void AlephVectorPETSc::
370writeToFile(
const String filename)
372 ARCANE_UNUSED(filename);
373 debug()<<
"\t\t[AlephVectorPETSc::writeToFile]";
379Real AlephVectorPETSc::
383 VecNorm(m_petsc_vector,NORM_INFINITY,&val);
398 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixPETSc] new AlephMatrixPETSc";
408 MatDestroy(&m_petsc_matrix);
410 KSPDestroy(&m_ksp_solver);
419void AlephMatrixPETSc::
420AlephMatrixCreate(
void)
425 for(
int iCpu=0;iCpu<m_kernel->size();++iCpu ){
426 if (m_kernel->rank()!=m_kernel->solverRanks(m_index)[iCpu])
continue;
427 if (ilower==-1) ilower=m_kernel->topology()->gathered_nb_row(iCpu);
428 iupper=m_kernel->topology()->gathered_nb_row(iCpu+1);
433 PetscInt* row_data = _toPetscInt(m_kernel->topology()->gathered_nb_row_elements().subView(ilower,size).data());
437 PetscInt M = m_kernel->topology()->nb_row_size();
440 MatCreate(PETSC_COMM_WORLD, &m_petsc_matrix);
441 MatSetSizes(m_petsc_matrix, m, n, M, N);
444 MatSetType(m_petsc_matrix, MATAIJ);
445 MatMPIAIJSetPreallocation(m_petsc_matrix, 0, row_data, 0, row_data);
448 MatSetFromOptions(m_petsc_matrix);
449#if PETSC_VERSION_LE(3,0,0)
453 PETSC_VERSION_MatCreate(MPI_COMM_SUB,
456 m_kernel->topology()->nb_row_size(),
457 m_kernel->topology()->nb_row_size(),
467 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE);
468 MatSetUp(m_petsc_matrix);
469 debug()<<
"\t\t[AlephMatrixPetsc::AlephMatrixCreate] PETSC MatrixCreate idx:"<<m_index
470 <<
", ("<<ilower<<
"->"<<(iupper-1)<<
")";
477void AlephMatrixPETSc::
478AlephMatrixSetFilled(
bool)
480 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixSetFilled] done";
481 MatSetOption(m_petsc_matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
488int AlephMatrixPETSc::
489AlephMatrixAssemble(
void)
491 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixAssemble] AlephMatrixAssemble";
492 MatAssemblyBegin(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
493 MatAssemblyEnd(m_petsc_matrix,MAT_FINAL_ASSEMBLY);
501void AlephMatrixPETSc::
504 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixFill] size="<<size;
505 for(
int i=0;i<size;i++){
507 MatSetValue(m_petsc_matrix, rows[i], cols[i], values[i], INSERT_VALUES);
509 debug()<<
"\t\t[AlephMatrixPETSc::AlephMatrixFill] done";
511 AlephMatrixAssemble();
519Real AlephMatrixPETSc::
524 throw FatalErrorException(
"LinftyNormVectorProductAndSub",
"error");
531bool AlephMatrixPETSc::
541 ARCANE_UNUSED(residual_norm);
542 ARCANE_UNUSED(params);
550int AlephMatrixPETSc::
560 KSPConvergedReason reason;
562 AlephVectorPETSc* x_petsc =
dynamic_cast<AlephVectorPETSc*
> (x->implementation());
563 AlephVectorPETSc* b_petsc =
dynamic_cast<AlephVectorPETSc*
> (b->implementation());
567 const bool is_parallel = m_kernel->subParallelMng(m_index)->isParallel();
569 Vec solution = x_petsc->m_petsc_vector;
570 Vec RHS = b_petsc->m_petsc_vector;
572 debug()<<
"[AlephMatrixSolve]";
574 KSPDestroy(&m_ksp_solver);
575 KSPCreate(MPI_COMM_SUB,&m_ksp_solver);
576#if PETSC_VERSION_GE(3,6,1)
577 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix);
579 KSPSetOperators(m_ksp_solver,m_petsc_matrix,m_petsc_matrix,SAME_NONZERO_PATTERN);
582 switch(solver_param->method()){
584 case TypesSolver::PCG : KSPSetType(m_ksp_solver,KSPCG);
break;
586 case TypesSolver::BiCGStab : KSPSetType(m_ksp_solver,KSPBCGS);
break;
589 case TypesSolver::BiCGStab2: KSPSetType(m_ksp_solver,KSPIBCGS);
break;
591 case TypesSolver::GMRES : KSPSetType(m_ksp_solver,KSPGMRES);
break;
592 default :
throw ArgumentException(
"AlephMatrixPETSc::AlephMatrixSolve",
"Unknown solver method");
595 KSPGetPC(m_ksp_solver,&prec);
597 switch(solver_param->precond()){
598 case TypesSolver::NONE : PCSetType(prec,PCNONE);
break;
600 case TypesSolver::DIAGONAL : PCSetType(prec,PCJACOBI);
break;
602 case TypesSolver::ILU:
606 PCSetType(prec,PCBJACOBI);
608 PCSetType(prec,PCILU);
611 case TypesSolver::IC:
614 PCSetType(prec,PCBJACOBI);
615 PetscInt nlocal, first;
618 KSPSetUp(m_ksp_solver);
619 PCBJacobiGetSubKSP(prec, &nlocal, &first, &subksp);
620 for (
int i = 0; i < nlocal; i++) {
621 KSPGetPC(subksp[i], &subpc);
622 PCSetType(subpc, PCICC);
626 PCSetType(prec,PCICC);
629 case TypesSolver::SPAIstat : PCSetType(prec,PCSPAI);
break;
633 case TypesSolver::AMG:
637 KSPSetType(m_ksp_solver,KSPFGMRES);
638 PCSetType(prec,PCMG);
640 PCMGSetLevels(prec,1,
641 (MPI_Comm*)(m_kernel->subParallelMng(m_index)->getMPICommunicator()));
643 PCMGSetType(prec,PC_MG_MULTIPLICATIVE);
645 PCMGSetCycleType(prec,PC_MG_CYCLE_V);
647#if PETSC_VERSION_GE(3,10,2)
648 PCMGSetNumberSmooth(prec, 1);
650 PCMGSetNumberSmoothDown(prec, 1);
651 PCMGSetNumberSmoothUp(prec, 1);
662 case TypesSolver::AINV :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement AINV indisponible");
663 case TypesSolver::SPAIdyn :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement SPAIdyn indisponible");
664 case TypesSolver::ILUp :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement ILUp indisponible");
665 case TypesSolver::POLY :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement POLY indisponible");
666 default :
throw ArgumentException(
"AlephMatrixPETSc",
"preconditionnement inconnu");
669 if (solver_param->xoUser())
670 KSPSetInitialGuessNonzero(m_ksp_solver,PETSC_TRUE);
672 KSPSetTolerances(m_ksp_solver,
673 solver_param->epsilon(),
676 solver_param->maxIter());
679 KSPSetFromOptions(m_ksp_solver);
682 PCSetFromOptions(prec) ;
685 KSPSetUp(m_ksp_solver);
690 KSPGetType(m_ksp_solver,&kt);
693 info(4) <<
"[AlephMatrixSolve] Will use KSP type :" << kt <<
" PC type : "<< pt
694 <<
" is_parallel=" << is_parallel;
697 debug()<<
"[AlephMatrixSolve] All set up, now solving";
698 petscCheck(KSPSolve(m_ksp_solver,RHS,solution));
699 debug()<<
"[AlephMatrixSolve] solved";
700 petscCheck(KSPGetConvergedReason(m_ksp_solver,&reason));
703#if !PETSC_VERSION_(3,0,0)
704 case(KSP_CONVERGED_RTOL_NORMAL):{
break;}
705 case(KSP_CONVERGED_ATOL_NORMAL):{
break;}
707 case(KSP_CONVERGED_RTOL):{
break;}
708 case(KSP_CONVERGED_ATOL):{
break;}
709 case(KSP_CONVERGED_ITS):{
break;}
710#if !PETSC_VERSION_GE(3,19,0)
711 case(KSP_CONVERGED_CG_NEG_CURVE):{
break;}
712 case(KSP_CONVERGED_CG_CONSTRAINED):{
break;}
714 case(KSP_CONVERGED_STEP_LENGTH):{
break;}
715 case(KSP_CONVERGED_HAPPY_BREAKDOWN):{
break;}
717 case(KSP_DIVERGED_NULL):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NULL");}
718 case(KSP_DIVERGED_ITS):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_ITS");}
719 case(KSP_DIVERGED_DTOL):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_DTOL");}
720 case(KSP_DIVERGED_BREAKDOWN):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_BREAKDOWN");}
721 case(KSP_DIVERGED_BREAKDOWN_BICG):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_BREAKDOWN_BICG");}
722 case(KSP_DIVERGED_NONSYMMETRIC):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NONSYMMETRIC");}
723 case(KSP_DIVERGED_INDEFINITE_PC):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_INDEFINITE_PC");}
724#if PETSC_VERSION_GE(3,6,1)
725 case(KSP_DIVERGED_NANORINF):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NANORINF");}
727 case(KSP_DIVERGED_NAN):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_NAN");}
729 case(KSP_DIVERGED_INDEFINITE_MAT):{
throw Exception(
"AlephMatrixPETSc::Solve",
"KSP_DIVERGED_INDEFINITE_MAT");}
731 case(KSP_CONVERGED_ITERATING):{
break;}
732 default:
throw Exception(
"AlephMatrixPETSc::Solve",
"");
735 KSPGetIterationNumber(m_ksp_solver,&its);
738 KSPGetResidualNorm(m_ksp_solver,&norm);
749void AlephMatrixPETSc::
750writeToFile(
const String filename)
752 ARCANE_UNUSED(filename);
#define ARCANE_CHECK_POINTER(ptr)
Macro retournant le pointeur ptr s'il est non nul ou lancant une exception s'il est nul.
#define ARCANE_FATAL(...)
Macro envoyant une exception FatalErrorException.
#define ARCANE_REGISTER_APPLICATION_FACTORY(aclass, ainterface, aname)
Enregistre un service de fabrique pour la classe aclass.
AbstractService(const ServiceBuildInfo &)
Constructeur à partir d'un ServiceBuildInfo.
Paramètres d'un système linéraire.
Vecteur d'un système linéaire.
Interface d'une fabrique d'implémentation pour Aleph.
virtual void * getMPICommunicator()=0
Adresse du communicateur MPI associé à ce gestionnaire.
Interface du gestionnaire de traces.
Structure contenant les informations pour créer un service.
Chaîne de caractères unicode.
TraceMessageDbg debug(Trace::eDebugLevel=Trace::Medium) const
Flot pour un message de debug.
TraceMessage info() const
Flot pour un message d'information.
Vecteur 1D de données avec sémantique par valeur (style STL).
-*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
Int32 Integer
Type représentant un entier.
double Real
Type représentant un réel.
int AlephInt
Type par défaut pour indexer les lignes et les colonnes des matrices et vecteurs.