26#include "./HypreComparer.h"
32#include "HYPRE_krylov.h"
34#include "HYPRE_parcsr_ls.h"
36#include <HYPRE_config.h>
38#if defined(HYPRE_EXAMPLE_USING_CUDA)
40#include <cuda_runtime.h>
42#ifndef HYPRE_USING_UNIFIED_MEMORY
43#error *** Running the examples on GPUs requires Unified Memory. Please reconfigure and rebuild with --enable-unified-memory ***
47gpu_malloc(
size_t size)
50 cudaMallocManaged(&ptr, size, cudaMemAttachGlobal);
55gpu_calloc(
size_t num,
size_t size)
58 cudaMallocManaged(&ptr, num * size, cudaMemAttachGlobal);
59 cudaMemset(ptr, 0, num * size);
63#define malloc(size) gpu_malloc(size)
64#define calloc(num, size) gpu_calloc(num, size)
65#define free(ptr) ( cudaFree(ptr), ptr = NULL )
72#include "arccore/base/Convert.h"
73#include "arccore/base/FatalErrorException.h"
74#include "arccore/alina/Profiler.h"
78int hypre_FlexGMRESModifyPCAMGExample(
void *precond_data,
int iterations,
79 double rel_residual_norm);
81#define my_min(a,b) (((a)<(b)) ? (a) : (b))
98 std::vector<ptrdiff_t>
const& _ptr,
99 std::vector<ptrdiff_t>
const& _col,
100 std::vector<double>
const& _val,
101 std::vector<double>
const& _rhs,
102 std::vector<double>& _x,
103 int argc,
char* argv[])
105 auto& prof = Alina::Profiler::globalProfiler();
106 auto t = prof.scoped_tic(
"Hypre");
108 std::cout <<
"DO_HYPRE nb_row=" << nb_row <<
"\n";
113 int local_size, extra;
119 HYPRE_ParCSRMatrix parcsr_A;
121 HYPRE_ParVector par_b;
123 HYPRE_ParVector par_x;
125 HYPRE_Solver solver, precond;
128 if (m_do_mpi_init_and_finalize) {
129 MPI_Init(&argc, &argv);
130 m_need_finalize =
true;
133 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
134 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
141#if defined(HYPRE_USING_GPU)
143 HYPRE_SetSpGemmUseVendor(0);
147 const int n = nb_row;
156 while (arg_index < argc) {
157 if (strcmp(argv[arg_index],
"-n") == 0) {
161 else if (strcmp(argv[arg_index],
"-solver") == 0) {
163 solver_id = atoi(argv[arg_index++]);
165 else if (strcmp(argv[arg_index],
"-print_system") == 0) {
169 else if (strcmp(argv[arg_index],
"-help") == 0) {
178 if ((print_usage) && (myid == 0)) {
180 printf(
"Usage: %s [<options>]\n", argv[0]);
182 printf(
" -n <n> : problem size in each direction (default: 33)\n");
183 printf(
" -solver <ID> : solver ID\n");
184 printf(
" 0 - AMG (default) \n");
185 printf(
" 1 - AMG-PCG\n");
186 printf(
" 8 - ParaSails-PCG\n");
187 printf(
" 50 - PCG\n");
188 printf(
" 61 - AMG-FlexGMRES\n");
189 printf(
" -vis : save the solution for GLVis visualization\n");
190 printf(
" -print_system : print the matrix and rhs\n");
200 std::vector<int> nb_value_per_row(n);
201 for (
int i = 0; i < n; ++i)
202 nb_value_per_row[i] =
static_cast<HYPRE_BigInt
>(_ptr[i + 1] - _ptr[i]);
206 std::vector<HYPRE_BigInt> hypre_column_index(_col.begin(), _col.end());
209 std::vector<HYPRE_Int> hypre_row_index(n);
210 for (
int i = 0; i < n; ++i) {
211 hypre_row_index[i] = i;
217 local_size = nb_row / num_procs;
218 extra = nb_row - local_size * num_procs;
220 ilower = local_size * myid;
221 ilower += my_min(myid, extra);
223 iupper = local_size * (myid + 1);
224 iupper += my_min(myid + 1, extra);
226 std::cout <<
"LOWER=" << ilower <<
" UPPER=" << iupper <<
"\n";
228 local_size = iupper - ilower + 1;
231 auto t = prof.scoped_tic(
"IJMatrix Create");
235 HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A);
238 HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR);
241 HYPRE_IJMatrixInitialize(A);
246 auto t = prof.scoped_tic(
"IJMatrix SetValues");
248 HYPRE_IJMatrixSetValues(A, n,
249 nb_value_per_row.data(),
250 hypre_row_index.data(),
251 hypre_column_index.data(),
256 auto t = prof.scoped_tic(
"IJMatrix Assemble");
258 HYPRE_IJMatrixAssemble(A);
275 HYPRE_IJMatrixGetObject(A, (
void**)&parcsr_A);
278 HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper, &b);
279 HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR);
280 HYPRE_IJVectorInitialize(b);
282 HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper, &x);
283 HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR);
284 HYPRE_IJVectorInitialize(x);
293 rows = (
int*)calloc(local_size,
sizeof(
int));
295 for (i = 0; i < local_size; i++) {
298 rows[i] = ilower + i;
301 HYPRE_IJVectorSetValues(b, local_size, rows, _rhs.data());
302 HYPRE_IJVectorSetValues(x, local_size, rows, _x.data());
309 HYPRE_IJVectorAssemble(b);
317 HYPRE_IJVectorGetObject(b, (
void**)&par_b);
319 HYPRE_IJVectorAssemble(x);
320 HYPRE_IJVectorGetObject(x, (
void**)&par_x);
325 HYPRE_IJMatrixPrint(A,
"IJ.out.A");
326 HYPRE_IJVectorPrint(b,
"IJ.out.b");
329 if (
auto v = Convert::Type<Int32>::tryParseFromEnvironment(
"ALINA_HYPRE_SOLVER",
true))
330 solver_id = v.value();
332 double solver_tolerance = 1.0e-8;
335 std::cout <<
"FINISH ASSEMBLING solver_id=" << solver_id <<
"\n";
337 if (solver_id == 0) {
338 auto t = prof.scoped_tic(
"HypreSolver AMG");
340 double final_res_norm;
343 HYPRE_BoomerAMGCreate(&solver);
346 HYPRE_BoomerAMGSetPrintLevel(solver, 3);
347 HYPRE_BoomerAMGSetOldDefault(solver);
348 HYPRE_BoomerAMGSetRelaxType(solver, 3);
349 HYPRE_BoomerAMGSetRelaxOrder(solver, 1);
350 HYPRE_BoomerAMGSetNumSweeps(solver, 1);
351 HYPRE_BoomerAMGSetMaxLevels(solver, 20);
352 HYPRE_BoomerAMGSetTol(solver, solver_tolerance);
356 auto t = prof.scoped_tic(
"AMG Setup");
357 HYPRE_BoomerAMGSetup(solver, parcsr_A, par_b, par_x);
360 auto t = prof.scoped_tic(
"AMG Solve");
361 HYPRE_BoomerAMGSolve(solver, parcsr_A, par_b, par_x);
365 HYPRE_BoomerAMGGetNumIterations(solver, &num_iterations);
366 HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, &final_res_norm);
369 printf(
"Iterations = %d\n", num_iterations);
370 printf(
"Final Relative Residual Norm = %e\n", final_res_norm);
375 HYPRE_BoomerAMGDestroy(solver);
378 else if (solver_id == 50) {
379 auto t = prof.scoped_tic(
"HypreSolver PCG");
381 double final_res_norm;
384 HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
387 HYPRE_PCGSetMaxIter(solver, 1000);
388 HYPRE_PCGSetTol(solver, solver_tolerance);
389 HYPRE_PCGSetTwoNorm(solver, 1);
390 HYPRE_PCGSetPrintLevel(solver, 2);
391 HYPRE_PCGSetLogging(solver, 1);
394 HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x);
395 HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);
398 HYPRE_PCGGetNumIterations(solver, &num_iterations);
399 HYPRE_PCGGetFinalRelativeResidualNorm(solver, &final_res_norm);
402 printf(
"Iterations = %d\n", num_iterations);
403 printf(
"Final Relative Residual Norm = %e\n", final_res_norm);
408 HYPRE_ParCSRPCGDestroy(solver);
411 else if (solver_id == 1) {
412 auto t = prof.scoped_tic(
"HypreSolver PCG-AMG");
414 double final_res_norm;
417 HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
420 HYPRE_PCGSetMaxIter(solver, 1000);
421 HYPRE_PCGSetTol(solver, solver_tolerance);
422 HYPRE_PCGSetTwoNorm(solver, 1);
423 HYPRE_PCGSetPrintLevel(solver, 2);
424 HYPRE_PCGSetLogging(solver, 1);
427 HYPRE_BoomerAMGCreate(&precond);
428 HYPRE_BoomerAMGSetPrintLevel(precond, 1);
429 HYPRE_BoomerAMGSetCoarsenType(precond, 6);
430 HYPRE_BoomerAMGSetOldDefault(precond);
431 HYPRE_BoomerAMGSetRelaxType(precond, 6);
432 HYPRE_BoomerAMGSetNumSweeps(precond, 1);
433 HYPRE_BoomerAMGSetTol(precond, 0.0);
434 HYPRE_BoomerAMGSetMaxIter(precond, 1);
437 HYPRE_PCGSetPrecond(solver, (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSolve,
438 (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSetup, precond);
442 HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x);
445 HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);
449 HYPRE_PCGGetNumIterations(solver, &num_iterations);
450 HYPRE_PCGGetFinalRelativeResidualNorm(solver, &final_res_norm);
453 printf(
"Iterations = %d\n", num_iterations);
454 printf(
"Final Relative Residual Norm = %e\n", final_res_norm);
459 HYPRE_ParCSRPCGDestroy(solver);
460 HYPRE_BoomerAMGDestroy(precond);
463 else if (solver_id == 8) {
464 auto t = prof.scoped_tic(
"HypreSolver PCG - Parasails");
466 double final_res_norm;
468 int sai_max_levels = 1;
469 double sai_threshold = 0.1;
470 double sai_filter = 0.05;
474 HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
477 HYPRE_PCGSetMaxIter(solver, 1000);
478 HYPRE_PCGSetTol(solver, solver_tolerance);
479 HYPRE_PCGSetTwoNorm(solver, 1);
480 HYPRE_PCGSetPrintLevel(solver, 2);
481 HYPRE_PCGSetLogging(solver, 1);
484 HYPRE_ParaSailsCreate(MPI_COMM_WORLD, &precond);
487 HYPRE_ParaSailsSetParams(precond, sai_threshold, sai_max_levels);
488 HYPRE_ParaSailsSetFilter(precond, sai_filter);
489 HYPRE_ParaSailsSetSym(precond, sai_sym);
490 HYPRE_ParaSailsSetLogging(precond, 3);
493 HYPRE_PCGSetPrecond(solver, (HYPRE_PtrToSolverFcn)HYPRE_ParaSailsSolve,
494 (HYPRE_PtrToSolverFcn)HYPRE_ParaSailsSetup, precond);
497 HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x);
498 HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);
501 HYPRE_PCGGetNumIterations(solver, &num_iterations);
502 HYPRE_PCGGetFinalRelativeResidualNorm(solver, &final_res_norm);
505 printf(
"Iterations = %d\n", num_iterations);
506 printf(
"Final Relative Residual Norm = %e\n", final_res_norm);
511 HYPRE_ParCSRPCGDestroy(solver);
512 HYPRE_ParaSailsDestroy(precond);
515 else if (solver_id == 61) {
516 auto t = prof.scoped_tic(
"HypreSolver Flexible GMRES - AMG");
518 double final_res_norm;
523 HYPRE_ParCSRFlexGMRESCreate(MPI_COMM_WORLD, &solver);
526 HYPRE_FlexGMRESSetKDim(solver, restart);
527 HYPRE_FlexGMRESSetMaxIter(solver, 1000);
528 HYPRE_FlexGMRESSetTol(solver, solver_tolerance);
529 HYPRE_FlexGMRESSetPrintLevel(solver, 2);
530 HYPRE_FlexGMRESSetLogging(solver, 1);
533 HYPRE_BoomerAMGCreate(&precond);
534 HYPRE_BoomerAMGSetPrintLevel(precond, 1);
535 HYPRE_BoomerAMGSetCoarsenType(precond, 6);
536 HYPRE_BoomerAMGSetOldDefault(precond);
537 HYPRE_BoomerAMGSetRelaxType(precond, 6);
538 HYPRE_BoomerAMGSetNumSweeps(precond, 1);
539 HYPRE_BoomerAMGSetTol(precond, 0.0);
540 HYPRE_BoomerAMGSetMaxIter(precond, 1);
543 HYPRE_FlexGMRESSetPrecond(solver, (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSolve,
544 (HYPRE_PtrToSolverFcn)HYPRE_BoomerAMGSetup, precond);
550 HYPRE_FlexGMRESSetModifyPC(solver, (HYPRE_PtrToModifyPCFcn)hypre_FlexGMRESModifyPCAMGExample);
555 auto t = prof.scoped_tic(
"FlexGMRES Setup");
556 HYPRE_ParCSRFlexGMRESSetup(solver, parcsr_A, par_b, par_x);
559 auto t = prof.scoped_tic(
"FlexGMRES Solve");
560 HYPRE_ParCSRFlexGMRESSolve(solver, parcsr_A, par_b, par_x);
564 HYPRE_FlexGMRESGetNumIterations(solver, &num_iterations);
565 HYPRE_FlexGMRESGetFinalRelativeResidualNorm(solver, &final_res_norm);
568 printf(
"Iterations = %d\n", num_iterations);
569 printf(
"Final Relative Residual Norm = %e\n", final_res_norm);
574 HYPRE_ParCSRFlexGMRESDestroy(solver);
575 HYPRE_BoomerAMGDestroy(precond);
579 ARCCORE_FATAL(
"Invalid solver id '{0}' specified.", solver_id);
584 HYPRE_IJVectorPrint(x,
"IJ.out.x");
587 HYPRE_IJMatrixDestroy(A);
588 HYPRE_IJVectorDestroy(b);
589 HYPRE_IJVectorDestroy(x);
608int hypre_FlexGMRESModifyPCAMGExample(
void* precond_data, [[maybe_unused]]
int iterations,
609 double rel_residual_norm)
612 if (rel_residual_norm > .1) {
613 HYPRE_BoomerAMGSetNumSweeps((HYPRE_Solver)precond_data, 10);
616 HYPRE_BoomerAMGSetNumSweeps((HYPRE_Solver)precond_data, 1);
#define ARCCORE_FATAL(...)
Macro throwing a FatalErrorException.
-- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature --