Merge pull request #236 from PFLAREProject/ilu_factors

Ilu factors

Author: stevendargaville

Makefile

@@ -162,7 +162,8 @@ export TEST_TARGETS = ex12f \
 		  ex12f_gmres_poly \
 		  mat_diag \
 		  adv_dg_upwind \
-		  ex6_two_airg
+		  ex6_two_airg \
+		  ilu_factors
 # Include kokkos examples
 ifeq ($(PETSC_HAVE_KOKKOS),1)
 export TEST_TARGETS := $(TEST_TARGETS) adv_1dk

include/pflare.h

@@ -45,6 +45,13 @@ PETSC_EXTERN void PCRegister_PFLARE();
 PETSC_EXTERN void compute_cf_splitting(Mat, int, double, int, int, int, double, IS*, IS*);
 PETSC_EXTERN void compute_diag_dom_submatrix(Mat, double, Mat*);
 
+/* Restrict input_mat onto output_mat's existing sparsity pattern.
+   With alpha_int = 1, performs output_mat += alpha * input_mat on output_mat's
+   pattern; with alpha_int = 0, performs output_mat = input_mat on output's pattern.
+   lump_int = 1 adds dropped entries to the diagonal. Auto-dispatches between the
+   CPU and Kokkos implementations based on the input matrix type. */
+PETSC_EXTERN void remove_from_sparse_match(Mat, Mat, int, int, PetscReal);
+
 /* Define PCPFLAREINV get routines */
 PETSC_EXTERN PetscErrorCode PCPFLAREINVGetPolyOrder(PC, PetscInt *);
 PETSC_EXTERN PetscErrorCode PCPFLAREINVGetSparsityOrder(PC, PetscInt *);

src/C_Fortran_Bindings.F90

@@ -7,6 +7,7 @@ module c_fortran_bindings
    use approx_inverse_setup, only: calculate_and_build_approximate_inverse, reset_inverse_mat
    use cf_splitting, only: compute_cf_splitting
    use matdiagdomsubmatrix, only: compute_diag_dom_submatrix
+   use petsc_helper, only: remove_from_sparse_match
    use air_data_type_routines, only: create_air_data
 
 #include "petsc/finclude/petscksp.h"
@@ -281,5 +282,43 @@ end subroutine compute_diag_dom_submatrix_c
 
    !------------------------------------------------------------------------------------------------------------------------
 
+   subroutine remove_from_sparse_match_c(input_mat_ptr, output_mat_ptr, lump_int, alpha_int, alpha) &
+         bind(C, name='remove_from_sparse_match_c')
+
+      ! Restrict input_mat onto output_mat's sparsity pattern. The underlying
+      ! Fortran remove_from_sparse_match auto-dispatches between the CPU and
+      ! Kokkos implementations based on the matrix type, so this C entry point
+      ! works for both paths.
+
+      ! ~~~~~~~~
+      integer(c_long_long), intent(in)    :: input_mat_ptr
+      integer(c_long_long), intent(inout) :: output_mat_ptr
+      integer(c_int), value, intent(in)   :: lump_int
+      integer(c_int), value, intent(in)   :: alpha_int
+      real(c_double), value, intent(in)   :: alpha
+
+      type(tMat) :: input_mat, output_mat
+      logical    :: lump
+      ! ~~~~~~~~
+
+      input_mat%v  = input_mat_ptr
+      output_mat%v = output_mat_ptr
+      lump         = (lump_int /= 0)
+
+      if (alpha_int /= 0) then
+         call remove_from_sparse_match(input_mat, output_mat, lump, alpha)
+      else
+         call remove_from_sparse_match(input_mat, output_mat, lump)
+      end if
+
+      ! The matrix's identity isn't replaced by remove_from_sparse_match (only
+      ! its values), but pass the handle back through anyway to match the
+      ! pattern used by the other C wrappers in this module.
+      output_mat_ptr = output_mat%v
+
+   end subroutine remove_from_sparse_match_c
+
+   !------------------------------------------------------------------------------------------------------------------------
+
 end module c_fortran_bindings
 

src/PCAIR.c

@@ -18,6 +18,8 @@ PETSC_EXTERN void compute_cf_splitting_c(Mat *input_mat, int symmetric_int,
    int ddc_its, double fraction_swap,
    IS *is_fine, IS *is_coarse);
 PETSC_EXTERN void compute_diag_dom_submatrix_c(Mat *input_mat, double max_dd_ratio, Mat *output_mat);
+PETSC_EXTERN void remove_from_sparse_match_c(Mat *input_mat, Mat *output_mat,
+   int lump_int, int alpha_int, PetscReal alpha);
 // Defined in PCAIR_C_Fortran_Bindings.F90 
 // External users should use the get/set routines without _c which have 
 // PetscErrorCode defined as return type, those routines are defined below this
@@ -208,6 +210,15 @@ PETSC_EXTERN void compute_diag_dom_submatrix(Mat input_mat, double max_dd_ratio,
    compute_diag_dom_submatrix_c(&input_mat, max_dd_ratio, output_mat);
 }
 
+// Restrict input_mat onto output_mat's existing sparsity pattern.
+// Auto-dispatches between the CPU and Kokkos implementations based on the
+// matrix type, so callers can mix kokkos and non-kokkos matrices freely.
+PETSC_EXTERN void remove_from_sparse_match(Mat input_mat, Mat output_mat,
+   int lump_int, int alpha_int, PetscReal alpha)
+{
+   remove_from_sparse_match_c(&input_mat, &output_mat, lump_int, alpha_int, alpha);
+}
+
 // Get routines
 
 PETSC_EXTERN PetscErrorCode PCAIRGetPrintStatsTimings(PC pc, PetscBool *input_bool)

tests/Makefile

@@ -83,8 +83,12 @@ run_tests_load_serial:
 	./ex6 -f data/mat_stream_2364 -pc_air_a_drop 1e-3 -pc_air_inverse_type power -ksp_max_it 5
 #
 	@echo ""
-	@echo "Test two concurrent AIRG solves don't corrupt per-level IS state"
+	@echo "Test two concurrent AIRG solves"
 	./ex6_two_airg -f data/mat_stream_2364
+#
+	@echo ""
+	@echo "Test AIRG can solve triangle factors from ILU factorisations"
+	./ilu_factors -f data/mat_stream_2364
 #
 	@echo ""
 	@echo "Test lAIR with GMRES polynomial smoothing for hyperbolic streaming problem"
@@ -186,8 +190,11 @@ run_tests_load_parallel:
 	@echo "Test AIRG with GMRES polynomials for hyperbolic streaming problem, matrix-free smoothing in C in parallel"
 	$(MPIEXEC) -n 2 ./ex6 -f data/mat_stream_2364 -pc_air_a_drop 1e-3 -pc_air_inverse_type power -pc_air_matrix_free_polys -ksp_max_it 5
 #
-	@echo "Test two concurrent AIRG solves don't corrupt per-level IS state in parallel"
+	@echo "Test two concurrent AIRG solves in parallel"
 	$(MPIEXEC) -n 2 ./ex6_two_airg -f data/mat_stream_2364
+#
+	@echo "Test AIRG can solve triangle factors from ILU factorisations in parallel"
+	$(MPIEXEC) -n 2 ./ilu_factors -f data/mat_stream_2364
 #
 	@echo "Test single level GMRES polynomial preconditioning for hyperbolic streaming problem in C in parallel"
 	$(MPIEXEC) -n 2 ./ex6 -f data/mat_stream_2364	-pc_type pflareinv -pc_pflareinv_type power -ksp_max_it 21