[petsc-users] MatCopy question

[Frank Bramkamp]

Dear PETSc team,


I would like to create two matrices. The first matrix is based on a compact stencil (COMPACT), and the second one has an extended stencil (WIDE)
So the non-zero elements of the  COMPACT matrix is contained in the WIDE matrix as well. But the WIDE matrix shall contain some additional elements.
I try to create the COMPACT matrix first, then copy it to the WIDE matrix and then add some more terms to the WIDE matrix.
We want to create a matrix that is based on a lower order approximation of our scheme which serves as basis matrix for preconditioning,
where we want to have a more accurate jacobian which has a larger stencil, that is based on the lower approximation plus some additional non-zero elements.
We want to avoid that we have to call the same routines twice to setup the more accurate matrix, but copy the elements from the lower order (compact) matrix.


I try to use MatCopy to copy the COMPACT matrix into the WIDE matrix. Somehow I cannot figure out to do this,
since there is always a problem wirh non conforming sizes

 [0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Mat A,Mat B: global dim (4,8) (4,8)


Here is a basic example that I try.
Adding additional elements to the WIDE matrix is disabled right now, since MatCopy already gives a problem.
I thought that in MatCopy(COMPACT, WIDE, DIFFERENT_NONZERO_PATTERN, ierr)
the option DIFFERENT_NONZERO_PATTERN should enable to copy different matrix sizes into each other.
Or is that not possible to do ?!

Here the COMPACT matrix is a 4x4 matrix and the WIDE matrix 8x8, for a basic test.

Thanks, Frank Bramkamp




program petsc_matrix_example
#include <petsc/finclude/petsc.h>
  use petsc
  implicit none

  PetscErrorCode :: ierr
  Mat :: COMPACT, WIDE
  PetscInt :: m, n, MM, NN, d_nz
  PetscInt :: i, j
  PetscScalar :: additional_value

! Initialize PETSc
  call PetscInitialize(PETSC_NULL_CHARACTER, ierr)

! Set dimensions for COMPACT matrix
  m = 4  ! Local rows
  n = 4  ! Local columns
  MM = 4  ! Global rows
  NN = 4  ! Global columns
  d_nz = 4 ! Estimated non-zeros per row (diagonal) o_nz = 0 ! Estimated non-zeros per row (off-diagonal)

! Create COMPACT matrix
  call MatCreateSeqAIJ(PETSC_COMM_SELF, MM, NN, PETSC_DECIDE, PETSC_NULL_INTEGER, COMPACT, ierr)

! Set some values in COMPACT matrix (example) do i = 0, M-1
  do i = 0, NN-1
     !do j = max(0, i-1), min(NN-1, i+1)
     do j = 0,NN-1
      call MatSetValue(COMPACT, i, j, 1.0_PETSC_REAL_KIND, INSERT_VALUES, ierr)
    end do
  end do

! Assemble COMPACT matrix
  call MatAssemblyBegin(COMPACT, MAT_FINAL_ASSEMBLY, ierr)
  call MatAssemblyEnd(COMPACT, MAT_FINAL_ASSEMBLY, ierr)

! Set dimensions for WIDE matrix
  n = 8  ! Increase local columns
  NN = 8  ! Increase global columns
  m  = 8
  MM = 8
  d_nz = 8
  
! Create WIDE matrix
  call MatCreateSeqAIJ(PETSC_COMM_SELF, MM, NN, PETSC_DECIDE, PETSC_NULL_INTEGER, WIDE, ierr)

  ! Assemble WIDE matrix: MUST WE ASSMBLE THE MATRIX BEFORE MatCOPY ?!
  call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr) 
  call MatAssemblyEnd(WIDE, MAT_FINAL_ASSEMBLY, ierr)
  
! Copy elements from COMPACT to WIDE
  call MatCopy(COMPACT, WIDE, DIFFERENT_NONZERO_PATTERN, ierr)



  !call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr) 
  
  ! Add additional elements to WIDE matrix
!!$  additional_value = 10.0
!!$
!!$  ! original
!!$  !do i = 0, MM-1
!!$  !  do j = NN/2, NN-1
!!$  !    call MatSetValue(WIDE, i, j, additional_value, INSERT_VALUES, ierr)
!!$  !  end do
!!$ !end do
!!$
!!$  do i = 0, NN-1
!!$    do j = max(0, i-1), min(NN-1, i+1)
!!$      call MatSetValue(WIDE, i, j, 10.0_PETSC_REAL_KIND, ADD_VALUES, ierr)
!!$    end do
!!$  end do
 

! Assemble WIDE matrix again
  call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr)
  call MatAssemblyEnd(WIDE, MAT_FINAL_ASSEMBLY, ierr)

! View matrices (optional)
  call PetscPrintf(PETSC_COMM_WORLD, "COMPACT matrix:\n", ierr)
  call MatView(COMPACT, PETSC_VIEWER_STDOUT_WORLD, ierr)
  call PetscPrintf(PETSC_COMM_WORLD, "WIDE matrix:\n", ierr)
  call MatView(WIDE, PETSC_VIEWER_STDOUT_WORLD, ierr)

! Clean up
  call MatDestroy(COMPACT, ierr)
  call MatDestroy(WIDE, ierr)

! Finalize PETSc
  call PetscFinalize(ierr)

end program petsc_matrix_example