Arcane  4.1.12.0
Developer documentation
Loading...
Searching...
No Matches
AlephSloop.cc
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2026 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// * Class AlephTopologySloop
39// ****************************************************************************
40class AlephTopologySloop : public IAlephTopology
41{
42 public:
43
44 AlephTopologySloop(ITraceMng* tm,
45 AlephKernel* kernel,
46 Integer index,
47 Integer nb_row_size)
48 : IAlephTopology(tm, kernel, index, nb_row_size)
49 , m_sloop_comminfo(NULL)
50 , m_world_comminfo(NULL)
51 , m_sloop_msg(NULL)
52 , m_sloop_topology(NULL)
53 {
54 ItacFunction(AlephTopologySloop);
55
56 if (!m_kernel->isParallel()) {
57 m_world_comminfo = new SLOOP::SLOOPSeqCommInfo();
58 }
59 else {
60 m_world_comminfo = new SLOOP::SLOOPMPICommInfo(MPI_COMM_WORLD);
61 }
62
63 if (!m_participating_in_solver) {
64 debug() << "\33[1;32m\t[AlephTopologySloop::AlephParallelInfoSloop] Not concerned with this solver, returning\33[0m";
65 return;
66 }
67
68 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] @" << this << "\33[0m";
69
70 if (!m_kernel->isParallel()) {
71 m_sloop_comminfo = new SLOOP::SLOOPSeqCommInfo();
72 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] SEQCommInfo @" << m_sloop_comminfo << "\33[0m";
73 }
74 else {
75 m_sloop_comminfo = new SLOOP::SLOOPMPICommInfo(*(SLOOP::SLOOP_Comm*)(m_kernel->subParallelMng(index)->getMPICommunicator()));
76 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] MPICommInfo @" << m_sloop_comminfo << "\33[0m";
77 }
78
79 m_sloop_msg = new SLOOP::SLOOPMsg(m_sloop_comminfo);
80 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] SLOOPMsg @" << m_sloop_msg << "\33[0m";
81
82 m_sloop_topology = new SLOOP::SLOOPTopology(nb_row_size,
83 m_kernel->topology()->rowLocalRange(index),
84 *m_sloop_comminfo, *m_sloop_msg);
85 debug() << "\33[1;32m\t\t[AlephTopologySloop::AlephTopologySloop] SLOOPTopology @" << m_sloop_topology << "\33[0m";
86 }
87 ~AlephTopologySloop()
88 {
89 debug() << "\33[1;5;32m\t\t\t[~AlephTopologySloop] deleting m_sloop_msg, m_sloop_comminfo & m_sloop_topology\33[0m";
90 delete m_sloop_msg;
91 delete m_sloop_comminfo;
92 delete m_sloop_topology;
93 }
94
95 public:
96
97 // We back up the current session and initialize the new session with our m_sloop_comminfo
98 void backupAndInitialize()
99 {
100 if (!m_participating_in_solver)
101 return;
102 SLOOP::SLOOPInitSession(m_sloop_comminfo, SLOOP::WITHOUT_EXCEPTIONS, SLOOP::OUTPUT_UNIQ);
103 }
104 void restore()
105 {
106 // In any case, everyone restores the WORLD session
107 SLOOP::SLOOPInitSession(m_world_comminfo, SLOOP::WITHOUT_EXCEPTIONS, SLOOP::OUTPUT_UNIQ);
108 }
109
110 public:
111
112 SLOOP::SLOOPCommInfo* m_sloop_comminfo;
113 SLOOP::SLOOPCommInfo* m_world_comminfo;
114 SLOOP::SLOOPMsg* m_sloop_msg;
115 SLOOP::SLOOPTopology* m_sloop_topology;
116};
117
118// ****************************************************************************
119// * AlephVectorSloop
120// ****************************************************************************
121class AlephVectorSloop
122: public IAlephVector
123{
124 public:
125
126 AlephVectorSloop(ITraceMng* tm,
127 AlephKernel* kernel,
128 Integer index)
129 : IAlephVector(tm, kernel, index)
130 , m_sloop_vector(NULL)
131 {
132 debug() << "\t\t[AlephVectorSloop::AlephVectorSloop] NEW AlephVectorSloop";
133 }
134
135 // ****************************************************************************
136 // * ~AlephVectorSloop
137 // ****************************************************************************
138 ~AlephVectorSloop()
139 {
140 debug() << "\33[1;5;32m\t\t\t[~AlephVectorSloop]\33[0m";
141 delete m_sloop_vector;
142 }
143
144 /******************************************************************************
145 * AlephVectorCreate
146 *****************************************************************************/
147 void AlephVectorCreate(void)
148 {
149 ItacFunction(AlephVectorSloop);
150 debug() << "\t[AlephVectorSloop::AlephVectorCreate] new SLOOP::SLOOPDistVector";
151 AlephTopologySloop* aleph_topology_sloop = dynamic_cast<AlephTopologySloop*>(m_kernel->getTopologyImplementation(m_index));
152 ARCANE_CHECK_POINTER(aleph_topology_sloop);
153 m_sloop_vector = new SLOOP::SLOOPDistVector(*aleph_topology_sloop->m_sloop_topology,
154 *aleph_topology_sloop->m_sloop_msg);
155 if (!m_sloop_vector)
156 throw FatalErrorException(A_FUNCINFO, " new SLOOPDistVector failed");
157 debug() << "\t[AlephVectorSloop::AlephVectorCreate] done";
158 }
159
160 // ****************************************************************************
161 // * AlephVectorSet
162 // ****************************************************************************
163 void AlephVectorSet(const double* bfr_val, const int* bfr_idx, Integer size)
164 {
165 debug() << "\t[AlephVectorSloop::AlephVectorSet]";
166 if (m_sloop_vector->locfill(bfr_val, bfr_idx, size))
167 throw FatalErrorException(A_FUNCINFO, "locfill() failed");
168 }
169
170 /******************************************************************************
171 * AlephVectorAssemble
172 *****************************************************************************/
173 int AlephVectorAssemble(void)
174 {
175 ItacFunction(AlephVectorSloop);
176 debug() << "\t\t[AlephVectorSloop::AlephVectorAssemble]";
177 return 0;
178 }
179
180 /******************************************************************************
181 * AlephVectorGet
182 *****************************************************************************/
183 void AlephVectorGet(double* bfr_val, const int* bfr_idx, Integer size)
184 {
185 // I am already retrieving the values from the Sloop vector.
186 debug() << "\t[AlephVectorSloop::AlephVectorGet]";
187 if (m_sloop_vector->get_locval(bfr_val, bfr_idx, size))
188 throw Exception("AlephVectorSloop::AlephVectorGet", "get_locval() failed");
189 }
190
191 /******************************************************************************
192 * writeToFile
193 *****************************************************************************/
194 void writeToFile(const String file_name)
195 {
196 ItacFunction(AlephVectorSloop);
197 m_sloop_vector->write_to_file(file_name.localstr());
198 }
199
200 public:
201
202 SLOOP::SLOOPDistVector* m_sloop_vector;
203};
204
205/******************************************************************************
206 AlephMatrixSloop
207*****************************************************************************/
208class AlephMatrixSloop : public IAlephMatrix
209{
210 public:
211
212 /******************************************************************************
213 * AlephMatrixSloop
214 *****************************************************************************/
215 AlephMatrixSloop(ITraceMng* tm,
216 AlephKernel* kernel,
217 Integer index)
218 : IAlephMatrix(tm, kernel, index)
219 , m_sloop_matrix(NULL)
220 {
221 debug() << "\t\t[AlephMatrixSloop::AlephMatrixSloop] NEW AlephMatrixSloop";
222 }
223
224 // ****************************************************************************
225 // * AlephMatrixSloop
226 // ****************************************************************************
227 ~AlephMatrixSloop()
228 {
229 debug() << "\33[1;5;32m\t\t\t[~AlephMatrixSloop]\33[0m";
230 delete m_sloop_matrix;
231 }
232
233 /******************************************************************************
234 * AlephMatrixCreate
235 *****************************************************************************/
236 void AlephMatrixCreate(void)
237 {
238 ItacFunction(AlephMatrixSloop);
239 debug() << "\t\t[AlephMatrixSloop::AlephMatrixCreate] create new SLOOP::SLOOPDistMatrix";
240 AlephTopologySloop* aleph_topology_sloop = dynamic_cast<AlephTopologySloop*>(m_kernel->getTopologyImplementation(m_index));
241 ARCANE_CHECK_POINTER(aleph_topology_sloop);
242 m_sloop_matrix = new SLOOP::SLOOPDistMatrix(*aleph_topology_sloop->m_sloop_topology,
243 *aleph_topology_sloop->m_sloop_msg,
244 true,
245 false);
246
247 if (!m_sloop_matrix)
248 throw Exception("AlephSolverMatrix::create", "new SLOOPDistMatrix() failed");
249
250 // Renumbering must be done before filling the matrix to be taken into account.
251 if (aleph_topology_sloop->m_sloop_topology->get_type() == SLOOP::contiguous)
252 m_sloop_matrix->set_renumbering_opt(SLOOP::SLOOPDistMatrix::interface, false);
253 m_sloop_matrix->set_renumbering_opt(SLOOP::SLOOPDistMatrix::processor, false);
254 m_sloop_matrix->set_renumbering_opt(SLOOP::SLOOPDistMatrix::interior, false);
255 m_sloop_matrix->init_length(m_kernel->topology()->gathered_nb_row_elements().unguardedBasePointer());
256 debug() << "\t\t[AlephMatrixSloop::AlephMatrixCreate] done";
257 }
258
259 /******************************************************************************
260 * AlephMatrixSetFilled
261 *****************************************************************************/
262 void AlephMatrixSetFilled(bool toggle)
263 {
264 ItacFunction(AlephMatrixSloop);
265 debug() << "\t\t[AlephMatrixSloop::AlephMatrixSetFilled]";
266 m_sloop_matrix->setfilled(toggle);
267 }
268
269 /******************************************************************************
270 * AlephMatrixAssemble
271 *****************************************************************************/
272 int AlephMatrixAssemble(void)
273 {
274 ItacFunction(AlephMatrixSloop);
275 debug() << "\t\t[AlephMatrixSloop::AlephMatrixAssemble]";
276 m_sloop_matrix->configure();
277 return 0;
278 }
279
280 /******************************************************************************
281 * AlephMatrixFill
282 *****************************************************************************/
283 void AlephMatrixFill(int size, int* rows, int* cols, double* values)
284 {
285 ItacFunction(AlephMatrixSloop);
286 debug() << "\t\t[AlephMatrixSloop::AlephMatrixFill] size=" << size;
287 m_sloop_matrix->locfill(values, rows, cols, size);
288 debug() << "\t\t[AlephMatrixSloop::AlephMatrixFill] done";
289 }
290
291 /******************************************************************************
292 * isAlreadySolved
293 *****************************************************************************/
294 bool isAlreadySolved(SLOOP::SLOOPDistVector* x,
295 SLOOP::SLOOPDistVector* b,
296 SLOOP::SLOOPDistVector* tmp,
297 Real* residual_norm,
298 AlephParams* params)
299 {
300 const Real res0 = b->norm_max();
301 const Real considered_as_null = params->minRHSNorm();
302 const bool convergence_analyse = params->convergenceAnalyse();
303
304 if (convergence_analyse)
305 debug() << "convergence analysis: max norm of the right-hand side res0: " << res0;
306
307 if (res0 < considered_as_null) {
308 x->fill(Real(0.0));
309 residual_norm[0] = res0;
310 if (convergence_analyse)
311 debug() << "convergence analysis: the right-hand side of the linear system is less than "
312 << considered_as_null;
313 return true;
314 }
315
316 if (params->xoUser()) {
317 // We test if b is already a solution to the system within epsilon tolerance
318 m_sloop_matrix->vector_multiply(*tmp, *x); // tmp=A*x
319 tmp->substract(*tmp, *b); // tmp=A*x-b
320 const Real residu = tmp->norm_max();
321 debug() << "[IAlephSloop::isAlreadySolved] residu=" << residu;
322
323 if (residu < considered_as_null) {
324 if (convergence_analyse) {
325 debug() << "convergence analysis: |Ax0-b| is less than " << considered_as_null;
326 debug() << "convergence analysis: x0 is already a solution to the system.";
327 }
328 residual_norm[0] = residu;
329 return true;
330 }
331 const Real relative_error = residu / res0;
332 if (convergence_analyse)
333 debug() << "convergence analysis: initial residual: " << residu
334 << " --- initial relative residual (residu/res0): " << residu / res0;
335
336 if (relative_error < (params->epsilon())) {
337 if (convergence_analyse)
338 debug() << "convergence analysis: X is already a solution to the system";
339 residual_norm[0] = residu;
340 return true;
341 }
342 }
343 return false;
344 }
345
346 /******************************************************************************
347 *****************************************************************************/
348 int AlephMatrixSolve(AlephVector* x,
349 AlephVector* b,
350 AlephVector* tmp,
351 Integer& nb_iteration,
352 Real* residual_norm,
353 AlephParams* solver_param)
354 {
355 ItacFunction(AlephMatrixSloop);
356 Integer status = 0;
357 debug() << "\t\t[AlephMatrixSloop::SloopSolve] getTopologyImplementation #" << m_index;
358 AlephTopologySloop* sloopParallelInfo = dynamic_cast<AlephTopologySloop*>(m_kernel->getTopologyImplementation(m_index));
359
360 AlephVectorSloop* x_sloop = dynamic_cast<AlephVectorSloop*>(x->implementation());
361 AlephVectorSloop* b_sloop = dynamic_cast<AlephVectorSloop*>(b->implementation());
362 AlephVectorSloop* tmp_sloop = dynamic_cast<AlephVectorSloop*>(tmp->implementation());
363
364 ARCANE_CHECK_POINTER(x_sloop);
365 ARCANE_CHECK_POINTER(b_sloop);
366 ARCANE_CHECK_POINTER(tmp_sloop);
367
368 SLOOP::SLOOPDistVector* solution = x_sloop->m_sloop_vector;
369 SLOOP::SLOOPDistVector* RHS = b_sloop->m_sloop_vector;
370 SLOOP::SLOOPDistVector* temp = tmp_sloop->m_sloop_vector;
371
372 if (isAlreadySolved(solution, RHS, temp, residual_norm, solver_param)) {
373 debug() << "\t[AlephMatrixSloop::AlephMatrixSolve] isAlreadySolved !";
374 nb_iteration = 0;
375 return 0;
376 }
377
378 ScopedPtrT<SLOOP::SLOOPStopCriteria> sc;
379 sc = createSloopStopCriteria(solver_param, *sloopParallelInfo->m_sloop_msg);
380
381 ScopedPtrT<SLOOP::SLOOPSolver> global_solver;
382 global_solver = createSloopSolver(solver_param, *sloopParallelInfo->m_sloop_msg);
383
384 ScopedPtrT<SLOOP::SLOOPPreconditioner> precond;
385 precond = createSloopPreconditionner(solver_param, *sloopParallelInfo->m_sloop_msg);
386
387 this->setSloopSolverParameters(solver_param, global_solver.get());
388 this->setSloopPreconditionnerParameters(solver_param, precond.get());
389
390 ScopedPtrT<SLOOP::SLOOPDistVector> diag;
391 const bool normalize = normalizeSolverMatrix(solver_param);
392 if (normalize) {
393 debug() << "\t\t[AlephMatrixSloop::SloopSolve] Normalize";
394 diag = new SLOOP::SLOOPDistVector(*sloopParallelInfo->m_sloop_topology, *sloopParallelInfo->m_sloop_msg);
395 global_solver->normalize(*m_sloop_matrix, *diag, *solution, *RHS);
396 }
397
398 const bool xo = solver_param->xoUser();
399 switch (solver_param->precond()) {
400 case TypesSolver::NONE:
401 // call without preconditioning: used in the case of SAMG and SuperLU solvers
402 if (xo) {
403 debug() << "\t\t[AlephMatrixSloop::SloopSolve] xo à true (without preconditioning)";
404 status = global_solver->solve(*m_sloop_matrix, *solution, *RHS, *sc);
405 }
406 else {
407 debug() << "\t\t[AlephMatrixSloop::SloopSolve] xo à false (without preconditioning)";
408 status = global_solver->solve_b(*m_sloop_matrix, *solution, *RHS, *sc);
409 }
410 break;
411 default:
412 // call with preconditioning
413 if (xo) {
414 debug() << "\t\t[AlephMatrixSloop::SloopSolve] xo à true (with preconditioning)";
415 status = global_solver->solve(*m_sloop_matrix, *solution, *RHS, *precond, *sc);
416 }
417 else {
418 debug() << "\t\t[AlephMatrixSloop::SloopSolve] xo à false (avec preconditionnement)";
419 status = global_solver->solve_b(*m_sloop_matrix, *solution, *RHS, *precond, *sc);
420 }
421 break;
422 }
423
424 if (normalize) {
425 debug() << "\t\t[AlephMatrixSloop::SloopSolve] INV-normalize";
426 global_solver->inv_normalize(*m_sloop_matrix, *diag, *solution, *RHS);
427 }
428
429 nb_iteration = global_solver->get_iteration();
430 residual_norm[0] = sc->get_criteria();
431 Integer max_iteration = global_solver->get_max_iteration();
432 residual_norm[3] = global_solver->get_stagnation();
433
434 if ((solver_param->getCriteriaStop() == TypesSolver::STAG) || (solver_param->getCriteriaStop() == TypesSolver::NIter)) {
435 // no iteration control in cases of stagnation criterion
436 // and the criterion on the number of iterations imposed by the solver
437 }
438 else {
439 if (nb_iteration == max_iteration && solver_param->stopErrorStrategy()) {
440 info() << "\n============================================================";
441 info() << "\nCette erreur est retournée après " << nb_iteration << "\n";
442 info() << "\nOn a atteind le nombre max d'itérations du solveur ";
443 info() << "\nIl est possible de demander au code de ne pas tenir compte de cette erreur.";
444 info() << "\nVoir la documentation du jeu de données concernant le service solveur.";
445 info() << "\n============================================================";
446 throw Exception("AlephMatrixSloop::SloopSolve", "On a atteind le nombre max d'itérations du solveur");
447 }
448 }
449 debug() << "\t\t[AlephMatrixSloop::SloopSolve] nbIteration=" << global_solver->get_iteration()
450 << ", criteria=" << residual_norm[0] << ", stagnation=" << residual_norm[3];
451 return status;
452 }
453
454 /******************************************************************************
455 * writeToFile
456 *****************************************************************************/
457 void writeToFile(const String file_name)
458 {
459 ItacFunction(AlephMatrixSloop);
460 m_sloop_matrix->write_to_file(file_name.localstr());
461 }
462
463 /******************************************************************************
464 * createSloopSolver
465 *****************************************************************************/
466 SLOOP::SLOOPSolver* createSloopSolver(AlephParams* solver_param, SLOOP::SLOOPMsg& info_sloop_msg)
467 {
468 const Integer output_level = solver_param->getOutputLevel();
469 ItacFunction(AlephMatrixSloop);
470 switch (solver_param->method()) {
471 case TypesSolver::PCG:
472 return new SLOOP::SLOOPCGSolver(info_sloop_msg, output_level);
473 case TypesSolver::BiCGStab:
474 return new SLOOP::SLOOPBiCGStabSolver(info_sloop_msg, output_level);
475 case TypesSolver::BiCGStab2:
476 return new SLOOP::SLOOPBiCGStabSolver(info_sloop_msg, output_level);
477 case TypesSolver::GMRES:
478 return new SLOOP::SLOOPGMRESSolver(info_sloop_msg, output_level);
479 case TypesSolver::SAMG:
480 return new SLOOP::SLOOPAMGSolver(info_sloop_msg, output_level);
481 case TypesSolver::QMR:
482 return new SLOOP::SLOOPQMRSolver(info_sloop_msg, output_level);
483 case TypesSolver::SuperLU: {
484 SLOOP::SLOOPSolver* solver = new SLOOP::SLOOPSuperLUSolver(info_sloop_msg, output_level);
485 solver->set_parameter(SLOOP::sv_residual_calculation, 1);
486 return solver;
487 }
488 default:
489 throw Exception("AlephMatrixSloop::createSloopSolver", "Type de solver non accessible pour la bibliothèque SLOOP");
490 }
491 return NULL;
492 }
493
494 /******************************************************************************
495 * createSloopPreconditionner
496 *****************************************************************************/
497 SLOOP::SLOOPPreconditioner* createSloopPreconditionner(AlephParams* solver_param,
498 SLOOP::SLOOPMsg& info_sloop_msg)
499 {
500 const Integer output_level = solver_param->getOutputLevel();
501 ItacFunction(AlephMatrixSloop);
502 switch (solver_param->precond()) {
503 case TypesSolver::AINV:
504 // For now, AINV(0) is used for symmetric matrices
505 // and SPAI for non-symmetric matrices
506 if (TypesSolver::PCG == solver_param->method())
507 return new SLOOP::SLOOPAINVPC(info_sloop_msg, output_level);
508 else
509 return new SLOOP::SLOOPMAINVPC(info_sloop_msg, output_level);
510 case TypesSolver::DIAGONAL:
511 return new SLOOP::SLOOPDiagPC(info_sloop_msg, output_level);
512 case TypesSolver::AMG:
513 return new SLOOP::SLOOPAMGPC(info_sloop_msg, output_level);
514 case TypesSolver::IC:
515 return new SLOOP::SLOOPCholPC(info_sloop_msg, output_level);
516 case TypesSolver::POLY:
517 return new SLOOP::SLOOPPolyPC(info_sloop_msg, output_level);
518 case TypesSolver::ILU:
519 return new SLOOP::SLOOPILUPC(info_sloop_msg, output_level);
520 case TypesSolver::ILUp:
521 return new SLOOP::SLOOPILUPC(info_sloop_msg, output_level);
522 case TypesSolver::SPAIstat:
523 return new SLOOP::SLOOPSPAIPC(info_sloop_msg, output_level);
524 case TypesSolver::SPAIdyn:
525 return new SLOOP::SLOOPSPAIPC(info_sloop_msg, output_level);
526 case TypesSolver::NONE:
527 return NULL;
528 default:
529 throw Exception("AlephMatrixSloop::createSloopPreconditionner",
530 "préconditionneur non accessible pour la bibliothèque SLOOP");
531 }
532 return NULL;
533 }
534
535 /******************************************************************************
536 * createSloopStopCriteria
537 *****************************************************************************/
538 SLOOP::SLOOPStopCriteria* createSloopStopCriteria(AlephParams* solver_param,
539 SLOOP::SLOOPMsg& info_sloop_msg)
540 {
541 ItacFunction(AlephMatrixSloop);
542 switch (solver_param->getCriteriaStop()) {
543 case TypesSolver::RR0:
544 return new SLOOP::SLOOPStopCriteriaRR0();
545 case TypesSolver::R:
546 return new SLOOP::SLOOPStopCriteriaR();
547 case TypesSolver::RCB:
548 return new SLOOP::SLOOPStopCriteriaRCB();
549 case TypesSolver::RBinf:
550 return new SLOOP::SLOOPStopCriteriaRBinf();
551 case TypesSolver::EpsA:
552 return new SLOOP::SLOOPStopCriteriaEpsA();
553 case TypesSolver::NIter:
554 return new SLOOP::SLOOPStopCriteriaNIter();
555 case TypesSolver::RR0inf:
556 return new SLOOP::SLOOPStopCriteriaRR0inf();
557 case TypesSolver::STAG:
558 return new SLOOP::SLOOPStopCriteriaSTAG();
559 case TypesSolver::RB:
560 default:
561 return new SLOOP::SLOOPStopCriteriaRB();
562 }
563 return NULL;
564 }
565
566 /******************************************************************************
567 *****************************************************************************/
568 void setSloopSolverParameters(AlephParams* solver_param,
569 SLOOP::SLOOPSolver* sloop_solver)
570 {
571 Integer error = 0;
572 Integer gamma = solver_param->gamma();
573 double alpha = solver_param->alpha();
574 ItacFunction(AlephMatrixSloop);
575 error += sloop_solver->set_parameter(SLOOP::sv_epsilon, solver_param->epsilon());
576 error += sloop_solver->set_parameter(SLOOP::sv_max_iter, solver_param->maxIter());
577 switch (solver_param->method()) {
578 case TypesSolver::PCG:
579 break; // TODOMB eigenvalues //error += sloop_solver->set_parameter(cg_spectrum_size, 50);
580 case TypesSolver::QMR:
581 break;
582 case TypesSolver::SuperLU:
583 break;
584 case TypesSolver::BiCGStab:
585 error += sloop_solver->set_parameter(SLOOP::bicg_dimension, 1);
586 break;
587 case TypesSolver::BiCGStab2:
588 error += sloop_solver->set_parameter(SLOOP::bicg_dimension, 2);
589 break;
590 case TypesSolver::GMRES:
591 error += sloop_solver->set_parameter(SLOOP::gmres_order, 20);
592 error += sloop_solver->set_parameter(SLOOP::gmres_type, SLOOP::ICGS);
593 break;
594 case TypesSolver::SAMG: {
595 error += sloop_solver->set_parameter(SLOOP::amg_buffer_size, 200);
596 // gamma maximum deraffinment levels
597 if (gamma == -1)
598 gamma = 50; // default value for the number of levels
599 error += sloop_solver->set_parameter(SLOOP::amg_level, gamma);
600 //(solver_param->gamma() == -1)?50:solver_param->gamma());//Works fine here
601 // default value for the influence parameter
602 if (alpha < 0.0)
603 alpha = 0.1;
604 error += sloop_solver->set_parameter(SLOOP::amg_alpha, alpha);
605 //(solver_param->alpha() < 0.0)?0.1:solver_param->alpha());// Works here too
606 // number of iterations of the AMG solver smoother
607 error += sloop_solver->set_parameter(SLOOP::amg_smoother_iter, solver_param->getAmgSmootherIter());
608 //1= V cycle, 2= W cycle, 3= FullMultigridV cycle (default 1)
609 error += sloop_solver->set_parameter(SLOOP::amg_iter, solver_param->getAmgCycle());
610 // definition of AMG deraffinment
611 SLOOP::SLOOPAMGCoarseningOption coarsening_option =
612 static_cast<SLOOP::SLOOPAMGCoarseningOption>(solver_param->getAmgCoarseningOption());
613 error += sloop_solver->set_parameter(SLOOP::amg_coarsening_option, coarsening_option);
614 //error += sloop_solver->set_parameter(SLOOP::amg_coarsening_option, solver_param->getAmgCoarseningOption());
615 // definition of the deraffinment solver
616 SLOOP::SLOOPAMGCoarseSolverOption coarse_solver_option =
617 static_cast<SLOOP::SLOOPAMGCoarseSolverOption>(solver_param->getAmgCoarseSolverOption());
618 error += sloop_solver->set_parameter(SLOOP::amg_coarse_solver_option, coarse_solver_option);
619 //error += sloop_solver->set_parameter(SLOOP::amg_coarse_solver_option, solver_param->getAmgCoarseSolverOption());
620 // definition of the AMG smoother
621 SLOOP::SLOOPAMGSmootherOption smoother_option =
622 static_cast<SLOOP::SLOOPAMGSmootherOption>(solver_param->getAmgSmootherOption());
623 error += sloop_solver->set_parameter(SLOOP::amg_smoother_option, smoother_option);
624 //error += sloop_solver->set_parameter(SLOOP::amg_smoother_option, solver_param->getAmgSmootherOption());
625 } break;
626 default:
627 throw Exception("AlephMatrixSloop::setSloopSolverParameters",
628 "Type de solver SLOOP non prévu dans la gestion des paramètres");
629 }
630 if (error)
631 throw Exception("AlephMatrixSloop::setSloopSolverParameters", "set_parameter() failed");
632 }
633
634 /******************************************************************************
635 *****************************************************************************/
636 void setSloopPreconditionnerParameters(AlephParams* solver_param,
637 SLOOP::SLOOPPreconditioner* preconditionner)
638 {
639 const TypesSolver::ePreconditionerMethod precond_method = solver_param->precond();
640 double alpha = solver_param->alpha();
641 int gamma = solver_param->gamma();
642 const String function_id = "SolverMatrixSloop::setSloopPreconditionnerParameters";
643 ItacFunction(AlephMatrixSloop);
644 // default value of the influence parameter
645 if (alpha < 0.0)
646 alpha = 0.1; // 0.25 ;
647 switch (precond_method) {
648 case TypesSolver::NONE:
649 break;
650 case TypesSolver::DIAGONAL:
651 break;
652 case TypesSolver::AMG: {
653 if (gamma == -1)
654 gamma = 50;
655 preconditionner->set_parameter(SLOOP::amg_level, gamma); //(solver_param->gamma()==-1)?50:solver_param->gamma());
656 preconditionner->set_parameter(SLOOP::amg_buffer_size, 200);
657 preconditionner->set_parameter(SLOOP::amg_solver_iter, solver_param->getAmgSolverIter());
658 preconditionner->set_parameter(SLOOP::amg_iter, solver_param->getAmgCycle());
659 preconditionner->set_parameter(SLOOP::amg_smoother_iter, solver_param->getAmgSmootherIter());
660 // smoother option
661 SLOOP::SLOOPAMGSmootherOption smoother_option =
662 static_cast<SLOOP::SLOOPAMGSmootherOption>(solver_param->getAmgSmootherOption());
663 // deraffinment option
664 SLOOP::SLOOPAMGCoarseningOption coarsening_option =
665 static_cast<SLOOP::SLOOPAMGCoarseningOption>(solver_param->getAmgCoarseningOption());
666 // choice of coarse system solver default (CG_coarse_solver) for a symmetric matrix
667 SLOOP::SLOOPAMGCoarseSolverOption coarse_solver_option =
668 static_cast<SLOOP::SLOOPAMGCoarseSolverOption>(solver_param->getAmgCoarseSolverOption());
669 // options for symmetric matrix
670 if (TypesSolver::PCG == solver_param->method()) {
671 // smoother control for symmetric matrix
672 switch (solver_param->getAmgSmootherOption()) {
673 case TypesSolver::CG_smoother:
674 case TypesSolver::Rich_IC_smoother:
675 case TypesSolver::Rich_AINV_smoother:
676 case TypesSolver::SymHybGSJ_smoother:
677 case TypesSolver::Rich_IC_block_smoother:
678 case TypesSolver::SymHybGSJ_block_smoother:
679 break;
680 default:
681 throw Exception("AlephMatrixSloop::setSloopPreconditionnerParameters",
682 "incorrect smoother choice for a symmetric matrix");
683 }
684 }
685 else { //options for non-symmetric matrices
686 // choice of smoother in the non-symmetric case
687 switch (solver_param->getAmgSmootherOption()) {
688 case TypesSolver::CG_smoother:
689 case TypesSolver::Rich_IC_smoother:
690 case TypesSolver::Rich_AINV_smoother:
691 case TypesSolver::Rich_IC_block_smoother:
692 case TypesSolver::SymHybGSJ_block_smoother:
693 throw Exception("AlephMatrixSloop::setSloopPreconditionnerParameters",
694 "incorrect smoother choice with a non-symmetric matrix ");
695 case TypesSolver::SymHybGSJ_smoother:
696 // we modify the default value for AMG
697 solver_param->setAmgSmootherOption(TypesSolver::HybGSJ_smoother);
698 break;
699 default:
700 break;
701 }
702 // coarse solver control
703 switch (solver_param->getAmgCoarseSolverOption()) {
704 case TypesSolver::CG_coarse_solver:
705 case TypesSolver::Cholesky_coarse_solver:
706 solver_param->setAmgCoarseSolverOption(TypesSolver::LU_coarse_solver);
707 break;
708 default:
709 solver_param->setAmgCoarseSolverOption(TypesSolver::BiCGStab_coarse_solver); // choice of coarse system solver
710 break;
711 }
712 } // end of symmetric matrix if
713 // definition of AMG smoother after control
714 preconditionner->set_parameter(SLOOP::amg_smoother_option, smoother_option); //solver_param->getAmgSmootherOption());
715 // definition of AMG deraffinment
716 preconditionner->set_parameter(SLOOP::amg_coarsening_option, coarsening_option); //solver_param->getAmgCoarseningOption());
717 // definition of the coarse system solver after control
718 preconditionner->set_parameter(SLOOP::amg_coarse_solver_option, coarse_solver_option); //solver_param->getAmgCoarseSolverOption());
719 } break;
720
721 case TypesSolver::POLY:
722 if (gamma == -1)
723 gamma = 3; // the value for the polynomial order
724 break;
725 case TypesSolver::AINV:
726 if (gamma == -1)
727 gamma = 0; // default value for the fill-in parameter -> AINV0
728 // the linear system must be normalized
729 break;
730 case TypesSolver::IC:
731 case TypesSolver::ILU:
732 if (gamma == -1)
733 gamma = 0;
734 break;
735 case TypesSolver::ILUp:
736 if (gamma == -1)
737 gamma = 1;
738 break;
739 case TypesSolver::SPAIstat:
740 preconditionner->set_parameter(SLOOP::spai_sparsity, SLOOP::StatSparsity);
741 preconditionner->set_parameter(SLOOP::spai_init_sparsity, SLOOP::PowerSparsity);
742 preconditionner->set_parameter(SLOOP::spai_power_level, 1);
743 preconditionner->set_parameter(SLOOP::spai_Amax_row_size, 30);
744 preconditionner->set_parameter(SLOOP::spai_A_drop_eps, 0.001);
745 break;
746 case TypesSolver::SPAIdyn:
747 preconditionner->set_parameter(SLOOP::spai_sparsity, SLOOP::DynSparsity);
748 preconditionner->set_parameter(SLOOP::spai_init_sparsity, SLOOP::DiagSparsity);
749 preconditionner->set_parameter(SLOOP::spai_Amax_row_size, 10);
750 preconditionner->set_parameter(SLOOP::spai_A_drop_eps, 0.001);
751 break;
752 default:
753 throw Exception("AlephMatrixSloop::setSloopPreconditionnerParameters", "Préconditionneur inconnu.");
754 }
755
756 // common initializations for preconditioners (except without precond)
757 switch (precond_method) {
758 case TypesSolver::NONE:
759 case TypesSolver::SPAIstat:
760 case TypesSolver::SPAIdyn:
761 break;
762
763 case TypesSolver::AMG:
764 preconditionner->set_parameter(SLOOP::amg_alpha, alpha);
765 preconditionner->set_parameter(SLOOP::amg_level, gamma);
766 preconditionner->set_parameter(SLOOP::amg_parallel_opt, 0);
767 // definition of the coarse solver stopping criterion (default RB)
768 preconditionner->set_parameter(SLOOP::amg_coarse_solver_sc_option, SLOOP::RB);
769 break;
770
771 default:
772 preconditionner->set_parameter(SLOOP::pc_alpha, alpha);
773 preconditionner->set_parameter(SLOOP::pc_nbelem, gamma);
774 preconditionner->set_parameter(SLOOP::pc_parallel_opt, 1);
775 preconditionner->set_parameter(SLOOP::pc_order, gamma);
776 break;
777 }
778 }
779
780 /******************************************************************************
781 *****************************************************************************/
782 bool normalizeSolverMatrix(AlephParams* solver_param)
783 {
784 ItacFunction(AlephMatrixSloop);
785 switch (solver_param->precond()) {
786 case TypesSolver::AINV:
787 //case TypesSolver::AMG:
788 case TypesSolver::SPAIstat:
789 case TypesSolver::SPAIdyn:
790 return true;
791 case TypesSolver::AMG:
792 case TypesSolver::NONE:
793 case TypesSolver::DIAGONAL:
794 case TypesSolver::IC:
795 case TypesSolver::POLY:
796 case TypesSolver::ILU:
797 case TypesSolver::ILUp:
798 return false;
799 default:
800 throw Exception("AlephMatrixSloop::normalizeSolverMatrix", "Préconditionneur inconnu.");
801 }
802 return false;
803 }
804
805 public:
806
807 SLOOP::SLOOPMatrix* m_sloop_matrix;
808};
809
810/*---------------------------------------------------------------------------*/
811/*---------------------------------------------------------------------------*/
812
813class SloopAlephFactoryImpl : public AbstractService
814, public IAlephFactoryImpl
815{
816 public:
817
818 SloopAlephFactoryImpl(const ServiceBuildInfo& sbi)
819 : AbstractService(sbi)
820 , m_IAlephVectors(0)
821 , m_IAlephMatrixs(0)
822 , m_IAlephTopologys(0)
823 {}
824 ~SloopAlephFactoryImpl()
825 {
826 debug() << "\33[1;32m[~SloopAlephFactoryImpl]\33[0m";
827 for (Integer i = 0, iMax = m_IAlephVectors.size(); i < iMax; ++i)
828 delete m_IAlephVectors.at(i);
829 for (Integer i = 0, iMax = m_IAlephMatrixs.size(); i < iMax; ++i)
830 delete m_IAlephMatrixs.at(i);
831 for (Integer i = 0, iMax = m_IAlephTopologys.size(); i < iMax; ++i)
832 delete m_IAlephTopologys.at(i);
833 }
834
835 public:
836
837 virtual void initialize() {}
838
839 virtual IAlephTopology* createTopology(ITraceMng* tm,
840 AlephKernel* kernel,
841 Integer index,
842 Integer nb_row_size)
843 {
844 IAlephTopology* new_topology = new AlephTopologySloop(tm, kernel, index, nb_row_size);
845 m_IAlephTopologys.add(new_topology);
846 return new_topology;
847 }
848
849 virtual IAlephVector* createVector(ITraceMng* tm,
850 AlephKernel* kernel,
851 Integer index)
852 {
853 IAlephVector* new_vector = new AlephVectorSloop(tm, kernel, index);
854 m_IAlephVectors.add(new_vector);
855 return new_vector;
856 }
857
858 virtual IAlephMatrix* createMatrix(ITraceMng* tm,
859 AlephKernel* kernel,
860 Integer index)
861 {
862 IAlephMatrix* new_matrix = new AlephMatrixSloop(tm, kernel, index);
863 m_IAlephMatrixs.add(new_matrix);
864 return new_matrix;
865 }
866
867 private:
868
869 UniqueArray<IAlephVector*> m_IAlephVectors;
870 UniqueArray<IAlephMatrix*> m_IAlephMatrixs;
871 UniqueArray<IAlephTopology*> m_IAlephTopologys;
872};
873
874ARCANE_REGISTER_APPLICATION_FACTORY(SloopAlephFactoryImpl, IAlephFactoryImpl, SloopAlephFactory);
875
876/*---------------------------------------------------------------------------*/
877/*---------------------------------------------------------------------------*/
878
879} // namespace Arcane
880
881/*---------------------------------------------------------------------------*/
882/*---------------------------------------------------------------------------*/
ARCANE_CHECK_POINTER
#define ARCANE_CHECK_POINTER(ptr)
Macro returning the pointer ptr if it is not null or throwing an exception if it is null.
Definition ArcaneGlobal.h:819
ARCANE_REGISTER_APPLICATION_FACTORY
#define ARCANE_REGISTER_APPLICATION_FACTORY(aclass, ainterface, aname)
Registers a factory service for the class aclass.
Definition ServiceFactory.h:636
Arcane::AbstractService::AbstractService
AbstractService(const ServiceBuildInfo &)
Constructor from a ServiceBuildInfo.
Definition AbstractService.cc:30
Arcane::AlephParams
Parameters of a linear system.
Definition AlephParams.h:32
Arcane::AlephVector
Vector of a linear system.
Definition AlephVector.h:33
Arcane::IAlephFactoryImpl
Interface of an implementation factory for Aleph.
Definition AlephInterface.h:188
Arcane::PtrT::get
T * get() const
Returns the object referenced by the instance.
Definition Ptr.h:122
Arcane::ScopedPtrT
Encapsulation of an automatically destructing pointer.
Definition ScopedPtr.h:44
Arcane::ServiceBuildInfo
Structure containing the information to create a service.
Definition ServiceBuildInfo.h:199
Arcane::String::localstr
const char * localstr() const
Returns the conversion of the instance into UTF-8 encoding.
Definition String.cc:229
Arcane::TraceAccessor::debug
TraceMessageDbg debug(Trace::eDebugLevel=Trace::Medium) const
Flow for a debug message.
Definition arccore/src/trace/arccore/trace/TraceAccessor.h:122
Arcane::TraceAccessor::info
TraceMessage info() const
Flow for an information message.
Definition TraceAccessor.cc:101
Arcane::TraceAccessor::error
TraceMessage error() const
Flow for an error message.
Definition TraceAccessor.cc:201
Arcane::UniqueArray
1D data vector with value semantics (STL style).
Definition arccore/src/common/arccore/common/Array.h:890
Arcane
-- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature --
Definition arcane/src/arcane/accelerator/AcceleratorGlobal.h:37
Arcane::Integer
Int32 Integer
Type representing an integer.
Definition ArccoreGlobal.h:293
Arcane::Real
double Real
Type representing a real number.
Definition ArccoreGlobal.h:275