From 1e3ae70aecb44ff926a7d9a7e107872ba94ae812 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 21 Jun 2023 11:38:37 -0700 Subject: [PATCH 1/7] Add doxygen string for HYPRE_ParCSRMatrixMatvec --- src/parcsr_mv/protos.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parcsr_mv/protos.h b/src/parcsr_mv/protos.h index e8705e0804..a5b41c8ebb 100644 --- a/src/parcsr_mv/protos.h +++ b/src/parcsr_mv/protos.h @@ -65,6 +65,9 @@ HYPRE_Int HYPRE_CSRMatrixToParCSRMatrix ( MPI_Comm comm, HYPRE_CSRMatrix A_CSR, HYPRE_BigInt *row_partitioning, HYPRE_BigInt *col_partitioning, HYPRE_ParCSRMatrix *matrix ); HYPRE_Int HYPRE_CSRMatrixToParCSRMatrix_WithNewPartitioning ( MPI_Comm comm, HYPRE_CSRMatrix A_CSR, HYPRE_ParCSRMatrix *matrix ); +/** + * Performs y <- alpha * A * x + beta * y + */ HYPRE_Int HYPRE_ParCSRMatrixMatvec ( HYPRE_Complex alpha, HYPRE_ParCSRMatrix A, HYPRE_ParVector x, HYPRE_Complex beta, HYPRE_ParVector y ); HYPRE_Int HYPRE_ParCSRMatrixMatvecT ( HYPRE_Complex alpha, HYPRE_ParCSRMatrix A, HYPRE_ParVector x, From af0ca63fba94a0d25ec189e6937d209bbe684772 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Tue, 20 Jun 2023 17:42:10 -0700 Subject: [PATCH 2/7] Make an example demonstrating AIR for advection --- src/examples/Makefile | 8 +- src/examples/ex19.c | 357 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 364 insertions(+), 1 deletion(-) create mode 100644 src/examples/ex19.c diff --git a/src/examples/Makefile b/src/examples/Makefile index 5c75153d73..ccb4612693 100644 --- a/src/examples/Makefile +++ b/src/examples/Makefile @@ -18,7 +18,7 @@ LINK_CUF = ${CUF} XL_DIR=$(dir $(shell which xlc)).. -HYPRE_DIR = ../hypre +HYPRE_DIR ?= ../hypre ######################################################################## # CUDA @@ -281,6 +281,12 @@ ex18: ex18.o ex18comp: ex18comp.o $(LINK_CC) -o $@ $^ $(LFLAGS) +######################################################################## +# Example 19 +######################################################################## +ex19: ex19.o + $(LINK_CC) -o $@ $^ $(LFLAGS) + ######################################################################## # Clean up ######################################################################## diff --git a/src/examples/ex19.c b/src/examples/ex19.c new file mode 100644 index 0000000000..a118578c15 --- /dev/null +++ b/src/examples/ex19.c @@ -0,0 +1,357 @@ +/****************************************************************************** + * Copyright (c) 1998 Lawrence Livermore National Security, LLC and other + * HYPRE Project Developers. See the top-level COPYRIGHT file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + ******************************************************************************/ + +/* + Example 19 + + Interface: Linear-Algebraic (IJ) + + Compile with: make ex19 + + Sample run: mpirun -np 4 ex19 + + Description: This example solves a 1D steady upwinded advection problem using + AIR. With one up sweep and one coarse sweep, AIR converges in one + iteration. If the transpose of the interpolation operator is used + instead (-restriction 0), the multigrid solve fails to converge +*/ + + +#include +#include +#include +#include +#include "_hypre_utilities.h" +#include "HYPRE_krylov.h" +#include "HYPRE.h" +#include "HYPRE_parcsr_ls.h" +#include "ex.h" + +#define my_min(a,b) (((a)<(b)) ? (a) : (b)) + +int main (int argc, char *argv[]) +{ + int i; + int myid, num_procs; + int N, n = 33; + + int ilower, iupper; + int local_size, extra; + + int print_system = 0; + + HYPRE_IJMatrix A; + HYPRE_ParCSRMatrix parcsr_A; + HYPRE_IJVector b; + HYPRE_ParVector par_b; + HYPRE_IJVector x; + HYPRE_ParVector par_x; + HYPRE_IJVector work; + HYPRE_ParVector par_work; + + HYPRE_Solver solver; + /* AMG */ + int num_iterations; + double final_res_norm; + HYPRE_Int interp_type = 100; /* 1-pt Interp */ + HYPRE_Int ns_coarse = 1, ns_down = 0, ns_up = 1; + HYPRE_Int **grid_relax_points = NULL; + HYPRE_Int coarse_threshold = 20; + HYPRE_Int agg_num_levels = 0; /* AIR does not support aggressive coarsening */ + HYPRE_Int restriction = 1; + + + /* Initialize MPI */ + MPI_Init(&argc, &argv); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + + /* Initialize HYPRE */ + HYPRE_Initialize(); + + /* Print GPU info */ + /* HYPRE_PrintDeviceInfo(); */ +#if defined(HYPRE_USING_GPU) + /* use vendor implementation for SpGEMM */ + HYPRE_SetSpGemmUseVendor(0); +#endif + + /* Parse command line */ + { + int arg_index = 0; + int print_usage = 0; + + while (arg_index < argc) + { + if ( strcmp(argv[arg_index], "-n") == 0 ) + { + arg_index++; + n = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-print_system") == 0 ) + { + arg_index++; + print_system = 1; + } + else if ( strcmp(argv[arg_index], "-help") == 0 ) + { + print_usage = 1; + break; + } + else if ( strcmp(argv[arg_index], "-ns_coarse") == 0) + { + arg_index++; + ns_coarse = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ns_down") == 0) + { + arg_index++; + ns_down = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ns_up") == 0) + { + arg_index++; + ns_up = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-restriction") == 0) + { + arg_index++; + restriction = atoi(argv[arg_index++]); + } + else + { + arg_index++; + } + } + + if ((print_usage) && (myid == 0)) + { + printf("\n"); + printf("Usage: %s []\n", argv[0]); + printf("\n"); + printf(" -n : problem size (default: 33)\n"); + printf(" -print_system : print the matrix, rhs, expected solution, and solution\n"); + printf(" -ns_coarse : the number of sweeps performed for the coarse cycle " + "(default: 1)\n"); + printf(" -ns_down : the number of sweeps performed for the down cycle " + "(default: 0)\n"); + printf(" -ns_up : the number of sweeps performed for the up cycle " + "(default: 1)\n"); + printf(" -restriction : defines the restriction operator, 0 for transpose of " + "interpolation, 1 for AIR-1, 2 for AIR-2 (default: 1)\n"); + printf("\n"); + } + + if (print_usage) + { + MPI_Finalize(); + return (0); + } + } + + /* Preliminaries: want at least one processor per row */ + N = n; /* global number of rows */ + + /* Each processor knows only of its own rows - the range is denoted by ilower + and upper. Here we partition the rows. We account for the fact that + N may not divide evenly by the number of processors. */ + local_size = N / num_procs; + extra = N - local_size * num_procs; + + ilower = local_size * myid; + ilower += my_min(myid, extra); + + iupper = local_size * (myid + 1); + iupper += my_min(myid + 1, extra); + iupper = iupper - 1; + + /* How many rows do I have? */ + local_size = iupper - ilower + 1; + + /* Create the matrix. + Note that this is a square matrix, so we indicate the row partition + size twice (since number of rows = number of cols) */ + HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A); + + /* Choose a parallel csr format storage (see the User's Manual) */ + HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR); + + /* Initialize before setting coefficients */ + HYPRE_IJMatrixInitialize(A); + + /* Set matrix entries corresponding to 1D upwind advection */ + { + int nnz; + /* OK to use constant-length arrays for CPUs + double values[2]; + int cols[2]; + */ + double *values = (double *) malloc(2 * sizeof(double)); + int *cols = (int *) malloc(2 * sizeof(int)); + int *tmp = (int *) malloc(2 * sizeof(int)); + + for (i = ilower; i <= iupper; i++) + { + if (i == 0) + { + nnz = 1; + cols[0] = 0; + values[0] = 1; + } + else + { + nnz = 2; + cols[0] = i - 1; + cols[1] = i; + values[0] = -1; + values[1] = 1; + } + + /* Set the values for row i */ + tmp[0] = nnz; + tmp[1] = i; + HYPRE_IJMatrixSetValues(A, 1, &tmp[0], &tmp[1], cols, values); + } + + free(values); + free(cols); + free(tmp); + } + + /* Assemble after setting the coefficients */ + HYPRE_IJMatrixAssemble(A); + + /* Note: for the testing of small problems, one may wish to read + in a matrix in IJ format (for the format, see the output files + from the -print_system option). + In this case, one would use the following routine: + HYPRE_IJMatrixRead( , MPI_COMM_WORLD, + HYPRE_PARCSR, &A ); + = IJ.A.out to read in what has been printed out + by -print_system (processor numbers are omitted). + A call to HYPRE_IJMatrixRead is an *alternative* to the + following sequence of HYPRE_IJMatrix calls: + Create, SetObjectType, Initialize, SetValues, and Assemble + */ + + + /* Get the parcsr matrix object to use */ + HYPRE_IJMatrixGetObject(A, (void**) &parcsr_A); + + + /* Create the solution, work vector, and b */ + HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper, &x); + HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR); + HYPRE_IJVectorInitialize(x); + + HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper, &work); + HYPRE_IJVectorSetObjectType(work, HYPRE_PARCSR); + HYPRE_IJVectorInitialize(work); + + HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper, &b); + HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR); + HYPRE_IJVectorInitialize(b); + + HYPRE_IJVectorGetObject(x, (void **) &par_x); + HYPRE_ParVectorSetRandomValues(par_x, 1); + + HYPRE_IJVectorGetObject(work, (void **) &par_work); + HYPRE_ParVectorSetRandomValues(par_work, 2); + + HYPRE_IJVectorGetObject(b, (void **) &par_b); + HYPRE_ParVectorSetConstantValues(par_b, 0); + HYPRE_ParCSRMatrixMatvec(1, parcsr_A, par_work, 0, par_b); + + + /* Print out the system - files names will be IJ.out.A.XXXXX + and IJ.out.b.XXXXX, where XXXXX = processor id */ + if (print_system) + { + HYPRE_IJMatrixPrint(A, "IJ.out.A"); + HYPRE_IJVectorPrint(work, "IJ.out.work"); + HYPRE_IJVectorPrint(b, "IJ.out.b"); + } + + /* this is a 2-D 4-by-k array using Double pointers */ + grid_relax_points = hypre_CTAlloc(HYPRE_Int*, 4, HYPRE_MEMORY_HOST); + grid_relax_points[0] = NULL; + grid_relax_points[1] = hypre_CTAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST); + grid_relax_points[2] = hypre_CTAlloc(HYPRE_Int, ns_up, HYPRE_MEMORY_HOST); + grid_relax_points[3] = hypre_CTAlloc(HYPRE_Int, ns_coarse, HYPRE_MEMORY_HOST); + /* down cycle: C */ + for (i = 0; i < ns_down; i++) + { + grid_relax_points[1][i] = 1; // C + } + /* up cycle: F */ + for (i=0; i Date: Wed, 21 Jun 2023 15:51:38 -0700 Subject: [PATCH 3/7] Use classical interpolation if not doing AIR --- src/examples/ex19.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/examples/ex19.c b/src/examples/ex19.c index a118578c15..5d73cbe3a3 100644 --- a/src/examples/ex19.c +++ b/src/examples/ex19.c @@ -128,6 +128,10 @@ int main (int argc, char *argv[]) } } + if (!restriction) + // Do classical modified interpolation if not doing AIR + interp_type = 0; + if ((print_usage) && (myid == 0)) { printf("\n"); @@ -303,7 +307,6 @@ int main (int argc, char *argv[]) /* Set some parameters (See Reference Manual for more parameters) */ HYPRE_BoomerAMGSetPrintLevel(solver, 3); /* print solve info + parameters */ - HYPRE_BoomerAMGSetOldDefault(solver); /* Falgout coarsening with modified classical interpolaiton */ HYPRE_BoomerAMGSetRelaxType(solver, 0); /* Jacobi */ HYPRE_BoomerAMGSetRelaxOrder(solver, 1); /* uses C/F relaxation */ HYPRE_BoomerAMGSetCycleNumSweeps(solver, ns_down, 1); From cf014b7a106a8aedee4052b948090c1a36de1ef0 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Sun, 25 Jun 2023 23:46:44 -0400 Subject: [PATCH 4/7] Remove comment from protos.h file --- src/parcsr_mv/protos.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/parcsr_mv/protos.h b/src/parcsr_mv/protos.h index a5b41c8ebb..e8705e0804 100644 --- a/src/parcsr_mv/protos.h +++ b/src/parcsr_mv/protos.h @@ -65,9 +65,6 @@ HYPRE_Int HYPRE_CSRMatrixToParCSRMatrix ( MPI_Comm comm, HYPRE_CSRMatrix A_CSR, HYPRE_BigInt *row_partitioning, HYPRE_BigInt *col_partitioning, HYPRE_ParCSRMatrix *matrix ); HYPRE_Int HYPRE_CSRMatrixToParCSRMatrix_WithNewPartitioning ( MPI_Comm comm, HYPRE_CSRMatrix A_CSR, HYPRE_ParCSRMatrix *matrix ); -/** - * Performs y <- alpha * A * x + beta * y - */ HYPRE_Int HYPRE_ParCSRMatrixMatvec ( HYPRE_Complex alpha, HYPRE_ParCSRMatrix A, HYPRE_ParVector x, HYPRE_Complex beta, HYPRE_ParVector y ); HYPRE_Int HYPRE_ParCSRMatrixMatvecT ( HYPRE_Complex alpha, HYPRE_ParCSRMatrix A, HYPRE_ParVector x, From 7e1bfe7f929e7f3985fc435914f39154f89cebfd Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Sun, 25 Jun 2023 23:49:05 -0400 Subject: [PATCH 5/7] Remove HYPRE_ data types --- src/examples/ex19.c | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/examples/ex19.c b/src/examples/ex19.c index 5d73cbe3a3..52b8c28dc5 100644 --- a/src/examples/ex19.c +++ b/src/examples/ex19.c @@ -54,16 +54,16 @@ int main (int argc, char *argv[]) HYPRE_ParVector par_work; HYPRE_Solver solver; - /* AMG */ - int num_iterations; - double final_res_norm; - HYPRE_Int interp_type = 100; /* 1-pt Interp */ - HYPRE_Int ns_coarse = 1, ns_down = 0, ns_up = 1; - HYPRE_Int **grid_relax_points = NULL; - HYPRE_Int coarse_threshold = 20; - HYPRE_Int agg_num_levels = 0; /* AIR does not support aggressive coarsening */ - HYPRE_Int restriction = 1; + /* AMG */ + int num_iterations; + double final_res_norm; + int interp_type = 100; /* 1-pt Interp */ + int ns_coarse = 1, ns_down = 0, ns_up = 1; + int **grid_relax_points = NULL; + int coarse_threshold = 20; + int agg_num_levels = 0; /* AIR does not support aggressive coarsening */ + int restriction = 1; /* Initialize MPI */ MPI_Init(&argc, &argv); @@ -281,7 +281,7 @@ int main (int argc, char *argv[]) } /* this is a 2-D 4-by-k array using Double pointers */ - grid_relax_points = hypre_CTAlloc(HYPRE_Int*, 4, HYPRE_MEMORY_HOST); + grid_relax_points = hypre_CTAlloc(int*, 4, HYPRE_MEMORY_HOST); grid_relax_points[0] = NULL; grid_relax_points[1] = hypre_CTAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST); grid_relax_points[2] = hypre_CTAlloc(HYPRE_Int, ns_up, HYPRE_MEMORY_HOST); From afd4ee5061530d20833617533c5f5d2fdccb33ee Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Sun, 25 Jun 2023 23:57:20 -0400 Subject: [PATCH 6/7] Style changes and update comment about grid_relax_points --- src/examples/ex19.c | 43 +++++++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/src/examples/ex19.c b/src/examples/ex19.c index 52b8c28dc5..91efe99b84 100644 --- a/src/examples/ex19.c +++ b/src/examples/ex19.c @@ -128,9 +128,11 @@ int main (int argc, char *argv[]) } } + /* Do classical modified interpolation if not doing AIR */ if (!restriction) - // Do classical modified interpolation if not doing AIR + { interp_type = 0; + } if ((print_usage) && (myid == 0)) { @@ -200,20 +202,20 @@ int main (int argc, char *argv[]) for (i = ilower; i <= iupper; i++) { - if (i == 0) - { - nnz = 1; - cols[0] = 0; - values[0] = 1; - } - else - { - nnz = 2; - cols[0] = i - 1; - cols[1] = i; - values[0] = -1; - values[1] = 1; - } + if (i == 0) + { + nnz = 1; + cols[0] = 0; + values[0] = 1; + } + else + { + nnz = 2; + cols[0] = i - 1; + cols[1] = i; + values[0] = -1; + values[1] = 1; + } /* Set the values for row i */ tmp[0] = nnz; @@ -270,7 +272,6 @@ int main (int argc, char *argv[]) HYPRE_ParVectorSetConstantValues(par_b, 0); HYPRE_ParCSRMatrixMatvec(1, parcsr_A, par_work, 0, par_b); - /* Print out the system - files names will be IJ.out.A.XXXXX and IJ.out.b.XXXXX, where XXXXX = processor id */ if (print_system) @@ -280,7 +281,8 @@ int main (int argc, char *argv[]) HYPRE_IJVectorPrint(b, "IJ.out.b"); } - /* this is a 2-D 4-by-k array using Double pointers */ + /* This is a 2-D 4-by-k array. Hypre takes ownership of this pointer, + so users do not need to free it */ grid_relax_points = hypre_CTAlloc(int*, 4, HYPRE_MEMORY_HOST); grid_relax_points[0] = NULL; grid_relax_points[1] = hypre_CTAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST); @@ -292,7 +294,7 @@ int main (int argc, char *argv[]) grid_relax_points[1][i] = 1; // C } /* up cycle: F */ - for (i=0; i Date: Mon, 26 Jun 2023 00:01:29 -0400 Subject: [PATCH 7/7] Replace hypre allocs with malloc --- src/examples/ex19.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/examples/ex19.c b/src/examples/ex19.c index 91efe99b84..2c100c242d 100644 --- a/src/examples/ex19.c +++ b/src/examples/ex19.c @@ -283,11 +283,11 @@ int main (int argc, char *argv[]) /* This is a 2-D 4-by-k array. Hypre takes ownership of this pointer, so users do not need to free it */ - grid_relax_points = hypre_CTAlloc(int*, 4, HYPRE_MEMORY_HOST); + grid_relax_points = (int**) malloc(4 * sizeof(int*)); grid_relax_points[0] = NULL; - grid_relax_points[1] = hypre_CTAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST); - grid_relax_points[2] = hypre_CTAlloc(HYPRE_Int, ns_up, HYPRE_MEMORY_HOST); - grid_relax_points[3] = hypre_CTAlloc(HYPRE_Int, ns_coarse, HYPRE_MEMORY_HOST); + grid_relax_points[1] = (int*) malloc(ns_down * sizeof(int)); + grid_relax_points[2] = (int*) malloc(ns_up * sizeof(int)); + grid_relax_points[3] = (int*) malloc(ns_coarse * sizeof(int)); /* down cycle: C */ for (i = 0; i < ns_down; i++) {