1: !     Test code contributed by Thibaut Appel <t.appel17@imperial.ac.uk>

  3: program test_assembly

  5: #include <petsc/finclude/petscmat.h>

  7:   use PetscMat
  8:   use ISO_Fortran_Env, only: real64

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

 17:   PetscInt cols(1), rows(1)
 18:   PetscScalar vals(1)

 20:   PetscCallA(PetscInitialize(ierr))

 22:   i0 = 0
 23:   i1 = 1

 25:   PetscCallA(MatCreate(PETSC_COMM_WORLD, L, ierr))
 26:   PetscCallA(MatSetType(L, MATAIJ, ierr))
 27:   PetscCallA(MatSetSizes(L, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))

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

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

 37:   ! assembling a diagonal matrix
 38:   do row = istart, iend - 1

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

 43:   end do

 45:   PetscCallA(MatAssemblyBegin(L, MAT_FINAL_ASSEMBLY, ierr))
 46:   PetscCallA(MatAssemblyEnd(L, MAT_FINAL_ASSEMBLY, ierr))

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

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

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

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

 64:   end do

 66:   PetscCallA(MatAssemblyBegin(L, MAT_FINAL_ASSEMBLY, ierr))
 67:   PetscCallA(MatAssemblyEnd(L, MAT_FINAL_ASSEMBLY, ierr))
 68:   PetscCallA(MatDestroy(L, ierr))

 70:   PetscCallA(PetscFinalize(ierr))

 72: end program test_assembly

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