diff --git a/src/C_PETSc_Interfaces.F90 b/src/C_PETSc_Interfaces.F90 index 0b0bdb06..6997f4cc 100644 --- a/src/C_PETSc_Interfaces.F90 +++ b/src/C_PETSc_Interfaces.F90 @@ -142,20 +142,8 @@ end subroutine vecscatter_mat_reverse_end_c end interface - interface - - subroutine MatSeqAIJGetArrayF90_mine(A_array, array) & - bind(c, name="MatSeqAIJGetArrayF90_mine") - use iso_c_binding - integer(c_long_long) :: A_array - type(c_ptr) :: array + interface - end subroutine MatSeqAIJGetArrayF90_mine - - end interface - - interface - subroutine allreducesum_petscint_mine(A_array, first_int, return_int) & bind(c, name="allreducesum_petscint_mine") use iso_c_binding diff --git a/src/C_PETSc_Routines.c b/src/C_PETSc_Routines.c index e3cacdf5..96d5ece8 100644 --- a/src/C_PETSc_Routines.c +++ b/src/C_PETSc_Routines.c @@ -104,15 +104,6 @@ PETSC_INTERN void boolscatter_mat_reverse_end_c(Mat *matrix, bool *local_vals, b PetscCallVoid(PetscSFReduceEnd(a->Mvctx, MPI_C_BOOL, nonlocal_vals, local_vals, MPI_LOR)); } -// Annoying as the fortran pointer returned by MatSeqAIJGetArrayF90 seems to be -// the wrong size and that can break things sometimes -PETSC_INTERN void MatSeqAIJGetArrayF90_mine(Mat *matrix, double **array) -{ - Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*matrix)->data); - *array = a->a; - return; -} - // Annoying as MPIU_INTEGER is defined in fortran, but not if we call that fortran // routine from C - so we have to do the allreduce here PETSC_INTERN void allreducesum_petscint_mine(Mat *matrix, PetscInt first_int, PetscInt *return_int) diff --git a/src/DDC_Module.F90 b/src/DDC_Module.F90 index c3da5828..91096d82 100644 --- a/src/DDC_Module.F90 +++ b/src/DDC_Module.F90 @@ -5,7 +5,7 @@ module ddc_module use petsc_helper, only: kokkos_debug, remove_small_from_sparse, MatCreateSubMatrixWrapper use c_petsc_interfaces, only: copy_cf_markers_d2h, copy_diag_dom_ratio_d2h, ddc_kokkos, & MatDiagDomRatio_kokkos, create_cf_is_kokkos, & - vecscatter_mat_begin_c, vecscatter_mat_end_c, vecscatter_mat_restore_c, MatSeqAIJGetArrayF90_mine + vecscatter_mat_begin_c, vecscatter_mat_end_c, vecscatter_mat_restore_c use pmisr_module, only: pmisr_existing_measure_cf_markers, pmisr_existing_measure_implicit_transpose use pflare_parameters, only: C_POINT, F_POINT use matdiagdom, only: MatDiagDomRatio diff --git a/src/Gmres_Poly.F90 b/src/Gmres_Poly.F90 index 1e6e7926..85c8a9a5 100644 --- a/src/Gmres_Poly.F90 +++ b/src/Gmres_Poly.F90 @@ -2,7 +2,7 @@ module gmres_poly use petscmat use sorting, only: intersect_pre_sorted_indices_only, merge_pre_sorted - use c_petsc_interfaces, only: MatSeqAIJGetArrayF90_mine, mat_mult_powers_share_sparsity_kokkos + use c_petsc_interfaces, only: mat_mult_powers_share_sparsity_kokkos use pflare_parameters, only: PFLAREINV_POWER, PFLAREINV_ARNOLDI, PFLAREINV_NEWTON, & PFLAREINV_NEWTON_NO_EXTRA, MF_VEC_DIAG, MF_VEC_RHS, MF_VEC_TEMP, & MF_VEC_TEMP_TWO, MF_VEC_TEMP_THREE @@ -895,7 +895,6 @@ subroutine mat_mult_powers_share_sparsity_cpu(matrix, poly_order, poly_sparsity_ type(tMat) :: Ad, Ao PetscInt, dimension(:), pointer :: colmap logical :: deallocate_submatrices = .FALSE. - type(c_ptr) :: vals_c_ptr type(tMat), dimension(size(coefficients)-1), target :: matrix_powers type(tMat), pointer :: mat_sparsity_match type(int_vec), dimension(:), allocatable :: symbolic_ones @@ -905,7 +904,7 @@ subroutine mat_mult_powers_share_sparsity_cpu(matrix, poly_order, poly_sparsity_ PetscReal, dimension(:), allocatable :: vals_power_temp, vals_previous_power_temp PetscInt, dimension(:), pointer :: submatrices_ia, submatrices_ja, cols_two_ptr, cols_ptr PetscReal, dimension(:), pointer :: vals_two_ptr, vals_ptr - real(c_double), pointer :: submatrices_vals(:) + PetscScalar, pointer :: submatrices_vals(:) logical :: reuse_triggered PetscBool :: symmetric = PETSC_FALSE, inodecompressed = PETSC_FALSE, done PetscInt, parameter :: one = 1, zero = 0 @@ -1078,14 +1077,10 @@ subroutine mat_mult_powers_share_sparsity_cpu(matrix, poly_order, poly_sparsity_ print *, "Pointers not set in call to MatGetRowIJF" call MPI_Abort(MPI_COMM_WORLD, MPI_ERR_OTHER, errorcode) end if - ! Returns the wrong size pointer and can break if that size goes negative?? - !call MatSeqAIJGetArrayF90(reuse_submatrices(1),submatrices_vals,ierr); - A_array = reuse_submatrices(1)%v - ! Now we must never overwrite the values in this pointer, and we must - ! never call restore on it, see comment on top of the commented out - ! MatSeqAIJRestoreArray below - call MatSeqAIJGetArrayF90_mine(A_array, vals_c_ptr) - call c_f_pointer(vals_c_ptr, submatrices_vals, shape=[size(submatrices_ja)]) + ! Read-only access to the value array - the read variant doesn't bump + ! the matrix object state, so passing in e.g. a pc->pmat won't trigger + ! a spurious pc re-setup on every iteration + call MatSeqAIJGetArrayRead(reuse_submatrices(1), submatrices_vals, ierr) ! ~~~~~~~~~~ @@ -1260,17 +1255,8 @@ subroutine mat_mult_powers_share_sparsity_cpu(matrix, poly_order, poly_sparsity_ deallocate(symbolic_vals, symbolic_ones) end do - call MatRestoreRowIJ(reuse_submatrices(1),shift,symmetric,inodecompressed,n,submatrices_ia,submatrices_ja,done,ierr) - ! We very deliberately don't call restorearray here! - ! There is no matseqaijgetarrayread or matseqaijrestorearrayread in Fortran - ! Those routines don't increment the PetscObjectStateGet which tells petsc - ! the mat has changed. Hence above we directly access the data pointer with - ! a call to MatSeqAIJGetArrayF90_mine and then never write into it - ! If we call the restorearrayf90, that does increment the object state - ! even though we only read from the array - ! That would mean if we pass in a pc->pmat for example, just setting up a pc - ! would trigger petsc setting up the pc on every iteration of the pc - ! call MatSeqAIJRestoreArray(reuse_submatrices(1),submatrices_vals,ierr); + call MatRestoreRowIJ(reuse_submatrices(1),shift,symmetric,inodecompressed,n,submatrices_ia,submatrices_ja,done,ierr) + call MatSeqAIJRestoreArrayRead(reuse_submatrices(1), submatrices_vals, ierr) ! ~~~~~~~~~~~ diff --git a/src/Gmres_Poly_Newton.F90 b/src/Gmres_Poly_Newton.F90 index 90e81562..29ebe50d 100644 --- a/src/Gmres_Poly_Newton.F90 +++ b/src/Gmres_Poly_Newton.F90 @@ -2,7 +2,7 @@ module gmres_poly_newton use petscmat use gmres_poly - use c_petsc_interfaces, only: MatSeqAIJGetArrayF90_mine, mat_mult_powers_share_sparsity_kokkos + use c_petsc_interfaces, only: mat_mult_powers_share_sparsity_kokkos #include "petsc/finclude/petscmat.h" @@ -1231,7 +1231,6 @@ subroutine mat_mult_powers_share_sparsity_newton_cpu(matrix, poly_order, poly_sp type(tMat) :: Ad, Ao, mat_sparsity_match, mat_product_save PetscInt, dimension(:), pointer :: colmap logical :: deallocate_submatrices = .FALSE. - type(c_ptr) :: vals_c_ptr type(int_vec), dimension(:), allocatable :: symbolic_ones type(real_vec), dimension(:), allocatable :: symbolic_vals integer(c_long_long) A_array @@ -1239,7 +1238,7 @@ subroutine mat_mult_powers_share_sparsity_newton_cpu(matrix, poly_order, poly_sp PetscReal, dimension(:), allocatable :: vals_power_temp, vals_previous_power_temp, temp PetscInt, dimension(:), pointer :: submatrices_ia, submatrices_ja, cols_two_ptr, cols_ptr PetscReal, dimension(:), pointer :: vals_two_ptr, vals_ptr - real(c_double), pointer :: submatrices_vals(:) + PetscScalar, pointer :: submatrices_vals(:) logical :: reuse_triggered PetscBool :: symmetric = PETSC_FALSE, inodecompressed = PETSC_FALSE, done PetscInt, parameter :: one = 1, zero = 0 @@ -1454,14 +1453,10 @@ subroutine mat_mult_powers_share_sparsity_newton_cpu(matrix, poly_order, poly_sp print *, "Pointers not set in call to MatGetRowIJF" call MPI_Abort(MPI_COMM_WORLD, MPI_ERR_OTHER, errorcode) end if - ! Returns the wrong size pointer and can break if that size goes negative?? - !call MatSeqAIJGetArrayF90(reuse_submatrices(1),submatrices_vals,ierr); - A_array = reuse_submatrices(1)%v - ! Now we must never overwrite the values in this pointer, and we must - ! never call restore on it, see comment on top of the commented out - ! MatSeqAIJRestoreArray below - call MatSeqAIJGetArrayF90_mine(A_array, vals_c_ptr) - call c_f_pointer(vals_c_ptr, submatrices_vals, shape=[size(submatrices_ja)]) + ! Read-only access to the value array - the read variant doesn't bump + ! the matrix object state, so passing in e.g. a pc->pmat won't trigger + ! a spurious pc re-setup on every iteration + call MatSeqAIJGetArrayRead(reuse_submatrices(1), submatrices_vals, ierr) ! ~~~~~~~~~~ @@ -1763,17 +1758,8 @@ subroutine mat_mult_powers_share_sparsity_newton_cpu(matrix, poly_order, poly_sp deallocate(symbolic_vals, symbolic_ones) end do - call MatRestoreRowIJ(reuse_submatrices(1),shift,symmetric,inodecompressed,n,submatrices_ia,submatrices_ja,done,ierr) - ! We very deliberately don't call restorearray here! - ! There is no matseqaijgetarrayread or matseqaijrestorearrayread in Fortran - ! Those routines don't increment the PetscObjectStateGet which tells petsc - ! the mat has changed. Hence above we directly access the data pointer with - ! a call to MatSeqAIJGetArrayF90_mine and then never write into it - ! If we call the restorearrayf90, that does increment the object state - ! even though we only read from the array - ! That would mean if we pass in a pc->pmat for example, just setting up a pc - ! would trigger petsc setting up the pc on every iteration of the pc - ! call MatSeqAIJRestoreArray(reuse_submatrices(1),submatrices_vals,ierr); + call MatRestoreRowIJ(reuse_submatrices(1),shift,symmetric,inodecompressed,n,submatrices_ia,submatrices_ja,done,ierr) + call MatSeqAIJRestoreArrayRead(reuse_submatrices(1), submatrices_vals, ierr) ! ~~~~~~~~~~~ diff --git a/src/MatDiagDom.F90 b/src/MatDiagDom.F90 index 99cacf8f..f76b2bb9 100644 --- a/src/MatDiagDom.F90 +++ b/src/MatDiagDom.F90 @@ -4,7 +4,7 @@ module matdiagdom use petscmat use petsc_helper, only: kokkos_debug use c_petsc_interfaces, only: copy_diag_dom_ratio_d2h, MatDiagDomRatio_kokkos, & - vecscatter_mat_begin_c, vecscatter_mat_end_c, vecscatter_mat_restore_c, MatSeqAIJGetArrayF90_mine + vecscatter_mat_begin_c, vecscatter_mat_end_c, vecscatter_mat_restore_c use pflare_parameters, only: C_POINT, F_POINT #include "petsc/finclude/petscmat.h" @@ -119,14 +119,15 @@ subroutine MatDiagDomRatio_cpu(input_mat, is_fine, cf_markers_local, diag_dom_ra PetscErrorCode :: ierr integer :: errorcode, comm_size MPIU_Comm :: MPI_COMM_MATRIX - integer(c_long_long) :: A_array, vec_long, Ad_array, Ao_array + integer(c_long_long) :: A_array, vec_long type(tMat) :: Ad, Ao type(tVec) :: cf_markers_vec PetscInt, dimension(:), pointer :: is_pointer => null(), colmap => null() PetscInt, dimension(:), pointer :: ad_ia => null(), ad_ja => null(), ao_ia => null(), ao_ja => null() - PetscReal, dimension(:), pointer :: ad_vals => null(), ao_vals => null(), cf_markers_nonlocal => null() + PetscScalar, dimension(:), pointer :: ad_vals => null(), ao_vals => null() + PetscReal, dimension(:), pointer :: cf_markers_nonlocal => null() PetscReal, dimension(:), allocatable, target :: cf_markers_local_real - type(c_ptr) :: ad_vals_c_ptr, ao_vals_c_ptr, cf_markers_nonlocal_ptr + type(c_ptr) :: cf_markers_nonlocal_ptr PetscInt :: shift = 0 PetscBool :: symmetric = PETSC_FALSE, inodecompressed = PETSC_FALSE, done PetscReal :: diag_val, off_diag_sum @@ -171,17 +172,13 @@ subroutine MatDiagDomRatio_cpu(input_mat, is_fine, cf_markers_local, diag_dom_ra end if end if - Ad_array = Ad%v - call MatSeqAIJGetArrayF90_mine(Ad_array, ad_vals_c_ptr) - call c_f_pointer(ad_vals_c_ptr, ad_vals, shape=[size(ad_ja)]) + call MatSeqAIJGetArrayRead(Ad, ad_vals, ierr) ! Off-diagonal rows require a halo exchange of cf markers. ! Start and finish the scatter once, then reuse the received nonlocal markers ! while looping over all local fine rows. if (mpi) then - Ao_array = Ao%v - call MatSeqAIJGetArrayF90_mine(Ao_array, ao_vals_c_ptr) - call c_f_pointer(ao_vals_c_ptr, ao_vals, shape=[size(ao_ja)]) + call MatSeqAIJGetArrayRead(Ao, ao_vals, ierr) allocate(cf_markers_local_real(local_rows)) if (local_rows > 0) cf_markers_local_real = dble(cf_markers_local(1:local_rows)) @@ -240,8 +237,10 @@ subroutine MatDiagDomRatio_cpu(input_mat, is_fine, cf_markers_local, diag_dom_ra end if ! Restore CSR pointers before returning. + call MatSeqAIJRestoreArrayRead(Ad, ad_vals, ierr) call MatRestoreRowIJ(Ad, shift, symmetric, inodecompressed, n_ad, ad_ia, ad_ja, done, ierr) if (mpi) then + call MatSeqAIJRestoreArrayRead(Ao, ao_vals, ierr) call MatRestoreRowIJ(Ao, shift, symmetric, inodecompressed, n_ao, ao_ia, ao_ja, done, ierr) end if