1#include <petsc/finclude/petscmat.h> 2program main 3 use petscmat 4 5 implicit none 6 7 Mat A 8 PetscInt, parameter :: n = 5, m = 5 9 PetscScalar, parameter :: two = 2.0, one = 1.0 10 PetscInt, pointer, dimension(:) :: dnnz, onnz 11 PetscInt :: i, rstart, rend, M1, N1 12 PetscErrorCode ierr 13 14 PetscCallA(PetscInitialize(ierr)) 15 16 allocate (dnnz(0:m - 1)) 17 allocate (onnz(0:m - 1)) 18 19 do i = 0, m - 1 20 dnnz(i) = 1 21 onnz(i) = 1 22 end do 23 24 PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, m, n, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DECIDE, dnnz, PETSC_DECIDE, onnz, A, ierr)) 25 PetscCallA(MatSetFromOptions(A, ierr)) 26 PetscCallA(MatSetUp(A, ierr)) 27 deallocate (dnnz) 28 deallocate (onnz) 29 30 !/* This assembly shrinks memory because we do not insert enough number of values */ 31 PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) 32 PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) 33 34 !/* MatResetPreallocation restores the memory required by users */ 35 PetscCallA(MatResetPreallocation(A, ierr)) 36 PetscCallA(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE, ierr)) 37 PetscCallA(MatGetOwnershipRange(A, rstart, rend, ierr)) 38 PetscCallA(MatGetSize(A, M1, N1, ierr)) 39 do i = rstart, rend - 1 40 PetscCallA(MatSetValue(A, i, i, two, INSERT_VALUES, ierr)) 41 if (rend < N1) PetscCallA(MatSetValue(A, i, rend, one, INSERT_VALUES, ierr)) 42 end do 43 PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) 44 PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) 45 PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr)) 46 PetscCallA(MatDestroy(A, ierr)) 47 PetscCallA(PetscFinalize(ierr)) 48 49end program 50 51!/*TEST 52! 53! test: 54! suffix: 1 55! output_file: output/ex4_1.out 56! 57! test: 58! suffix: 2 59! nsize: 2 60! output_file: output/ex4_2.out 61! 62!TEST*/ 63