Finish translation of example 5.

examples/ex5.jl

@@ -19,14 +19,14 @@
 using MPI
 using HYPRE.LibHYPRE # LibHYPRE submodule contains the raw C-interface generated by Clang.jl
 
-function main()
+function main(argc, argv)
 
     # Initialize MPI
     # MPI_Init(Ref{Cint}(length(ARGS)), ARGS)
     MPI.Init()
-    comm = MPI.COMM_WORLD
+    MPI_COMM_WORLD = MPI.COMM_WORLD
-    myid = MPI.Comm_rank(comm)
+    myid = MPI.Comm_rank(MPI_COMM_WORLD)
-    num_procs = MPI.Comm_size(comm)
+    num_procs = MPI.Comm_size(MPI_COMM_WORLD)
 
     # Initialize HYPRE
     HYPRE_Init()
@@ -36,7 +36,49 @@ function main()
     solver_id = 0
     print_system = false
 
-    # TODO: Parse command line
+    # Parse command line
+    arg_index = 1
+    print_usage = false
+    let
+        while arg_index <= argc
+            if argv[arg_index] == "-n"
+                n = parse(Int, argv[arg_index += 1])
+                arg_index += 1
+            elseif argv[arg_index] == "-solver"
+                solver_id = parse(Int, argv[arg_index += 1])
+                arg_index += 1
+            elseif argv[arg_index] == "-print_system"
+                print_system = true
+                arg_index += 1
+            elseif argv[arg_index] == "-help"
+                print_usage = true
+                break
+            else
+                arg_index += 1
+            end
+        end
+
+        if print_usage && myid == 0
+            println()
+            println("Usage: julia ex5.jl [<options>]")
+            println()
+            println("  -n <n>              : problem size in each direction (default: 33)")
+            println("  -solver <ID>        : solver ID")
+            println("                        0  - AMG (default)")
+            println("                        1  - AMG-PCG")
+            println("                        8  - ParaSails-PCG")
+            println("                        50 - PCG")
+            println("                        61 - AMG-FlexGMRES")
+            # println("  -vis                : save the solution for GLVis visualization")
+            println("  -print_system       : print the matrix and rhs")
+            println()
+        end
+
+        if print_usage
+            MPI.Finalize()
+            return 0
+        end
+    end
 
     # Preliminaries: want at least one processor per row
     if n * n < num_procs; n = trunc(Int, sqrt(n)) + 1; end
@@ -64,7 +106,7 @@ function main()
     # 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)
+    HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, Aref)
     A = Aref[]
 
     # Choose a parallel csr format storage (see the User's Manual)
@@ -82,9 +124,9 @@ function main()
     # 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)
+        values = Vector{Cdouble}(undef, 5)
-        cols = Vector{HYPRE_Int}(undef, 5)
+        cols = Vector{Cint}(undef, 5)
-        tmp = Vector{HYPRE_Int}(undef, 2)
+        tmp = Vector{Cint}(undef, 2)
 
         for i in ilower:iupper
             nnz = 1
@@ -150,24 +192,25 @@ function main()
     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)
+    HYPRE_IJVectorCreate(MPI_COMM_WORLD, 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)
+    HYPRE_IJVectorCreate(MPI_COMM_WORLD, 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)
+        rhs_values = Vector{Cdouble}(undef, local_size)
-        x_values = Vector{Float64}(undef, local_size)
+        x_values = Vector{Cdouble}(undef, local_size)
-        rows = Vector{HYPRE_Int}(undef, local_size)
+        rows = Vector{Cint}(undef, local_size)
 
         for i in 1:local_size
             rhs_values[i] = h2
@@ -179,6 +222,13 @@ function main()
         HYPRE_IJVectorSetValues(x, local_size, rows, x_values)
     end
 
+    # As with the matrix, for testing purposes, one may wish to read in a rhs:
+    # HYPRE_IJVectorRead( <filename>, MPI_COMM_WORLD,
+    #                           HYPRE_PARCSR, &b );
+    # as an alternative to the
+    # following sequence of HYPRE_IJVectors calls:
+    # Create, SetObjectType, Initialize, SetValues, and Assemble
+
     HYPRE_IJVectorAssemble(b)
     par_b_ref = Ref{Ptr{Cvoid}}(C_NULL)
     HYPRE_IJVectorGetObject(b, par_b_ref)
@@ -198,6 +248,10 @@ function main()
 
     # Choose a solver and solve the system
     solver_ref = Ref{HYPRE_Solver}(C_NULL)
+    precond_ref = Ref{HYPRE_Solver}(C_NULL)
+
+    num_iterations = Ref{Cint}()
+    final_res_norm = Ref{Cdouble}()
 
     # AMG
     if solver_id == 0
@@ -219,8 +273,6 @@ function main()
         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
@@ -232,10 +284,178 @@ function main()
 
         # Destroy solver
         HYPRE_BoomerAMGDestroy(solver)
-    else
+
+    # PCG
+    elseif solver_id == 50
+        # Create solver
+        HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, solver_ref)
+        solver = solver_ref[]
+
+        # Set some parameters (See Reference Manual for more parameters)
+        HYPRE_PCGSetMaxIter(solver, 1000) # max iterations
+        HYPRE_PCGSetTol(solver, 1e-7) # conv. tolerance
+        HYPRE_PCGSetTwoNorm(solver, 1) # use the two norm as the stopping criteria
+        HYPRE_PCGSetPrintLevel(solver, 2) # prints out the iteration info
+        HYPRE_PCGSetLogging(solver, 1) # needed to get run info later
+
+        # Now setup and solve!
+        HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x)
+        HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x)
+
+        # Run info - needed logging turned on
+        HYPRE_PCGGetNumIterations(solver, num_iterations)
+        HYPRE_PCGGetFinalRelativeResidualNorm(solver, final_res_norm)
         if myid == 0
-            println("Invalid sovler id specified.")
+            println()
+            println("Iterations = $(num_iterations[])")
+            println("Final Relative Residual Norm = $(final_res_norm[])")
+            println()
         end
+
+        # Destroy solver
+        HYPRE_ParCSRPCGDestroy(solver)
+
+    # PCG with AMG preconditioner
+    elseif solver_id == 1
+        # Create solver
+        HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, solver_ref)
+        solver = solver_ref[]
+
+        # Set some parameters (See Reference Manual for more parameters)
+        HYPRE_PCGSetMaxIter(solver, 1000) # max iterations
+        HYPRE_PCGSetTol(solver, 1e-7) # conv. tolerance
+        HYPRE_PCGSetTwoNorm(solver, 1) # use the two norm as the stopping criteria
+        HYPRE_PCGSetPrintLevel(solver, 2) # print solve info
+        HYPRE_PCGSetLogging(solver, 1) # needed to get run info later
+
+        # Now set up the AMG preconditioner and specify any parameters
+        HYPRE_BoomerAMGCreate(precond_ref)
+        precond = precond_ref[]
+        HYPRE_BoomerAMGSetPrintLevel(precond, 1) # print amg solution info
+        HYPRE_BoomerAMGSetCoarsenType(precond, 6)
+        HYPRE_BoomerAMGSetOldDefault(precond)
+        HYPRE_BoomerAMGSetRelaxType(precond, 6) # Sym G.S./Jacobi hybrid
+        HYPRE_BoomerAMGSetNumSweeps(precond, 1)
+        HYPRE_BoomerAMGSetTol(precond, 0.0) # conv. tolerance zero
+        HYPRE_BoomerAMGSetMaxIter(precond, 1) # do only one iteration!
+
+        # Set the PCG preconditioner
+        HYPRE_PCGSetPrecond(solver, HYPRE_BoomerAMGSolve, HYPRE_BoomerAMGSetup, precond)
+
+        # Now setup and solve!
+        HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x)
+        HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x)
+
+        # Run info - needed logging turned on
+        HYPRE_PCGGetNumIterations(solver, num_iterations)
+        HYPRE_PCGGetFinalRelativeResidualNorm(solver, final_res_norm)
+        if myid == 0
+            println()
+            println("Iterations = $(num_iterations[])")
+            println("Final Relative Residual Norm = $(final_res_norm[])")
+            println()
+        end
+
+        # Destroy solver and preconditioner
+        HYPRE_ParCSRPCGDestroy(solver)
+        HYPRE_BoomerAMGDestroy(precond)
+
+    # PCG with Parasails Preconditioner
+    elseif solver_id == 8
+        # Create solver
+        HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, solver_ref)
+        solver = solver_ref[]
+
+        # Set some parameters (See Reference Manual for more parameters)
+        HYPRE_PCGSetMaxIter(solver, 1000) # max iterations
+        HYPRE_PCGSetTol(solver, 1e-7) # conv. tolerance
+        HYPRE_PCGSetTwoNorm(solver, 1) # use the two norm as the stopping criteria
+        HYPRE_PCGSetPrintLevel(solver, 2) # print solve info
+        HYPRE_PCGSetLogging(solver, 1) # needed to get run info later
+
+        # Now set up the ParaSails preconditioner and specify any parameters
+        HYPRE_ParaSailsCreate(MPI_COMM_WORLD, precond_ref)
+        precond = precond_ref[]
+
+        # Set some parameters (See Reference Manual for more parameters)
+        sai_max_levels = 1
+        sai_threshold = 0.1
+        sai_filter = 0.05
+        sai_sym = 1
+        HYPRE_ParaSailsSetParams(precond, sai_threshold, sai_max_levels)
+        HYPRE_ParaSailsSetFilter(precond, sai_filter)
+        HYPRE_ParaSailsSetSym(precond, sai_sym)
+        HYPRE_ParaSailsSetLogging(precond, 3)
+
+        # Set the PCG preconditioner
+        HYPRE_PCGSetPrecond(solver, HYPRE_ParaSailsSolve, HYPRE_ParaSailsSetup, precond)
+
+        # Now setup and solve!
+        HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x)
+        HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x)
+
+        # Run info - needed logging turned on
+        HYPRE_PCGGetNumIterations(solver, num_iterations)
+        HYPRE_PCGGetFinalRelativeResidualNorm(solver, final_res_norm)
+        if myid == 0
+            println()
+            println("Iterations = $(num_iterations[])")
+            println("Final Relative Residual Norm = $(final_res_norm[])")
+            println()
+        end
+
+        # Destory solver and preconditioner
+        HYPRE_ParCSRPCGDestroy(solver)
+        HYPRE_ParaSailsDestroy(precond)
+
+    # Flexible GMRES with AMG Preconditioner
+    elseif solver_id == 61
+
+        # Create solver
+        HYPRE_ParCSRFlexGMRESCreate(MPI_COMM_WORLD, solver_ref)
+        solver = solver_ref[]
+
+        # Set some parameters (See Reference Manual for more parameters)
+        HYPRE_FlexGMRESSetKDim(solver, 30) # restart
+        HYPRE_FlexGMRESSetMaxIter(solver, 1000) # max iterations
+        HYPRE_FlexGMRESSetTol(solver, 1e-7) # conv. tolerance
+        HYPRE_FlexGMRESSetPrintLevel(solver, 2) # print solve info
+        HYPRE_FlexGMRESSetLogging(solver, 1) # needed to get run info later
+
+        # Now set up the AMG preconditioner and specify any parameters
+        HYPRE_BoomerAMGCreate(precond_ref)
+        precond = precond_ref[]
+        HYPRE_BoomerAMGSetPrintLevel(precond, 1) # print amg solution info
+        HYPRE_BoomerAMGSetCoarsenType(precond, 6)
+        HYPRE_BoomerAMGSetOldDefault(precond)
+        HYPRE_BoomerAMGSetRelaxType(precond, 6) # Sym G.S./Jacobi hybrid
+        HYPRE_BoomerAMGSetNumSweeps(precond, 1)
+        HYPRE_BoomerAMGSetTol(precond, 0.0) # conv. tolerance zero
+        HYPRE_BoomerAMGSetMaxIter(precond, 1) # do only one iteration!
+
+        # Set the FlexGMRES preconditioner
+        HYPRE_FlexGMRESSetPrecond(solver, HYPRE_BoomerAMGSolve, HYPRE_BoomerAMGSetup, precond)
+
+        # Now setup and solve!
+        HYPRE_ParCSRFlexGMRESSetup(solver, parcsr_A, par_b, par_x)
+        HYPRE_ParCSRFlexGMRESSolve(solver, parcsr_A, par_b, par_x)
+
+        # Run info - needed logging turned on
+        HYPRE_FlexGMRESGetNumIterations(solver, num_iterations)
+        HYPRE_FlexGMRESGetFinalRelativeResidualNorm(solver, final_res_norm)
+        if myid == 0
+            println()
+            println("Iterations = $(num_iterations[])")
+            println("Final Relative Residual Norm = $(final_res_norm[])")
+            println()
+        end
+
+        # Destory solver and preconditioner
+        HYPRE_ParCSRFlexGMRESDestroy(solver)
+        HYPRE_BoomerAMGDestroy(precond)
+
+    else
+        if myid == 0; println("Invalid solver id specified."); end
     end
 
     # Clean up
@@ -254,5 +474,5 @@ end
 
 # Run it
 if abspath(PROGRAM_FILE) == @__FILE__
-    main()
+    main(length(ARGS), ARGS)
 end