From 450162d178c8f1ccb1178e3482a3858d3dc3e3f2 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Wed, 15 Jun 2022 18:35:31 +0200 Subject: [PATCH] Add translation of example 5 from hypre source code. --- examples/ex5.jl | 258 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 258 insertions(+) create mode 100644 examples/ex5.jl diff --git a/examples/ex5.jl b/examples/ex5.jl new file mode 100644 index 0000000..deeaff8 --- /dev/null +++ b/examples/ex5.jl @@ -0,0 +1,258 @@ +# Example translated from C to Julia based on this original source (MIT license): +# https://github.com/hypre-space/hypre/blob/ac9d7d0d7b43cd3d0c7f24ec5d64b58fbf900097/src/examples/ex5.c + +# Example 5 +# +# Interface: Linear-Algebraic (IJ) +# +# Sample run: mpirun -np 4 julia ex5.jl +# +# Description: This example solves the 2-D Laplacian problem with zero boundary +# conditions on an n x n grid. The number of unknowns is N=n^2. +# The standard 5-point stencil is used, and we solve for the +# interior nodes only. +# +# This example solves the same problem as Example 3. Available +# solvers are AMG, PCG, and PCG with AMG or Parasails +# preconditioners. + +using MPI +using HYPRE.LibHYPRE # LibHYPRE submodule contains the raw C-interface generated by Clang.jl + +function main() + + # Initialize MPI + # MPI_Init(Ref{Cint}(length(ARGS)), ARGS) + MPI.Init() + comm = MPI.COMM_WORLD + myid = MPI.Comm_rank(comm) + num_procs = MPI.Comm_size(comm) + + # Initialize HYPRE + HYPRE_Init() + + # Default problem parameters + n = 33 + solver_id = 0 + print_system = false + + # TODO: Parse command line + + # Preliminaries: want at least one processor per row + if n * n < num_procs; n = trunc(Int, sqrt(n)) + 1; end + N = n * n # global number of rows + h = 1.0 / (n + 1) # mesh size + h2 = h * h + + # 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 += min(myid, extra) + + iupper = local_size * (myid + 1) + iupper += 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) + Aref = Ref{HYPRE_IJMatrix}(C_NULL) + HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, Aref) + A = Aref[] + + # Choose a parallel csr format storage (see the User's Manual) + HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR) + + # Initialize before setting coefficients + HYPRE_IJMatrixInitialize(A) + + # Now go through my local rows and set the matrix entries. + # Each row has at most 5 entries. For example, if n=3: + + # A = [M -I 0; -I M -I; 0 -I M] + # M = [4 -1 0; -1 4 -1; 0 -1 4] + + # Note that here we are setting one row at a time, though + # one could set all the rows together (see the User's Manual). + let + values = Vector{Float64}(undef, 5) + cols = Vector{HYPRE_Int}(undef, 5) + tmp = Vector{HYPRE_Int}(undef, 2) + + for i in ilower:iupper + nnz = 1 + + # The left identity block:position i-n + if (i - n) >= 0 + cols[nnz] = i - n + values[nnz] = -1.0 + nnz += 1 + end + + # The left -1: position i-1 + if i % n != 0 + cols[nnz] = i - 1 + values[nnz] = -1.0 + nnz += 1 + end + + # Set the diagonal: position i + cols[nnz] = i + values[nnz] = 4.0 + nnz += 1 + + # The right -1: position i+1 + if (i + 1) % n != 0 + cols[nnz] = i + 1 + values[nnz] = -1.0 + nnz += 1 + end + + # The right identity block:position i+n + if (i + n) < N + cols[nnz] = i + n + values[nnz] = -1.0 + nnz += 1 + end + + # Set the values for row i + tmp[1] = nnz - 1 + tmp[2] = i + + HYPRE_IJMatrixSetValues(A, 1, Ref(tmp[1]), Ref(tmp[2]), cols, values) + end + end + + # 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 + parcsr_A_ref = Ref{Ptr{Cvoid}}(C_NULL) + HYPRE_IJMatrixGetObject(A, parcsr_A_ref) + parcsr_A = convert(Ptr{HYPRE_ParCSRMatrix}, parcsr_A_ref[]) + + # Create the rhs and solution + b_ref = Ref{HYPRE_IJVector}(C_NULL) + HYPRE_IJVectorCreate(comm, ilower, iupper, b_ref) + b = b_ref[] + HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR) + HYPRE_IJVectorInitialize(b) + + x_ref = Ref{HYPRE_IJVector}(C_NULL) + HYPRE_IJVectorCreate(comm, ilower, iupper, x_ref) + x = x_ref[] + HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR) + HYPRE_IJVectorInitialize(x) + + # Set the rhs values to h^2 and the solution to zero + let + rhs_values = Vector{Float64}(undef, local_size) + x_values = Vector{Float64}(undef, local_size) + rows = Vector{HYPRE_Int}(undef, local_size) + + for i in 1:local_size + rhs_values[i] = h2 + x_values[i] = 0.0 + rows[i] = ilower + i - 1 + end + + HYPRE_IJVectorSetValues(b, local_size, rows, rhs_values) + HYPRE_IJVectorSetValues(x, local_size, rows, x_values) + end + + HYPRE_IJVectorAssemble(b) + par_b_ref = Ref{Ptr{Cvoid}}(C_NULL) + HYPRE_IJVectorGetObject(b, par_b_ref) + par_b = convert(Ptr{HYPRE_ParVector}, par_b_ref[]) + + HYPRE_IJVectorAssemble(x) + par_x_ref = Ref{Ptr{Cvoid}}(C_NULL) + HYPRE_IJVectorGetObject(x, par_x_ref) + par_x = convert(Ptr{HYPRE_ParVector}, par_x_ref[]) + + # 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(b, "IJ.out.b") + end + + # Choose a solver and solve the system + solver_ref = Ref{HYPRE_Solver}(C_NULL) + + # AMG + if solver_id == 0 + # Create solver + HYPRE_BoomerAMGCreate(solver_ref) + solver = solver_ref[] + + # 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, 3) # G-S/Jacobi hybrid relaxation + HYPRE_BoomerAMGSetRelaxOrder(solver, 1) # uses C/F relaxation + HYPRE_BoomerAMGSetNumSweeps(solver, 1) # Sweeeps on each level + HYPRE_BoomerAMGSetMaxLevels(solver, 20) # maximum number of levels + HYPRE_BoomerAMGSetTol(solver, 1e-7) # conv. tolerance + + # Now setup and solve! + HYPRE_BoomerAMGSetup(solver, parcsr_A, par_b, par_x) + HYPRE_BoomerAMGSolve(solver, parcsr_A, par_b, par_x) + + # Run info - needed logging turned on + num_iterations = Ref{HYPRE_Int}() + final_res_norm = Ref{Float64}() + HYPRE_BoomerAMGGetNumIterations(solver, num_iterations) + HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, final_res_norm) + if myid == 0 + println() + println("Iterations = $(num_iterations[])") + println("Final Relative Residual Norm = $(final_res_norm[])") + println() + end + + # Destroy solver + HYPRE_BoomerAMGDestroy(solver) + else + if myid == 0 + println("Invalid sovler id specified.") + end + end + + # Clean up + HYPRE_IJMatrixDestroy(A) + HYPRE_IJVectorDestroy(b) + HYPRE_IJVectorDestroy(x) + + # Finalize HYPRE + HYPRE_Finalize() + + # Finalize MPI + MPI.Finalize() + + return 0 +end + +# Run it +if abspath(PROGRAM_FILE) == @__FILE__ + main() +end