Newton assemble

Makefile

@@ -138,7 +138,7 @@ OBJS := $(OBJS) $(SRCDIR)/PETSc_Helper.o \
 		  $(SRCDIR)/PCPFLAREINV.o	
 
 # Define a variable containing all the tests
-export TEST_TARGETS = ex12f ex6f ex6f_getcoeffs ex6 adv_1d adv_diff_2d ex6_cf_splitting adv_diff_cg_supg matrandom matrandom_check_reset
+export TEST_TARGETS = ex12f ex6f ex6f_getcoeffs ex6 adv_1d adv_diff_2d ex6_cf_splitting adv_diff_cg_supg matrandom matrandom_check_reset ex12f_gmres_poly mat_diag
 # Include kokkos examples
 ifeq ($(PETSC_HAVE_KOKKOS),1)
 export TEST_TARGETS := $(TEST_TARGETS) adv_1dk

README.md

@@ -37,8 +37,8 @@ PCPFLAREINV contains methods for computing approximate inverses, most of which c
    | ------------- | -- | ------------- | -- |
    | power  |  PFLAREINV_POWER  | GMRES polynomial, applied as a mononomial, with coefficients computed with a power basis  | Yes |
    | arnoldi  |  PFLAREINV_ARNOLDI  | GMRES polynomial, applied as a mononomial, with coefficients computed with an Arnoldi method  | Yes |
-   | newton  |  PFLAREINV_NEWTON  | GMRES polynomial, applied as a Newton polynomial, with roots computed with an Arnoldi method and with extra roots added for stability  | Yes |
-   | newton_no_extra  |  PFLAREINV_NEWTON_NO_EXTRA  | GMRES polynomial, applied as a Newton polynomial, with roots computed with an Arnoldi method and with no extra roots added   | Yes |
+   | newton  |  PFLAREINV_NEWTON  | GMRES polynomial, applied as a Newton polynomial, with roots computed with an Arnoldi method and with extra roots added for stability  | Matrix-free: Yes Assembled: No |
+   | newton_no_extra  |  PFLAREINV_NEWTON_NO_EXTRA  | GMRES polynomial, applied as a Newton polynomial, with roots computed with an Arnoldi method and with no extra roots added   | Matrix-free: Yes Assembled: No |
    | neumann  |  PFLAREINV_NEUMANN  | Neumann polynomial  | Yes |
    | sai  |  PFLAREINV_SAI  | Sparse approximate inverse  | No |
    | isai  |  PFLAREINV_ISAI  | Incomplete sparse approximate inverse (equivalent to a one-level RAS)  | No |

dockerfiles/Dockerfile_intel

@@ -28,6 +28,9 @@ RUN source /opt/intel/oneapi/setvars.sh && \
     sed -i '/^run_tests_load_parallel:/,/^run_tests_no_load_serial:/s/-pc_pflareinv_type power/-pc_pflareinv_type arnoldi/g' tests/Makefile && \
     sed -i '/^run_tests_no_load_parallel:/,/^run_tests_medium_serial:/s/-pc_air_inverse_type power/-pc_air_inverse_type arnoldi/g' tests/Makefile && \
     sed -i '/^run_tests_no_load_parallel:/,/^run_tests_medium_serial:/s/-pc_pflareinv_type power/-pc_pflareinv_type arnoldi/g' tests/Makefile && \
+    echo "Disabling power basis in parallel test ex12f_gmres_poly in tests/Makefile" && \
+    sed -i '/^run_tests_load_parallel:/,/^run_tests_no_load_serial:/s/ex12f_gmres_poly/ex12f_gmres_poly -no_power/g' tests/Makefile && \
+    sed -i '/^run_tests_no_load_parallel:/,/^run_tests_medium_serial:/s/ex12f_gmres_poly/ex12f_gmres_poly -no_power/g' tests/Makefile && \    
     make -j2 && make -j2 check && \
     make -j2 tests
 

src/AIR_MG_Setup.F90

@@ -424,7 +424,7 @@ subroutine setup_air_pcmg(amat, pmat, air_data, pcmg_input)
          end if
 
          ! Convert Aff to a matdiagonal type
-         ! Haven't rewritten sai to take advantage of matdiagonal
+         ! Haven't rewritten some inverse types to take advantage of matdiagonal
          if (aff_diag .AND. &
                   inverse_type_aff /= PFLAREINV_SAI .AND. &
                   inverse_type_aff /= PFLAREINV_ISAI) then

src/Approx_Inverse_Setup.F90

@@ -70,7 +70,12 @@ subroutine calculate_and_build_approximate_inverse(matrix, inverse_type, &
             allocate(coefficients(poly_order + 1, 1))
          end if
       else
-         coefficients => coefficients_stack
+         if (inverse_type == PFLAREINV_NEWTON .OR. inverse_type == PFLAREINV_NEWTON_NO_EXTRA) then
+            ! Newton basis needs storage for real and imaginary roots
+            allocate(coefficients(poly_order + 1, 2))
+         else         
+            coefficients => coefficients_stack
+         end if
       end if
 
       ! This is diabolical - In petsc 3.22, they changed the way to test for 
@@ -303,13 +308,9 @@ subroutine finish_approximate_inverse(matrix, inverse_type, &
       ! Gmres polynomial with newton basis
       else if (inverse_type == PFLAREINV_NEWTON .OR. inverse_type == PFLAREINV_NEWTON_NO_EXTRA) then
 
-         if (.NOT. matrix_free) then
-            print *, "GMRES polynomial with Newton basis must be applied matrix-free"
-            call MPI_Abort(MPI_COMM_WORLD, MPI_ERR_OTHER, errorcode)
-         end if
-
          call build_gmres_polynomial_newton_inverse(matrix, poly_order, &
                            coefficients, &
+                           inverse_sparsity_order, matrix_free, reuse_mat, reuse_submatrices, &
                            inv_matrix)         
 
       ! Neumann polynomial

src/Gmres_Poly.F90

@@ -851,7 +851,6 @@ subroutine mat_mult_powers_share_sparsity_cpu(matrix, poly_order, poly_sparsity_
 
       ! Compute matrix powers c = coeff(1) * I + coeff(2) * A + coeff(3) * A^2 + coeff(4) * A^3 + ... 
       ! where a c and the powers all share the same sparsity as the power input in poly_sparsity_order
-      ! Assuming cmat has not been built/allocated
       ! This also finishes the async comms required to compute the gmres poly coefficients if buffers%request is allocated
 
       ! ~~~~~~~~~~
@@ -1655,7 +1654,7 @@ subroutine build_gmres_polynomial_inverse_0th_order_sparsity(matrix, poly_order,
       integer :: order
       PetscErrorCode :: ierr      
       logical :: reuse_triggered
-      type(tVec) :: rhs_copy, diag_vec, power_vec
+      type(tVec) :: inv_vec, diag_vec, power_vec
       ! ~~~~~~
 
       reuse_triggered = .NOT. PetscObjectIsNull(inv_matrix)
@@ -1666,9 +1665,9 @@ subroutine build_gmres_polynomial_inverse_0th_order_sparsity(matrix, poly_order,
 
       ! This stores D^order
       if (.NOT. reuse_triggered) then
-         call VecDuplicate(diag_vec, rhs_copy, ierr)
+         call VecDuplicate(diag_vec, inv_vec, ierr)
       else
-         call MatDiagonalGetDiagonal(inv_matrix, rhs_copy, ierr)
+         call MatDiagonalGetDiagonal(inv_matrix, inv_vec, ierr)
       end if
       call VecCopy(diag_vec, power_vec, ierr)
 
@@ -1678,22 +1677,22 @@ subroutine build_gmres_polynomial_inverse_0th_order_sparsity(matrix, poly_order,
       call finish_gmres_polynomial_coefficients_power(poly_order, buffers, coefficients)         
 
       ! Set: alpha_0 * I
-      call VecSet(rhs_copy, coefficients(1), ierr)
+      call VecSet(inv_vec, coefficients(1), ierr)
 
       ! Calculate: alpha_0 * I + alpha_1 * D + alpha_2 * D^2
       do order = 1, poly_order
-         call VecAXPY(rhs_copy, coefficients(order+1), power_vec, ierr)
+         call VecAXPY(inv_vec, coefficients(order+1), power_vec, ierr)
          ! Compute power_vec = power_vec * D
          if (order /= poly_order) call VecPointwiseMult(power_vec, power_vec, diag_vec, ierr)
       end do
 
       ! We may be reusing with the same sparsity
       if (.NOT. reuse_triggered) then
-         ! The matrix takes ownership of rhs_copy and increases ref counter
-         call MatCreateDiagonal(rhs_copy, inv_matrix, ierr)
-         call VecDestroy(rhs_copy, ierr)
+         ! The matrix takes ownership of inv_vec and increases ref counter
+         call MatCreateDiagonal(inv_vec, inv_matrix, ierr)
+         call VecDestroy(inv_vec, ierr)
       else
-         call MatDiagonalRestoreDiagonal(inv_matrix, rhs_copy, ierr)
+         call MatDiagonalRestoreDiagonal(inv_matrix, inv_vec, ierr)
       end if  
 
       call VecDestroy(diag_vec, ierr)