196class AlephMatrixSloop:
public IAlephMatrix{
204 Integer index):IAlephMatrix(tm,kernel,index),
205 m_sloop_matrix(NULL){
206 debug()<<
"\t\t[AlephMatrixSloop::AlephMatrixSloop] NEW AlephMatrixSloop";
213 debug() <<
"\33[1;5;32m\t\t\t[~AlephMatrixSloop]\33[0m";
214 delete m_sloop_matrix;
222 void AlephMatrixCreate(
void)
224 ItacFunction(AlephMatrixSloop);
225 debug()<<
"\t\t[AlephMatrixSloop::AlephMatrixCreate] create new SLOOP::SLOOPDistMatrix";
228 m_sloop_matrix =
new SLOOP::SLOOPDistMatrix(*aleph_topology_sloop->m_sloop_topology,
229 *aleph_topology_sloop->m_sloop_msg,
233 if (!m_sloop_matrix)
throw Exception(
"AlephSolverMatrix::create",
"new SLOOPDistMatrix() failed");
236 if (aleph_topology_sloop->m_sloop_topology->get_type()==SLOOP::contiguous)
237 m_sloop_matrix->set_renumbering_opt(SLOOP::SLOOPDistMatrix::interface,
false);
238 m_sloop_matrix->set_renumbering_opt(SLOOP::SLOOPDistMatrix::processor,
false);
239 m_sloop_matrix->set_renumbering_opt(SLOOP::SLOOPDistMatrix::interior,
false);
240 m_sloop_matrix->init_length(m_kernel->topology()->gathered_nb_row_elements().unguardedBasePointer());
241 debug()<<
"\t\t[AlephMatrixSloop::AlephMatrixCreate] done";
248 void AlephMatrixSetFilled(
bool toggle){
249 ItacFunction(AlephMatrixSloop);
250 debug()<<
"\t\t[AlephMatrixSloop::AlephMatrixSetFilled]";
251 m_sloop_matrix->setfilled(toggle);
258 int AlephMatrixAssemble(
void){
259 ItacFunction(AlephMatrixSloop);
260 debug()<<
"\t\t[AlephMatrixSloop::AlephMatrixAssemble]";
261 m_sloop_matrix->configure();
269 void AlephMatrixFill(
int size,
int *rows,
int *cols,
double *values){
270 ItacFunction(AlephMatrixSloop);
271 debug()<<
"\t\t[AlephMatrixSloop::AlephMatrixFill] size="<<size;
272 m_sloop_matrix->locfill(values, rows, cols, size);
273 debug()<<
"\t\t[AlephMatrixSloop::AlephMatrixFill] done";
280 bool isAlreadySolved(SLOOP::SLOOPDistVector* x,
281 SLOOP::SLOOPDistVector* b,
282 SLOOP::SLOOPDistVector* tmp,
285 const Real res0 = b->norm_max();
286 const Real considered_as_null = params->minRHSNorm();
287 const bool convergence_analyse = params->convergenceAnalyse();
289 if (convergence_analyse)
290 debug() <<
"analyse convergence : norme max du second membre res0 : " << res0;
292 if (res0 < considered_as_null) {
294 residual_norm[0]= res0;
295 if (convergence_analyse)
296 debug() <<
"analyse convergence : le second membre du système linéaire est inférieur à "
297 << considered_as_null;
301 if (params->xoUser()) {
303 m_sloop_matrix->vector_multiply(*tmp,*x);
304 tmp->substract(*tmp,*b);
305 const Real residu= tmp->norm_max();
306 debug() <<
"[IAlephSloop::isAlreadySolved] residu="<<residu;
308 if (residu < considered_as_null) {
309 if (convergence_analyse) {
310 debug() <<
"analyse convergence : |Ax0-b| est inférieur à " << considered_as_null;
311 debug() <<
"analyse convergence : x0 est déjà solution du système.";
313 residual_norm[0] = residu;
316 const Real relative_error = residu / res0;
317 if (convergence_analyse)
318 debug() <<
"analyse convergence : résidu initial : " << residu
319 <<
" --- residu relatif initial (residu/res0) : " << residu / res0;
321 if (relative_error < (params->epsilon())) {
322 if (convergence_analyse)
323 debug() <<
"analyse convergence : X est déjà solution du système";
324 residual_norm[0] = residu;
343 ItacFunction(AlephMatrixSloop);
345 debug() <<
"\t\t[AlephMatrixSloop::SloopSolve] getTopologyImplementation #"<<m_index;
356 SLOOP::SLOOPDistVector* solution = x_sloop->m_sloop_vector;
357 SLOOP::SLOOPDistVector* RHS = b_sloop->m_sloop_vector;
358 SLOOP::SLOOPDistVector* temp = tmp_sloop->m_sloop_vector;
360 if (isAlreadySolved(solution,RHS,temp,residual_norm,solver_param)){
361 debug() <<
"\t[AlephMatrixSloop::AlephMatrixSolve] isAlreadySolved !";
367 sc=createSloopStopCriteria(solver_param, *sloopParallelInfo->m_sloop_msg);
370 global_solver=createSloopSolver(solver_param, *sloopParallelInfo->m_sloop_msg);
373 precond=createSloopPreconditionner(solver_param, *sloopParallelInfo->m_sloop_msg);
375 this->setSloopSolverParameters(solver_param, global_solver.
get());
376 this->setSloopPreconditionnerParameters(solver_param, precond.
get());
379 const bool normalize = normalizeSolverMatrix(solver_param);
381 debug() <<
"\t\t[AlephMatrixSloop::SloopSolve] Normalize";
382 diag =
new SLOOP::SLOOPDistVector(*sloopParallelInfo->m_sloop_topology, *sloopParallelInfo->m_sloop_msg);
383 global_solver->normalize(*m_sloop_matrix, *diag, *solution, *RHS);
386 const bool xo = solver_param->xoUser();
387 switch (solver_param->precond()) {
388 case TypesSolver::NONE:
391 debug() <<
"\t\t[AlephMatrixSloop::SloopSolve] xo à true (sans preconditionnement)";
392 status = global_solver->solve(*m_sloop_matrix, *solution, *RHS, *sc);
394 debug() <<
"\t\t[AlephMatrixSloop::SloopSolve] xo à false (sans preconditionnement)";
395 status = global_solver->solve_b(*m_sloop_matrix, *solution, *RHS, *sc);
401 debug() <<
"\t\t[AlephMatrixSloop::SloopSolve] xo à true (avec preconditionnement)";
402 status = global_solver->solve(*m_sloop_matrix, *solution, *RHS, *precond, *sc);
404 debug() <<
"\t\t[AlephMatrixSloop::SloopSolve] xo à false (avec preconditionnement)";
405 status = global_solver->solve_b(*m_sloop_matrix, *solution, *RHS, *precond, *sc);
411 debug() <<
"\t\t[AlephMatrixSloop::SloopSolve] INV-normalize";
412 global_solver->inv_normalize(*m_sloop_matrix, *diag, *solution, *RHS);
415 nb_iteration = global_solver->get_iteration();
416 residual_norm[0] = sc->get_criteria();
417 Integer max_iteration= global_solver->get_max_iteration();
418 residual_norm[3] = global_solver->get_stagnation();
420 if ((solver_param->getCriteriaStop()==TypesSolver::STAG)||(solver_param->getCriteriaStop()==TypesSolver::NIter)){
424 if (nb_iteration == max_iteration && solver_param->stopErrorStrategy()){
425 info() <<
"\n============================================================";
426 info() <<
"\nCette erreur est retournée après " << nb_iteration <<
"\n";
427 info() <<
"\nOn a atteind le nombre max d'itérations du solveur ";
428 info() <<
"\nIl est possible de demander au code de ne pas tenir compte de cette erreur.";
429 info() <<
"\nVoir la documentation du jeu de données concernant le service solveur.";
430 info() <<
"\n============================================================";
431 throw Exception(
"AlephMatrixSloop::SloopSolve",
"On a atteind le nombre max d'itérations du solveur");
434 debug() <<
"\t\t[AlephMatrixSloop::SloopSolve] nbIteration="<< global_solver->get_iteration()
435 <<
", criteria=" << residual_norm[0] <<
", stagnation=" << residual_norm[3];
443 void writeToFile(
const String file_name){
444 ItacFunction(AlephMatrixSloop);
445 m_sloop_matrix->write_to_file(file_name.
localstr());
452 SLOOP::SLOOPSolver* createSloopSolver(
AlephParams* solver_param, SLOOP::SLOOPMsg& info_sloop_msg){
453 const Integer output_level = solver_param->getOutputLevel();
454 ItacFunction(AlephMatrixSloop);
455 switch (solver_param->method()) {
456 case TypesSolver::PCG:
return new SLOOP::SLOOPCGSolver(info_sloop_msg, output_level);
457 case TypesSolver::BiCGStab:
return new SLOOP::SLOOPBiCGStabSolver(info_sloop_msg, output_level);
458 case TypesSolver::BiCGStab2:
return new SLOOP::SLOOPBiCGStabSolver(info_sloop_msg, output_level);
459 case TypesSolver::GMRES:
return new SLOOP::SLOOPGMRESSolver(info_sloop_msg, output_level);
460 case TypesSolver::SAMG:
return new SLOOP::SLOOPAMGSolver(info_sloop_msg, output_level);
461 case TypesSolver::QMR:
return new SLOOP::SLOOPQMRSolver(info_sloop_msg, output_level);
462 case TypesSolver::SuperLU:{
463 SLOOP::SLOOPSolver *solver =
new SLOOP::SLOOPSuperLUSolver(info_sloop_msg, output_level);
464 solver->set_parameter(SLOOP::sv_residual_calculation, 1);
467 default:
throw Exception(
"AlephMatrixSloop::createSloopSolver",
"Type de solver non accessible pour la bibliothèque SLOOP");
476 SLOOP::SLOOPPreconditioner* createSloopPreconditionner(
AlephParams* solver_param,
477 SLOOP::SLOOPMsg& info_sloop_msg){
478 const Integer output_level = solver_param->getOutputLevel();
479 ItacFunction(AlephMatrixSloop);
480 switch (solver_param->precond()) {
481 case TypesSolver::AINV:
484 if (TypesSolver::PCG == solver_param->method())
485 return new SLOOP::SLOOPAINVPC(info_sloop_msg, output_level);
487 return new SLOOP::SLOOPMAINVPC(info_sloop_msg, output_level);
488 case TypesSolver::DIAGONAL:
return new SLOOP::SLOOPDiagPC(info_sloop_msg, output_level);
489 case TypesSolver::AMG:
return new SLOOP::SLOOPAMGPC(info_sloop_msg, output_level);
490 case TypesSolver::IC:
return new SLOOP::SLOOPCholPC(info_sloop_msg, output_level);
491 case TypesSolver::POLY:
return new SLOOP::SLOOPPolyPC(info_sloop_msg, output_level);
492 case TypesSolver::ILU:
return new SLOOP::SLOOPILUPC(info_sloop_msg, output_level);
493 case TypesSolver::ILUp:
return new SLOOP::SLOOPILUPC(info_sloop_msg, output_level);
494 case TypesSolver::SPAIstat:
return new SLOOP::SLOOPSPAIPC(info_sloop_msg, output_level);
495 case TypesSolver::SPAIdyn:
return new SLOOP::SLOOPSPAIPC(info_sloop_msg, output_level);
496 case TypesSolver::NONE:
return NULL;
497 default:
throw Exception(
"AlephMatrixSloop::createSloopPreconditionner",
498 "préconditionneur non accessible pour la bibliothèque SLOOP");
507 SLOOP::SLOOPStopCriteria * createSloopStopCriteria(
AlephParams* solver_param,
508 SLOOP::SLOOPMsg& info_sloop_msg) {
509 ItacFunction(AlephMatrixSloop);
510 switch(solver_param->getCriteriaStop()){
511 case TypesSolver::RR0:
return new SLOOP::SLOOPStopCriteriaRR0();
512 case TypesSolver::R:
return new SLOOP::SLOOPStopCriteriaR();
513 case TypesSolver::RCB:
return new SLOOP::SLOOPStopCriteriaRCB();
514 case TypesSolver::RBinf:
return new SLOOP::SLOOPStopCriteriaRBinf();
515 case TypesSolver::EpsA:
return new SLOOP::SLOOPStopCriteriaEpsA();
516 case TypesSolver::NIter:
return new SLOOP::SLOOPStopCriteriaNIter();
517 case TypesSolver::RR0inf:
return new SLOOP::SLOOPStopCriteriaRR0inf();
518 case TypesSolver::STAG:
return new SLOOP::SLOOPStopCriteriaSTAG();
519 case TypesSolver::RB:
520 default:
return new SLOOP::SLOOPStopCriteriaRB();
528 void setSloopSolverParameters(
AlephParams* solver_param,
529 SLOOP::SLOOPSolver* sloop_solver){
531 Integer gamma = solver_param->gamma();
532 double alpha = solver_param->alpha();
533 ItacFunction(AlephMatrixSloop);
534 error+=sloop_solver->set_parameter(SLOOP::sv_epsilon, solver_param->epsilon());
535 error+=sloop_solver->set_parameter(SLOOP::sv_max_iter, solver_param->maxIter());
536 switch (solver_param->method()) {
537 case TypesSolver::PCG:
break;
538 case TypesSolver::QMR:
break;
539 case TypesSolver::SuperLU:
break;
540 case TypesSolver::BiCGStab:
error += sloop_solver->set_parameter(SLOOP::bicg_dimension, 1);
break;
541 case TypesSolver::BiCGStab2:
error += sloop_solver->set_parameter(SLOOP::bicg_dimension, 2);
break;
542 case TypesSolver::GMRES:
543 error += sloop_solver->set_parameter(SLOOP::gmres_order, 20);
544 error += sloop_solver->set_parameter(SLOOP::gmres_type, SLOOP::ICGS);
546 case TypesSolver::SAMG: {
547 error += sloop_solver->set_parameter(SLOOP::amg_buffer_size, 200);
549 if (gamma == -1) gamma = 50;
550 error += sloop_solver->set_parameter(SLOOP::amg_level, gamma);
553 if (alpha < 0.0) alpha = 0.1;
554 error += sloop_solver->set_parameter(SLOOP::amg_alpha, alpha);
557 error += sloop_solver->set_parameter(SLOOP::amg_smoother_iter, solver_param->getAmgSmootherIter());
559 error += sloop_solver->set_parameter(SLOOP::amg_iter, solver_param->getAmgCycle());
561 SLOOP::SLOOPAMGCoarseningOption coarsening_option =
562 static_cast<SLOOP::SLOOPAMGCoarseningOption
> (solver_param->getAmgCoarseningOption());
563 error += sloop_solver->set_parameter(SLOOP::amg_coarsening_option, coarsening_option);
566 SLOOP::SLOOPAMGCoarseSolverOption coarse_solver_option =
567 static_cast<SLOOP::SLOOPAMGCoarseSolverOption
> (solver_param->getAmgCoarseSolverOption());
568 error += sloop_solver->set_parameter(SLOOP::amg_coarse_solver_option, coarse_solver_option);
571 SLOOP::SLOOPAMGSmootherOption smoother_option =
572 static_cast<SLOOP::SLOOPAMGSmootherOption
> (solver_param->getAmgSmootherOption());
573 error += sloop_solver->set_parameter(SLOOP::amg_smoother_option, smoother_option);
577 default:
throw Exception(
"AlephMatrixSloop::setSloopSolverParameters",
578 "Type de solver SLOOP non prévu dans la gestion des paramètres");
581 throw Exception(
"AlephMatrixSloop::setSloopSolverParameters",
"set_parameter() failed");
587 void setSloopPreconditionnerParameters(
AlephParams* solver_param,
588 SLOOP::SLOOPPreconditioner* preconditionner){
589 const TypesSolver::ePreconditionerMethod precond_method = solver_param->precond();
590 double alpha = solver_param->alpha();
591 int gamma = solver_param->gamma();
592 const String function_id =
"SolverMatrixSloop::setSloopPreconditionnerParameters";
593 ItacFunction(AlephMatrixSloop);
595 if (alpha < 0.0) alpha = 0.1;
596 switch (precond_method) {
597 case TypesSolver::NONE:
break;
598 case TypesSolver::DIAGONAL:
break;
599 case TypesSolver::AMG: {
600 if (gamma == -1) gamma = 50;
601 preconditionner->set_parameter(SLOOP::amg_level, gamma);
602 preconditionner->set_parameter(SLOOP::amg_buffer_size, 200);
603 preconditionner->set_parameter(SLOOP::amg_solver_iter, solver_param->getAmgSolverIter());
604 preconditionner->set_parameter(SLOOP::amg_iter, solver_param->getAmgCycle());
605 preconditionner->set_parameter(SLOOP::amg_smoother_iter, solver_param->getAmgSmootherIter());
607 SLOOP::SLOOPAMGSmootherOption smoother_option =
608 static_cast<SLOOP::SLOOPAMGSmootherOption
> (solver_param->getAmgSmootherOption());
610 SLOOP::SLOOPAMGCoarseningOption coarsening_option =
611 static_cast<SLOOP::SLOOPAMGCoarseningOption
> (solver_param->getAmgCoarseningOption());
613 SLOOP::SLOOPAMGCoarseSolverOption coarse_solver_option =
614 static_cast<SLOOP::SLOOPAMGCoarseSolverOption
> (solver_param->getAmgCoarseSolverOption());
616 if (TypesSolver::PCG == solver_param->method()) {
618 switch (solver_param->getAmgSmootherOption()) {
619 case TypesSolver::CG_smoother:
620 case TypesSolver::Rich_IC_smoother:
621 case TypesSolver::Rich_AINV_smoother:
622 case TypesSolver::SymHybGSJ_smoother:
623 case TypesSolver::Rich_IC_block_smoother:
624 case TypesSolver::SymHybGSJ_block_smoother:
626 default:
throw Exception(
"AlephMatrixSloop::setSloopPreconditionnerParameters",
627 "choix du smoother incorrect pour une matrice symetrique");
631 switch (solver_param->getAmgSmootherOption()) {
632 case TypesSolver::CG_smoother:
633 case TypesSolver::Rich_IC_smoother:
634 case TypesSolver::Rich_AINV_smoother:
635 case TypesSolver::Rich_IC_block_smoother:
636 case TypesSolver::SymHybGSJ_block_smoother:
637 throw Exception(
"AlephMatrixSloop::setSloopPreconditionnerParameters",
638 "choix du smoother incorrect avec une matrice non-symetrique ");
639 case TypesSolver::SymHybGSJ_smoother:
641 solver_param->setAmgSmootherOption(TypesSolver::HybGSJ_smoother);
646 switch (solver_param->getAmgCoarseSolverOption()) {
647 case TypesSolver::CG_coarse_solver:
648 case TypesSolver::Cholesky_coarse_solver:
649 solver_param->setAmgCoarseSolverOption(TypesSolver::LU_coarse_solver);
651 default:solver_param->setAmgCoarseSolverOption(TypesSolver::BiCGStab_coarse_solver);
656 preconditionner->set_parameter(SLOOP::amg_smoother_option, smoother_option);
658 preconditionner->set_parameter(SLOOP::amg_coarsening_option, coarsening_option);
660 preconditionner->set_parameter(SLOOP::amg_coarse_solver_option, coarse_solver_option);
664 case TypesSolver::POLY:
665 if (gamma == -1) gamma = 3;
667 case TypesSolver::AINV:
668 if (gamma == -1) gamma = 0;
671 case TypesSolver::IC:
672 case TypesSolver::ILU:
if (gamma == -1) gamma = 0;
break;
673 case TypesSolver::ILUp:
if (gamma == -1) gamma = 1;
break;
674 case TypesSolver::SPAIstat:
675 preconditionner->set_parameter(SLOOP::spai_sparsity,SLOOP::StatSparsity);
676 preconditionner->set_parameter(SLOOP::spai_init_sparsity,SLOOP::PowerSparsity);
677 preconditionner->set_parameter(SLOOP::spai_power_level, 1);
678 preconditionner->set_parameter(SLOOP::spai_Amax_row_size, 30);
679 preconditionner->set_parameter(SLOOP::spai_A_drop_eps, 0.001);
681 case TypesSolver::SPAIdyn:
682 preconditionner->set_parameter(SLOOP::spai_sparsity, SLOOP::DynSparsity);
683 preconditionner->set_parameter(SLOOP::spai_init_sparsity, SLOOP::DiagSparsity);
684 preconditionner->set_parameter(SLOOP::spai_Amax_row_size, 10);
685 preconditionner->set_parameter(SLOOP::spai_A_drop_eps, 0.001);
688 throw Exception(
"AlephMatrixSloop::setSloopPreconditionnerParameters",
"Préconditionneur inconnu.");
692 switch (precond_method) {
693 case TypesSolver::NONE:
694 case TypesSolver::SPAIstat:
695 case TypesSolver::SPAIdyn:
698 case TypesSolver::AMG:
699 preconditionner->set_parameter(SLOOP::amg_alpha, alpha);
700 preconditionner->set_parameter(SLOOP::amg_level, gamma);
701 preconditionner->set_parameter(SLOOP::amg_parallel_opt, 0);
703 preconditionner->set_parameter(SLOOP::amg_coarse_solver_sc_option,SLOOP::RB);
707 preconditionner->set_parameter(SLOOP::pc_alpha, alpha);
708 preconditionner->set_parameter(SLOOP::pc_nbelem, gamma);
709 preconditionner->set_parameter(SLOOP::pc_parallel_opt, 1);
710 preconditionner->set_parameter(SLOOP::pc_order, gamma);
719 bool normalizeSolverMatrix(
AlephParams* solver_param){
720 ItacFunction(AlephMatrixSloop);
721 switch (solver_param->precond()){
722 case TypesSolver::AINV:
724 case TypesSolver::SPAIstat:
725 case TypesSolver::SPAIdyn:
return true;
726 case TypesSolver::AMG:
727 case TypesSolver::NONE:
728 case TypesSolver::DIAGONAL:
729 case TypesSolver::IC:
730 case TypesSolver::POLY:
731 case TypesSolver::ILU:
732 case TypesSolver::ILUp:
return false;
733 default:
throw Exception(
"AlephMatrixSloop::normalizeSolverMatrix",
"Préconditionneur inconnu.");
738 SLOOP::SLOOPMatrix* m_sloop_matrix;