Add HYPRE interop

CSparse.Benchmark/Examples/Double/TestHypre.cs

@@ -0,0 +1,157 @@
+
+namespace CSparse.Double.Examples
+{
+    using System;
+    using CSparse.Interop.Hypre;
+
+    static class TestHypre
+    {
+        public static void Run()
+        {
+            NativeMethods.HYPRE_Init();
+
+            int n = 33;
+            var h = 1.0 / (n + 1);
+
+            var A = CSparse.Double.Generate.Laplacian(n, n);
+            var b = Vector.Create(A.RowCount, h * h);
+
+            using (var hypreA = new HypreMatrix<double>(A, true))
+            {
+                Check(TestAMG(hypreA, b), "AMG");
+                Check(TestPCG(hypreA, b), "PCG");
+                Check(TestPCG_AMG(hypreA, b), "PCG + AMG");
+                Check(TestPCG_ParaSails(hypreA, b), "PCG + ParaSails");
+                Check(TestFlexGMRES_AMG(hypreA, b), "FlexGMRES + AMG");
+            }
+
+            NativeMethods.HYPRE_Finalize();
+        }
+
+        private static void Check(HypreResult result, string name)
+        {
+            Console.WriteLine(name);
+            Console.WriteLine("  Iterations = {0}", result.NumIterations);
+            Console.WriteLine("  Final Relative Residual Norm = {0}", result.RelResidualNorm);
+            Console.WriteLine();
+        }
+
+        private static HypreResult TestAMG(HypreMatrix<double> A, double[] _b)
+        {
+            var _x = Vector.Create(A.RowCount, 0.0);
+
+            using var x = new HypreVector<double>(_x);
+            using var b = new HypreVector<double>(_b);
+
+            using var solver = new BoomerAMG<double>();
+
+            solver.SetOldDefault(); // Falgout coarsening with modified classical interpolation
+            solver.SetRelaxType(3); // G-S/Jacobi hybrid relaxation
+            solver.SetRelaxOrder(1); // uses C/F relaxation
+            solver.SetNumSweeps(1);  // Sweeps on each level
+            solver.SetMaxLevels(20);
+            solver.SetTol(1e-7);
+            //solver.PrintLevel = 3;
+
+            return solver.Solve(A, x, b);
+        }
+
+        private static HypreResult TestPCG(HypreMatrix<double> A, double[] _b)
+        {
+            var _x = Vector.Create(A.RowCount, 0.0);
+
+            using var x = new HypreVector<double>(_x);
+            using var b = new HypreVector<double>(_b);
+
+            using var solver = new PCG<double>();
+
+            solver.SetMaxIter(1000);
+            solver.SetTol(1e-7);
+            solver.SetTwoNorm(1); // use the two norm as the stopping criteria
+            //solver.PrintLevel = 2;
+            //solver.Logging = 1;
+
+            return solver.Solve(A, x, b);
+        }
+
+        private static HypreResult TestPCG_AMG(HypreMatrix<double> A, double[] _b)
+        {
+            var _x = Vector.Create(A.RowCount, 0.0);
+
+            using var x = new HypreVector<double>(_x);
+            using var b = new HypreVector<double>(_b);
+
+            using var precond = new BoomerAMG<double>();
+
+            //precond.PrintLevel = 1;
+            precond.SetCoarsenType(6);
+            precond.SetOldDefault();
+            precond.SetRelaxType(6); // Sym G.S./Jacobi hybrid
+            precond.SetNumSweeps(1);
+            precond.SetTol(0.0);
+            precond.SetMaxIter(1); // do only one iteration!
+
+            using var solver = new PCG<double>(precond);
+
+            solver.SetMaxIter(1000);
+            solver.SetTol(1e-7);
+            solver.SetTwoNorm(1); // use the two norm as the stopping criteria
+            //solver.PrintLevel = 2;
+            //solver.Logging = 1;
+
+            return solver.Solve(A, x, b);
+        }
+
+        private static HypreResult TestPCG_ParaSails(HypreMatrix<double> A, double[] _b)
+        {
+            var _x = Vector.Create(A.RowCount, 0.0);
+
+            using var x = new HypreVector<double>(_x);
+            using var b = new HypreVector<double>(_b);
+
+            using var precond = new ParaSails<double>(true);
+
+            precond.SetParams(0.1, 1);
+            precond.SetFilter(0.05);
+            //precond.Logging = 3;
+
+            using var solver = new PCG<double>(precond);
+
+            solver.SetMaxIter(1000);
+            solver.SetTol(1e-7);
+            solver.SetTwoNorm(1); // use the two norm as the stopping criteria
+            //solver.PrintLevel = 2;
+            //solver.Logging = 1;
+
+            return solver.Solve(A, x, b);
+        }
+
+        private static HypreResult TestFlexGMRES_AMG(HypreMatrix<double> A, double[] _b)
+        {
+            var _x = Vector.Create(A.RowCount, 0.0);
+
+            using var x = new HypreVector<double>(_x);
+            using var b = new HypreVector<double>(_b);
+
+            using var precond = new BoomerAMG<double>();
+
+            //precond.PrintLevel = 1;
+            precond.SetCoarsenType(6);
+            precond.SetOldDefault();
+            precond.SetRelaxType(6); // Sym G.S./Jacobi hybrid
+            precond.SetNumSweeps(1);
+            precond.SetTol(0.0); // conv. tolerance zero
+            precond.SetMaxIter(1); // do only one iteration!
+
+            using var solver = new FlexGMRES<double>(precond);
+
+            solver.SetKDim(30);
+            solver.SetMaxIter(1000);
+            solver.SetTol(1e-7);
+            //solver.PrintLevel = 2;
+            //solver.Logging = 1;
+
+            return solver.Solve(A, x, b);
+        }
+    }
+}

CSparse.Interop/CSparse.Interop.csproj

@@ -30,5 +30,9 @@
   <ItemGroup>
     <PackageReference Include="CSparse" Version="4.2.0" />
   </ItemGroup>
+
+  <ItemGroup>
+    <Folder Include="Double\Hypre\" />
+  </ItemGroup>
 
 </Project>

CSparse.Interop/Interop/Hypre/BiCGSTAB.cs

@@ -0,0 +1,68 @@
+using System;
+
+namespace CSparse.Interop.Hypre
+{
+    public class BiCGSTAB<T> : HypreSolver<T>
+        where T : struct, IEquatable<T>, IFormattable
+    {
+        public BiCGSTAB()
+        {
+            NativeMethods.HYPRE_ParCSRBiCGSTABCreate(Constants.MPI_COMM_WORLD, out solver);
+        }
+
+        public BiCGSTAB(IHyprePreconditioner<T> precond)
+            : this()
+        {
+            precond.Bind(this);
+        }
+
+        internal override void SetPreconditioner(
+            HYPRE_PtrToParSolverFcn precond,
+            HYPRE_PtrToParSolverFcn precond_setup,
+            HYPRE_Solver precond_solver)
+        {
+            NativeMethods.HYPRE_ParCSRBiCGSTABSetPrecond(solver, precond, precond_setup, precond_solver);
+        }
+
+        public override HypreResult Solve(HypreMatrix<T> A, HypreVector<T> x, HypreVector<T> b)
+        {
+            NativeMethods.HYPRE_BiCGSTABSetPrintLevel(solver, PrintLevel);
+            NativeMethods.HYPRE_BiCGSTABSetLogging(solver, Logging);
+
+            NativeMethods.HYPRE_BiCGSTABSetTol(solver, 1e-7);
+            //-HYPRE_BiCGSTABSetAbsoluteTol
+            //-HYPRE_BiCGSTABSetConvergenceFactorTol
+            //-HYPRE_BiCGSTABSetStopCrit
+            //-HYPRE_BiCGSTABSetMinIter
+            NativeMethods.HYPRE_BiCGSTABSetMaxIter(solver, 1000);
+
+            var parcsr_A = A.GetObject();
+            var par_b = b.GetObject();
+            var par_x = x.GetObject();
+
+            NativeMethods.HYPRE_ParCSRBiCGSTABSetup(solver, parcsr_A, par_b, par_x);
+            NativeMethods.HYPRE_ParCSRBiCGSTABSolve(solver, parcsr_A, par_b, par_x);
+
+            x.Synchronize();
+
+            HypreResult result;
+
+            NativeMethods.HYPRE_BiCGSTABGetNumIterations(solver, out result.NumIterations);
+            //NativeMethods.HYPRE_BiCGSTABGetConverged(solver, out result.Converged);
+            NativeMethods.HYPRE_BiCGSTABGetFinalRelativeResidualNorm(solver, out result.RelResidualNorm);
+
+            result.Converged = 1; // result.NumIterations < max_iter;
+
+            return result;
+        }
+
+        protected override void Dispose(bool disposing)
+        {
+            if (solver.Ptr != IntPtr.Zero)
+            {
+                NativeMethods.HYPRE_ParCSRBiCGSTABDestroy(solver);
+                solver.Ptr = IntPtr.Zero;
+            }
+        }
+    }
+}

CSparse.Interop/Interop/Hypre/BoomerAMG.cs

@@ -0,0 +1,93 @@
+
+namespace CSparse.Interop.Hypre
+{
+    using System;
+
+    public class BoomerAMG<T> : HypreSolver<T>, IHyprePreconditioner<T>
+        where T : struct, IEquatable<T>, IFormattable
+    {
+        public BoomerAMG()
+        {
+            NativeMethods.HYPRE_BoomerAMGCreate(out solver);
+        }
+
+        public void Bind(HypreSolver<T> solver)
+        {
+            solver.SetPreconditioner(NativeMethods.HYPRE_BoomerAMGSolve, NativeMethods.HYPRE_BoomerAMGSetup, this.solver);
+        }
+
+        public void SetTol(double tol)
+        {
+            NativeMethods.HYPRE_BoomerAMGSetTol(solver, tol);
+        }
+
+        public void SetMaxLevels(int max_levels)
+        {
+            NativeMethods.HYPRE_BoomerAMGSetMaxLevels(solver, max_levels);
+        }
+
+        public void SetNumSweeps(int num_sweeps)
+        {
+            NativeMethods.HYPRE_BoomerAMGSetNumSweeps(solver, num_sweeps);
+        }
+
+        public void SetRelaxOrder(int relax_order)
+        {
+            NativeMethods.HYPRE_BoomerAMGSetRelaxOrder(solver, relax_order);
+        }
+
+        public void SetRelaxType(int relax_type)
+        {
+            NativeMethods.HYPRE_BoomerAMGSetRelaxType(solver, relax_type);
+        }
+
+        public void SetCoarsenType(int coarsen)
+        {
+            NativeMethods.HYPRE_BoomerAMGSetCoarsenType(solver, coarsen);
+        }
+
+        public void SetOldDefault()
+        {
+            NativeMethods.HYPRE_BoomerAMGSetOldDefault(solver);
+        }
+
+        public void SetMaxIter(int max_iter)
+        {
+            NativeMethods.HYPRE_BoomerAMGSetMaxIter(solver, max_iter);
+        }
+
+        public override HypreResult Solve(HypreMatrix<T> A, HypreVector<T> x, HypreVector<T> b)
+        {
+            NativeMethods.HYPRE_BoomerAMGSetPrintLevel(solver, PrintLevel);
+            NativeMethods.HYPRE_BoomerAMGSetLogging(solver, Logging);
+
+            var par_A = A.GetObject();
+            var par_b = b.GetObject();
+            var par_x = x.GetObject();
+
+            NativeMethods.HYPRE_BoomerAMGSetup(solver, par_A, par_b, par_x);
+            NativeMethods.HYPRE_BoomerAMGSolve(solver, par_A, par_b, par_x);
+
+            x.Synchronize();
+
+            HypreResult result;
+
+            NativeMethods.HYPRE_BoomerAMGGetNumIterations(solver, out result.NumIterations);
+            //NativeMethods.HYPRE_BoomerAMGGetConverged(solver, out result.Converged);
+            NativeMethods.HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, out result.RelResidualNorm);
+
+            result.Converged = 1;
+
+            return result;
+        }
+
+        protected override void Dispose(bool disposing)
+        {
+            if (solver.Ptr != IntPtr.Zero)
+            {
+                NativeMethods.HYPRE_BoomerAMGDestroy(solver);
+                solver.Ptr = IntPtr.Zero;
+            }
+        }
+    }
+}