Arcane  v3.16.2.0
Documentation développeur
Chargement...
Recherche...
Aucune correspondance
AlephSloop.cc
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2023 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7/*---------------------------------------------------------------------------*/
8/* AlephSloop.cc (C) 2010-2023 */
9/* */
10/*---------------------------------------------------------------------------*/
11/*---------------------------------------------------------------------------*/
12
13#define HASMPI
14#define SLOOP_MATH
15#define OMPI_SKIP_MPICXX
16#define MPICH_SKIP_MPICXX
17#define PARALLEL_SLOOP
18#include <mpi.h>
19#include "SLOOP.h"
20
21#ifndef ItacFunction
22#define ItacFunction(x)
23#endif
24
25#include "arcane/utils/FatalErrorException.h"
26#include "arcane/aleph/AlephArcane.h"
27
28/*---------------------------------------------------------------------------*/
29/*---------------------------------------------------------------------------*/
30
31namespace Arcane
32{
33
34/*---------------------------------------------------------------------------*/
35/*---------------------------------------------------------------------------*/
36
37
38// ****************************************************************************
39// * Class AlephTopologySloop
40// ****************************************************************************
41class AlephTopologySloop: public IAlephTopology{
42public:
43 AlephTopologySloop(ITraceMng* tm,
44 AlephKernel *kernel,
45 Integer index,
46 Integer nb_row_size):
47 IAlephTopology(tm, kernel, index, nb_row_size),
48 m_sloop_comminfo(NULL),
49 m_world_comminfo(NULL),
50 m_sloop_msg(NULL),
51 m_sloop_topology(NULL){
52 ItacFunction(AlephTopologySloop);
53
54 if (!m_kernel->isParallel()){
55 m_world_comminfo=new SLOOP::SLOOPSeqCommInfo();
56 }else{
57 m_world_comminfo=new SLOOP::SLOOPMPICommInfo(MPI_COMM_WORLD);
58 }
59
60 if (!m_participating_in_solver){
61 debug() << "\33[1;32m\t[AlephTopologySloop::AlephParallelInfoSloop] Not concerned with this solver, returning\33[0m";
62 return;
63 }
64
65 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] @"<<this<<"\33[0m";
66
67 if (!m_kernel->isParallel()){
68 m_sloop_comminfo=new SLOOP::SLOOPSeqCommInfo();
69 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] SEQCommInfo @"<<m_sloop_comminfo<<"\33[0m";
70 }else{
71 m_sloop_comminfo=new SLOOP::SLOOPMPICommInfo(*(SLOOP::SLOOP_Comm*)
72 (m_kernel->subParallelMng(index)->getMPICommunicator()));
73 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] MPICommInfo @"<<m_sloop_comminfo<<"\33[0m";
74 }
75
76 m_sloop_msg=new SLOOP::SLOOPMsg(m_sloop_comminfo);
77 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] SLOOPMsg @"<<m_sloop_msg<<"\33[0m";
78
79 m_sloop_topology=new SLOOP::SLOOPTopology(nb_row_size,
80 m_kernel->topology()->rowLocalRange(index),
81 *m_sloop_comminfo, *m_sloop_msg);
82 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] SLOOPTopology @"<<m_sloop_topology<<"\33[0m";
83 }
84 ~AlephTopologySloop(){
85 debug() << "\33[1;5;32m\t\t\t[~AlephTopologySloop] deleting m_sloop_msg, m_sloop_comminfo & m_sloop_topology\33[0m";
86 delete m_sloop_msg;
87 delete m_sloop_comminfo;
88 delete m_sloop_topology;
89 }
90public:
91 // On backup la session courante et on initialise la nouvelle session avec notre m_sloop_comminfo
92 void backupAndInitialize(){
93 if (!m_participating_in_solver) return;
94 SLOOP::SLOOPInitSession(m_sloop_comminfo,SLOOP::WITHOUT_EXCEPTIONS,SLOOP::OUTPUT_UNIQ);
95 }
96 void restore(){
97 // Dans tous les cas, tout le monde restore la WORLD session
98 SLOOP::SLOOPInitSession(m_world_comminfo,SLOOP::WITHOUT_EXCEPTIONS,SLOOP::OUTPUT_UNIQ);
99
100 }
101public:
102 SLOOP::SLOOPCommInfo *m_sloop_comminfo;
103 SLOOP::SLOOPCommInfo *m_world_comminfo;
104 SLOOP::SLOOPMsg *m_sloop_msg;
105 SLOOP::SLOOPTopology *m_sloop_topology;
106};
107
108
109// ****************************************************************************
110// * AlephVectorSloop
111// ****************************************************************************
112class AlephVectorSloop
113: public IAlephVector
114{
115 public:
116 AlephVectorSloop(ITraceMng* tm,
117 AlephKernel *kernel,
118 Integer index):IAlephVector(tm,kernel,index),
119 m_sloop_vector(NULL){
120 debug()<<"\t\t[AlephVectorSloop::AlephVectorSloop] NEW AlephVectorSloop";
121 }
122
123 // ****************************************************************************
124 // * ~AlephVectorSloop
125 // ****************************************************************************
126 ~AlephVectorSloop(){
127 debug() << "\33[1;5;32m\t\t\t[~AlephVectorSloop]\33[0m";
128 delete m_sloop_vector;
129 }
130
131 /******************************************************************************
132 * AlephVectorCreate
133 *****************************************************************************/
134 void AlephVectorCreate(void)
135 {
136 ItacFunction(AlephVectorSloop);
137 debug()<<"\t[AlephVectorSloop::AlephVectorCreate] new SLOOP::SLOOPDistVector";
138 AlephTopologySloop* aleph_topology_sloop = dynamic_cast<AlephTopologySloop*>(m_kernel->getTopologyImplementation(m_index));
139 ARCANE_CHECK_POINTER(aleph_topology_sloop);
140 m_sloop_vector = new SLOOP::SLOOPDistVector(*aleph_topology_sloop->m_sloop_topology,
141 *aleph_topology_sloop->m_sloop_msg);
142 if (!m_sloop_vector)
143 throw FatalErrorException(A_FUNCINFO, " new SLOOPDistVector failed");
144 debug()<<"\t[AlephVectorSloop::AlephVectorCreate] done";
145 }
146
147 // ****************************************************************************
148 // * AlephVectorSet
149 // ****************************************************************************
150 void AlephVectorSet(const double *bfr_val, const int *bfr_idx, Integer size)
151 {
152 debug()<<"\t[AlephVectorSloop::AlephVectorSet]";
153 if (m_sloop_vector->locfill(bfr_val, bfr_idx, size))
154 throw FatalErrorException(A_FUNCINFO, "locfill() failed");
155 }
156
157/******************************************************************************
158 * AlephVectorAssemble
159 *****************************************************************************/
160int AlephVectorAssemble(void){
161 ItacFunction(AlephVectorSloop);
162 debug()<<"\t\t[AlephVectorSloop::AlephVectorAssemble]";
163 return 0;
164}
165
166
167/******************************************************************************
168 * AlephVectorGet
169 *****************************************************************************/
170void AlephVectorGet(double *bfr_val, const int *bfr_idx, Integer size){
171 // Je récupère déjà les valeurs depuis le vecteur Sloop
172 debug()<<"\t[AlephVectorSloop::AlephVectorGet]";
173 if (m_sloop_vector->get_locval(bfr_val, bfr_idx, size))
174 throw Exception("AlephVectorSloop::AlephVectorGet", "get_locval() failed");
175}
176
177
178/******************************************************************************
179 * writeToFile
180 *****************************************************************************/
181 void writeToFile(const String file_name){
182 ItacFunction(AlephVectorSloop);
183 m_sloop_vector->write_to_file(file_name.localstr());
184 }
185
186public:
187 SLOOP::SLOOPDistVector* m_sloop_vector;
188};
189
190
191
192
193/******************************************************************************
194 AlephMatrixSloop
195*****************************************************************************/
196class AlephMatrixSloop: public IAlephMatrix{
197public:
198
199/******************************************************************************
200 * AlephMatrixSloop
201 *****************************************************************************/
202 AlephMatrixSloop(ITraceMng* tm,
203 AlephKernel *kernel,
204 Integer index):IAlephMatrix(tm,kernel,index),
205 m_sloop_matrix(NULL){
206 debug()<<"\t\t[AlephMatrixSloop::AlephMatrixSloop] NEW AlephMatrixSloop";
207 }
208
209// ****************************************************************************
210// * AlephMatrixSloop
211// ****************************************************************************
212 ~AlephMatrixSloop(){
213 debug() << "\33[1;5;32m\t\t\t[~AlephMatrixSloop]\33[0m";
214 delete m_sloop_matrix;
215 }
216
217
218
219 /******************************************************************************
220 * AlephMatrixCreate
221 *****************************************************************************/
222 void AlephMatrixCreate(void)
223 {
224 ItacFunction(AlephMatrixSloop);
225 debug()<<"\t\t[AlephMatrixSloop::AlephMatrixCreate] create new SLOOP::SLOOPDistMatrix";
226 AlephTopologySloop* aleph_topology_sloop = dynamic_cast<AlephTopologySloop*>(m_kernel->getTopologyImplementation(m_index));
227 ARCANE_CHECK_POINTER(aleph_topology_sloop);
228 m_sloop_matrix = new SLOOP::SLOOPDistMatrix(*aleph_topology_sloop->m_sloop_topology,
229 *aleph_topology_sloop->m_sloop_msg,
230 true,
231 false);
232
233 if (!m_sloop_matrix) throw Exception("AlephSolverMatrix::create","new SLOOPDistMatrix() failed");
234
235 // la renumerotation doit etre faite avant le remplissage de la matrice pour etre prise en compte
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";
242 }
243
244
245/******************************************************************************
246 * AlephMatrixSetFilled
247 *****************************************************************************/
248 void AlephMatrixSetFilled(bool toggle){
249 ItacFunction(AlephMatrixSloop);
250 debug()<<"\t\t[AlephMatrixSloop::AlephMatrixSetFilled]";
251 m_sloop_matrix->setfilled(toggle);
252 }
253
254
255/******************************************************************************
256 * AlephMatrixAssemble
257 *****************************************************************************/
258 int AlephMatrixAssemble(void){
259 ItacFunction(AlephMatrixSloop);
260 debug()<<"\t\t[AlephMatrixSloop::AlephMatrixAssemble]";
261 m_sloop_matrix->configure();
262 return 0;
263 }
264
265
266/******************************************************************************
267 * AlephMatrixFill
268 *****************************************************************************/
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";
274 }
275
276
277/******************************************************************************
278 * isAlreadySolved
279 *****************************************************************************/
280 bool isAlreadySolved(SLOOP::SLOOPDistVector* x,
281 SLOOP::SLOOPDistVector* b,
282 SLOOP::SLOOPDistVector* tmp,
283 Real* residual_norm,
284 AlephParams* params) {
285 const Real res0 = b->norm_max();
286 const Real considered_as_null = params->minRHSNorm();
287 const bool convergence_analyse = params->convergenceAnalyse();
288
289 if (convergence_analyse)
290 debug() << "analyse convergence : norme max du second membre res0 : " << res0;
291
292 if (res0 < considered_as_null) {
293 x->fill(Real(0.0));
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;
298 return true;
299 }
300
301 if (params->xoUser()) {
302 // on test si b est déjà solution du système à epsilon près
303 m_sloop_matrix->vector_multiply(*tmp,*x); // tmp=A*x
304 tmp->substract(*tmp,*b); // tmp=A*x-b
305 const Real residu= tmp->norm_max();
306 debug() << "[IAlephSloop::isAlreadySolved] residu="<<residu;
307
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.";
312 }
313 residual_norm[0] = residu;
314 return true;
315 }
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;
320
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;
325 return true;
326 }
327 }
328 return false;
329 }
330
331
332
333
334 /******************************************************************************
335 *****************************************************************************/
336 int AlephMatrixSolve(AlephVector* x,
337 AlephVector* b,
338 AlephVector* tmp,
339 Integer& nb_iteration,
340 Real* residual_norm,
341 AlephParams* solver_param)
342 {
343 ItacFunction(AlephMatrixSloop);
344 Integer status = 0;
345 debug() << "\t\t[AlephMatrixSloop::SloopSolve] getTopologyImplementation #"<<m_index;
346 AlephTopologySloop* sloopParallelInfo = dynamic_cast<AlephTopologySloop*>(m_kernel->getTopologyImplementation(m_index));
347
348 AlephVectorSloop* x_sloop = dynamic_cast<AlephVectorSloop*> (x->implementation());
349 AlephVectorSloop* b_sloop = dynamic_cast<AlephVectorSloop*> (b->implementation());
350 AlephVectorSloop* tmp_sloop = dynamic_cast<AlephVectorSloop*> (tmp->implementation());
351
352 ARCANE_CHECK_POINTER(x_sloop);
353 ARCANE_CHECK_POINTER(b_sloop);
354 ARCANE_CHECK_POINTER(tmp_sloop);
355
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;
359
360 if (isAlreadySolved(solution,RHS,temp,residual_norm,solver_param)){
361 debug() << "\t[AlephMatrixSloop::AlephMatrixSolve] isAlreadySolved !";
362 nb_iteration = 0;
363 return 0;
364 }
365
366 ScopedPtrT<SLOOP::SLOOPStopCriteria> sc;
367 sc=createSloopStopCriteria(solver_param, *sloopParallelInfo->m_sloop_msg);
368
369 ScopedPtrT<SLOOP::SLOOPSolver> global_solver;
370 global_solver=createSloopSolver(solver_param, *sloopParallelInfo->m_sloop_msg);
371
372 ScopedPtrT<SLOOP::SLOOPPreconditioner> precond;
373 precond=createSloopPreconditionner(solver_param, *sloopParallelInfo->m_sloop_msg);
374
375 this->setSloopSolverParameters(solver_param, global_solver.get());
376 this->setSloopPreconditionnerParameters(solver_param, precond.get());
377
378 ScopedPtrT<SLOOP::SLOOPDistVector> diag;
379 const bool normalize = normalizeSolverMatrix(solver_param);
380 if (normalize) {
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);
384 }
385
386 const bool xo = solver_param->xoUser();
387 switch (solver_param->precond()) {
388 case TypesSolver::NONE:
389 // appel sans preconditionnement : utilise dans le CAS des solveurs SAMG et SuperLU
390 if (xo){
391 debug() << "\t\t[AlephMatrixSloop::SloopSolve] xo à true (sans preconditionnement)";
392 status = global_solver->solve(*m_sloop_matrix, *solution, *RHS, *sc);
393 }else{
394 debug() << "\t\t[AlephMatrixSloop::SloopSolve] xo à false (sans preconditionnement)";
395 status = global_solver->solve_b(*m_sloop_matrix, *solution, *RHS, *sc);
396 }
397 break;
398 default:
399 // appel avec preconditionnement
400 if (xo){
401 debug() << "\t\t[AlephMatrixSloop::SloopSolve] xo à true (avec preconditionnement)";
402 status = global_solver->solve(*m_sloop_matrix, *solution, *RHS, *precond, *sc);
403 } else{
404 debug() << "\t\t[AlephMatrixSloop::SloopSolve] xo à false (avec preconditionnement)";
405 status = global_solver->solve_b(*m_sloop_matrix, *solution, *RHS, *precond, *sc);
406 }
407 break;
408 }
409
410 if (normalize) {
411 debug() << "\t\t[AlephMatrixSloop::SloopSolve] INV-normalize";
412 global_solver->inv_normalize(*m_sloop_matrix, *diag, *solution, *RHS);
413 }
414
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();
419
420 if ((solver_param->getCriteriaStop()==TypesSolver::STAG)||(solver_param->getCriteriaStop()==TypesSolver::NIter)){
421 // pas de controle des iterations dans les cas du critere de stagnation
422 // et du critere sur le nombre d'iterations impose du solveur
423 }else{
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");
432 }
433 }
434 debug() << "\t\t[AlephMatrixSloop::SloopSolve] nbIteration="<< global_solver->get_iteration()
435 << ", criteria=" << residual_norm[0] << ", stagnation=" << residual_norm[3];
436 return status;
437 }
438
439
440/******************************************************************************
441 * writeToFile
442 *****************************************************************************/
443 void writeToFile(const String file_name){
444 ItacFunction(AlephMatrixSloop);
445 m_sloop_matrix->write_to_file(file_name.localstr());
446 }
447
448
449/******************************************************************************
450 * createSloopSolver
451 *****************************************************************************/
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);
465 return solver;
466 }
467 default: throw Exception("AlephMatrixSloop::createSloopSolver", "Type de solver non accessible pour la bibliothèque SLOOP");
468 }
469 return NULL;
470 }
471
472
473/******************************************************************************
474 * createSloopPreconditionner
475 *****************************************************************************/
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:
482 // Pour le moment utilisation AINV(0) pour des matrices symetriques
483 // et SPAI pour les matrices non symetriques
484 if (TypesSolver::PCG == solver_param->method())
485 return new SLOOP::SLOOPAINVPC(info_sloop_msg, output_level);
486 else
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");
499 }
500 return NULL;
501 }
502
503
504/******************************************************************************
505 * createSloopStopCriteria
506 *****************************************************************************/
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();
521 }
522 return NULL;
523 }
524
525
526/******************************************************************************
527 *****************************************************************************/
528 void setSloopSolverParameters(AlephParams* solver_param,
529 SLOOP::SLOOPSolver* sloop_solver){
530 Integer error=0;
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; // TODOMB valeurs propres //error += sloop_solver->set_parameter(cg_spectrum_size, 50);
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);
545 break;
546 case TypesSolver::SAMG: {
547 error += sloop_solver->set_parameter(SLOOP::amg_buffer_size, 200);
548 // gamma niveaux maximum de deraffinement
549 if (gamma == -1) gamma = 50; // valeur par defaut du nombre de niveaux
550 error += sloop_solver->set_parameter(SLOOP::amg_level, gamma);
551 //(solver_param->gamma() == -1)?50:solver_param->gamma());//Works fine here
552 // valeur par defaut du parametre d'influence
553 if (alpha < 0.0) alpha = 0.1;
554 error += sloop_solver->set_parameter(SLOOP::amg_alpha, alpha);
555 //(solver_param->alpha() < 0.0)?0.1:solver_param->alpha());// Works here too
556 // nombre d'itérations du lisseur du solveur AMG
557 error += sloop_solver->set_parameter(SLOOP::amg_smoother_iter, solver_param->getAmgSmootherIter());
558 //1= cycle en V , 2=cycle en W, 3=cycle en FullMultigridV (defaut 1)
559 error += sloop_solver->set_parameter(SLOOP::amg_iter, solver_param->getAmgCycle());
560 // definition du deraffinement AMG
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);
564 //error += sloop_solver->set_parameter(SLOOP::amg_coarsening_option, solver_param->getAmgCoarseningOption());
565 // definition du solveur du deraffinement
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);
569 //error += sloop_solver->set_parameter(SLOOP::amg_coarse_solver_option, solver_param->getAmgCoarseSolverOption());
570 // definition du lisseur AMG
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);
574 //error += sloop_solver->set_parameter(SLOOP::amg_smoother_option, solver_param->getAmgSmootherOption());
575 }
576 break;
577 default: throw Exception("AlephMatrixSloop::setSloopSolverParameters",
578 "Type de solver SLOOP non prévu dans la gestion des paramètres");
579 }
580 if (error)
581 throw Exception("AlephMatrixSloop::setSloopSolverParameters", "set_parameter() failed");
582 }
583
584
585/******************************************************************************
586 *****************************************************************************/
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);
594 //valeur par defaut du parametre d'influence
595 if (alpha < 0.0) alpha = 0.1; // 0.25 ;
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);//(solver_param->gamma()==-1)?50:solver_param->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());
606 // option du lisseur
607 SLOOP::SLOOPAMGSmootherOption smoother_option =
608 static_cast<SLOOP::SLOOPAMGSmootherOption> (solver_param->getAmgSmootherOption());
609 // option du deraffinement
610 SLOOP::SLOOPAMGCoarseningOption coarsening_option =
611 static_cast<SLOOP::SLOOPAMGCoarseningOption> (solver_param->getAmgCoarseningOption());
612 // choix du solveur du systeme grossier defaut (CG_coarse_solver) pour une matrice symetrique
613 SLOOP::SLOOPAMGCoarseSolverOption coarse_solver_option =
614 static_cast<SLOOP::SLOOPAMGCoarseSolverOption> (solver_param->getAmgCoarseSolverOption());
615 // options pour matrice symetrique
616 if (TypesSolver::PCG == solver_param->method()) {
617 // controle du smoother pour matrice symetrique
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:
625 break;
626 default: throw Exception("AlephMatrixSloop::setSloopPreconditionnerParameters",
627 "choix du smoother incorrect pour une matrice symetrique");
628 }
629 }else{ //options pour matrices non symetrique
630 // choix du lisseur dans le cas non-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:
640 // on modifie la valeur mise par defaut pour AMG
641 solver_param->setAmgSmootherOption(TypesSolver::HybGSJ_smoother);
642 break;
643 default: break;
644 }
645 // controle du solveur
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);
650 break;
651 default:solver_param->setAmgCoarseSolverOption(TypesSolver::BiCGStab_coarse_solver); // choix du solveur du systeme grossier
652 break;
653 }
654 } // fin du if matrice symetrique
655 // definition du lisseur AMG apres controle
656 preconditionner->set_parameter(SLOOP::amg_smoother_option, smoother_option);//solver_param->getAmgSmootherOption());
657 // definition du deraffinement AMG
658 preconditionner->set_parameter(SLOOP::amg_coarsening_option, coarsening_option);//solver_param->getAmgCoarseningOption());
659 // definition du solveur du systeme grossier apres controle
660 preconditionner->set_parameter(SLOOP::amg_coarse_solver_option, coarse_solver_option);//solver_param->getAmgCoarseSolverOption());
661 }
662 break;
663
664 case TypesSolver::POLY:
665 if (gamma == -1) gamma = 3; // la valeur pour l'ordre du polynome
666 break;
667 case TypesSolver::AINV:
668 if (gamma == -1) gamma = 0; // valeur par defaut pour le parametre de remplissage -> AINV0
669 // le système linéaire doit être normalisé
670 break;
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);
680 break;
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);
686 break;
687 default:
688 throw Exception("AlephMatrixSloop::setSloopPreconditionnerParameters", "Préconditionneur inconnu.");
689 }
690
691 // initialisations communes aux preconditionneurs ( sauf sans precond )
692 switch (precond_method) {
693 case TypesSolver::NONE:
694 case TypesSolver::SPAIstat:
695 case TypesSolver::SPAIdyn:
696 break;
697
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);
702 // definition du critere d'arret du solveur grossier (default RB)
703 preconditionner->set_parameter(SLOOP::amg_coarse_solver_sc_option,SLOOP::RB);
704 break;
705
706 default:
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);
711 break;
712 }
713
714 }
715
716
717/******************************************************************************
718 *****************************************************************************/
719 bool normalizeSolverMatrix(AlephParams* solver_param){
720 ItacFunction(AlephMatrixSloop);
721 switch (solver_param->precond()){
722 case TypesSolver::AINV:
723 //case TypesSolver::AMG:
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.");
734 }
735 return false;
736 }
737public:
738 SLOOP::SLOOPMatrix* m_sloop_matrix;
739};
740
741/*---------------------------------------------------------------------------*/
742/*---------------------------------------------------------------------------*/
743
744class SloopAlephFactoryImpl: public AbstractService,
745 public IAlephFactoryImpl{
746public:
747 SloopAlephFactoryImpl(const ServiceBuildInfo& sbi) :
748 AbstractService(sbi),
749 m_IAlephVectors(0),
750 m_IAlephMatrixs(0),
751 m_IAlephTopologys(0){}
752 ~SloopAlephFactoryImpl(){
753 debug() << "\33[1;32m[~SloopAlephFactoryImpl]\33[0m";
754 for(Integer i=0,iMax=m_IAlephVectors.size(); i<iMax; ++i)
755 delete m_IAlephVectors.at(i);
756 for(Integer i=0,iMax=m_IAlephMatrixs.size(); i<iMax; ++i)
757 delete m_IAlephMatrixs.at(i);
758 for(Integer i=0,iMax=m_IAlephTopologys.size(); i<iMax; ++i)
759 delete m_IAlephTopologys.at(i);
760 }
761public:
762 virtual void initialize() {}
763
764 virtual IAlephTopology* createTopology(ITraceMng* tm,
765 AlephKernel* kernel,
766 Integer index,
767 Integer nb_row_size){
768 IAlephTopology *new_topology=new AlephTopologySloop(tm, kernel, index, nb_row_size);
769 m_IAlephTopologys.add(new_topology);
770 return new_topology;
771 }
772
773 virtual IAlephVector* createVector(ITraceMng* tm,
774 AlephKernel* kernel,
775 Integer index){
776 IAlephVector *new_vector=new AlephVectorSloop(tm,kernel,index);
777 m_IAlephVectors.add(new_vector);
778 return new_vector;
779 }
780
781 virtual IAlephMatrix* createMatrix(ITraceMng* tm,
782 AlephKernel* kernel,
783 Integer index){
784 IAlephMatrix *new_matrix=new AlephMatrixSloop(tm,kernel,index);
785 m_IAlephMatrixs.add(new_matrix);
786 return new_matrix;
787 }
788private:
789 UniqueArray<IAlephVector*> m_IAlephVectors;
790 UniqueArray<IAlephMatrix*> m_IAlephMatrixs;
791 UniqueArray<IAlephTopology*> m_IAlephTopologys;
792};
793
794ARCANE_REGISTER_APPLICATION_FACTORY(SloopAlephFactoryImpl,IAlephFactoryImpl,SloopAlephFactory);
795
796/*---------------------------------------------------------------------------*/
797/*---------------------------------------------------------------------------*/
798
799}
800
801/*---------------------------------------------------------------------------*/
802/*---------------------------------------------------------------------------*/
ARCANE_CHECK_POINTER
#define ARCANE_CHECK_POINTER(ptr)
Macro retournant le pointeur ptr s'il est non nul ou lancant une exception s'il est nul.
Definition ArcaneGlobal.h:823
ARCANE_REGISTER_APPLICATION_FACTORY
#define ARCANE_REGISTER_APPLICATION_FACTORY(aclass, ainterface, aname)
Enregistre un service de fabrique pour la classe aclass.
Definition ServiceFactory.h:572
Arcane::AbstractService::AbstractService
AbstractService(const ServiceBuildInfo &)
Constructeur à partir d'un ServiceBuildInfo.
Definition AbstractService.cc:29
Arcane::AlephParams
Paramètres d'un système linéraire.
Definition AlephParams.h:34
Arcane::AlephVector
Vecteur d'un système linéaire.
Definition AlephVector.h:33
Arcane::Exception
Classe de base d'une exception.
Definition arccore/src/base/arccore/base/Exception.h:48
Arcane::FatalErrorException
Exception lorsqu'une erreur fatale est survenue.
Definition arccore/src/base/arccore/base/FatalErrorException.h:32
Arcane::IAlephFactoryImpl
Interface d'une fabrique d'implémentation pour Aleph.
Definition AlephInterface.h:189
Arcane::ITraceMng
Interface du gestionnaire de traces.
Definition arccore/src/trace/arccore/trace/ITraceMng.h:156
Arcane::PtrT::get
T * get() const
Retourne l'objet référé par l'instance.
Definition Ptr.h:121
Arcane::ScopedPtrT
Encapsulation d'un pointeur qui se détruit automatiquement.
Definition ScopedPtr.h:44
Arcane::ServiceBuildInfo
Structure contenant les informations pour créer un service.
Definition ServiceBuildInfo.h:198
Arcane::String
Chaîne de caractères unicode.
Definition arccore/src/base/arccore/base/String.h:70
Arcane::String::localstr
const char * localstr() const
Retourne la conversion de l'instance dans l'encodage UTF-8.
Definition String.cc:227
Arcane::TraceAccessor::debug
TraceMessageDbg debug(Trace::eDebugLevel=Trace::Medium) const
Flot pour un message de debug.
Definition arccore/src/trace/arccore/trace/TraceAccessor.h:120
Arcane::TraceAccessor::info
TraceMessage info() const
Flot pour un message d'information.
Definition TraceAccessor.cc:101
Arcane::TraceAccessor::error
TraceMessage error() const
Flot pour un message d'erreur.
Definition TraceAccessor.cc:201
Arcane::UniqueArray
Vecteur 1D de données avec sémantique par valeur (style STL).
Definition arccore/src/collections/arccore/collections/Array.h:1844
Arcane
-*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
Definition AcceleratorGlobal.h:36
Arcane::Integer
Int32 Integer
Type représentant un entier.
Definition ArccoreGlobal.h:240
Arcane::Real
double Real
Type représentant un réel.
Definition ArccoreGlobal.h:223