1! Test code contributed by Thibaut Appel <t.appel17@imperial.ac.uk> 2#include <petsc/finclude/petscmat.h> 3program test_assembly 4 5 use petscmat 6 use ISO_Fortran_Env, only: real64 7 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 14 15 PetscInt cols(1), rows(1) 16 PetscScalar vals(1) 17 18 PetscCallA(PetscInitialize(ierr)) 19 20 i0 = 0 21 i1 = 1 22 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)) 26 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)) 32 33 PetscCallA(MatGetOwnershipRange(L, istart, iend, ierr)) 34 35 ! assembling a diagonal matrix 36 do row = istart, iend - 1 37 38 cols = [row]; vals = [one]; rows = [row] 39 PetscCallA(MatSetValues(L, i1, rows, i1, cols, vals, ADD_VALUES, ierr)) 40 41 end do 42 43 PetscCallA(MatAssemblyBegin(L, MAT_FINAL_ASSEMBLY, ierr)) 44 PetscCallA(MatAssemblyEnd(L, MAT_FINAL_ASSEMBLY, ierr)) 45 46 PetscCallA(MatSetOption(L, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr)) 47 48 !PetscCallA(MatZeroEntries(L,ierr)) 49 50 ! assembling a diagonal matrix, adding a zero value to non-diagonal part 51 do row = istart, iend - 1 52 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)) 61 62 end do 63 64 PetscCallA(MatAssemblyBegin(L, MAT_FINAL_ASSEMBLY, ierr)) 65 PetscCallA(MatAssemblyEnd(L, MAT_FINAL_ASSEMBLY, ierr)) 66 PetscCallA(MatDestroy(L, ierr)) 67 68 PetscCallA(PetscFinalize(ierr)) 69 70end program test_assembly 71 72!/*TEST 73! 74! build: 75! requires: complex 76! 77! test: 78! nsize: 2 79! output_file: output/empty.out 80! 81!TEST*/ 82