1: !     Test code contributed by Thibaut Appel <t.appel17@imperial.ac.uk>
  2: #include <petsc/finclude/petscmat.h>
  3: program test_assembly

  5:   use petscmat
  6:   use ISO_Fortran_Env, only: real64

  8:   implicit none
  9:   PetscInt, parameter :: wp = real64, n = 10
 10:   PetscScalar, parameter :: zero = 0.0, one = 1.0
 11:   Mat      :: L
 12:   PetscInt :: istart, iend, row, i1, i0
 13:   PetscErrorCode :: ierr

 15:   PetscInt cols(1), rows(1)
 16:   PetscScalar vals(1)

 18:   PetscCallA(PetscInitialize(ierr))

 20:   i0 = 0
 21:   i1 = 1

 23:   PetscCallA(MatCreate(PETSC_COMM_WORLD, L, ierr))
 24:   PetscCallA(MatSetType(L, MATAIJ, ierr))
 25:   PetscCallA(MatSetSizes(L, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))

 27:   PetscCallA(MatSeqAIJSetPreallocation(L, i1, PETSC_NULL_INTEGER_ARRAY, ierr))
 28:   PetscCallA(MatMPIAIJSetPreallocation(L, i1, PETSC_NULL_INTEGER_ARRAY, i0, PETSC_NULL_INTEGER_ARRAY, ierr)) ! No allocated non-zero in off-diagonal part
 29:   PetscCallA(MatSetOption(L, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE, ierr))
 30:   PetscCallA(MatSetOption(L, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE, ierr))
 31:   PetscCallA(MatSetOption(L, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE, ierr))

 33:   PetscCallA(MatGetOwnershipRange(L, istart, iend, ierr))

 35:   ! assembling a diagonal matrix
 36:   do row = istart, iend - 1

 38:     cols = [row]; vals = [one]; rows = [row]
 39:     PetscCallA(MatSetValues(L, i1, rows, i1, cols, vals, ADD_VALUES, ierr))

 41:   end do

 43:   PetscCallA(MatAssemblyBegin(L, MAT_FINAL_ASSEMBLY, ierr))
 44:   PetscCallA(MatAssemblyEnd(L, MAT_FINAL_ASSEMBLY, ierr))

 46:   PetscCallA(MatSetOption(L, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr))

 48:   !PetscCallA(MatZeroEntries(L,ierr))

 50:   ! assembling a diagonal matrix, adding a zero value to non-diagonal part
 51:   do row = istart, iend - 1

 53:     if (row == 0) then
 54:       cols = [n - 1]
 55:       vals = [zero]
 56:       rows = [row]
 57:       PetscCallA(MatSetValues(L, i1, rows, i1, cols, vals, ADD_VALUES, ierr))
 58:     end if
 59:     cols = [row]; vals = [one]; rows = [row]
 60:     PetscCallA(MatSetValues(L, i1, rows, i1, cols, vals, ADD_VALUES, ierr))

 62:   end do

 64:   PetscCallA(MatAssemblyBegin(L, MAT_FINAL_ASSEMBLY, ierr))
 65:   PetscCallA(MatAssemblyEnd(L, MAT_FINAL_ASSEMBLY, ierr))
 66:   PetscCallA(MatDestroy(L, ierr))

 68:   PetscCallA(PetscFinalize(ierr))

 70: end program test_assembly

 72: !/*TEST
 73: !
 74: !   build:
 75: !      requires: complex
 76: !
 77: !   test:
 78: !      nsize: 2
 79: !      output_file: output/empty.out
 80: !
 81: !TEST*/