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