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