Skip to content

MatSeqAIJGetArrayRead fortran #241

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
stevendargaville merged 1 commit into from
Jun 17, 2026
Merged

MatSeqAIJGetArrayRead fortran #241

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 1 addition & 13 deletions src/C_PETSc_Interfaces.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
9 changes: 0 additions & 9 deletions src/C_PETSc_Routines.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/DDC_Module.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
30 changes: 8 additions & 22 deletions src/Gmres_Poly.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)

! ~~~~~~~~~~

Expand Down Expand Up @@ -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)

! ~~~~~~~~~~~

Expand Down
30 changes: 8 additions & 22 deletions src/Gmres_Poly_Newton.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -1231,15 +1231,14 @@ 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
MPIU_Comm :: MPI_COMM_MATRIX
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
Expand Down Expand Up @@ -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)

! ~~~~~~~~~~

Expand Down Expand Up @@ -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)

! ~~~~~~~~~~~

Expand Down
19 changes: 9 additions & 10 deletions src/MatDiagDom.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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

Expand Down
Loading