1! 2! 3! This program demonstrates use of PETSc dense matrices. 4! 5#include <petsc/finclude/petscsys.h> 6#include <petsc/finclude/petscmat.h> 7program main 8 use petscsys 9 implicit none 10 11 PetscErrorCode ierr 12 13 PetscCallA(PetscInitialize(ierr)) 14 15! Demo of PETSc-allocated dense matrix storage 16 call Demo1() 17 18! Demo of user-allocated dense matrix storage 19 call Demo2() 20 21 PetscCallA(PetscFinalize(ierr)) 22end 23 24! ----------------------------------------------------------------- 25! 26! Demo1 - This subroutine demonstrates the use of PETSc-allocated dense 27! matrix storage. Here MatDenseGetArray() is used for direct access to the 28! array that stores the dense matrix. 29! 30! Note the use of PETSC_NULL_SCALAR_ARRAY in MatCreateSeqDense() to indicate that no 31! storage is being provided by the user. (Do NOT pass a zero in that 32! location.) 33! 34subroutine Demo1() 35 use petscmat 36 implicit none 37 38 Mat A 39 PetscInt n, m 40 PetscErrorCode ierr 41 PetscScalar, pointer :: aa(:, :) 42 43 n = 4 44 m = 5 45 46! Create matrix 47 48 PetscCall(MatCreate(PETSC_COMM_SELF, A, ierr)) 49 PetscCall(MatSetSizes(A, m, n, m, n, ierr)) 50 PetscCall(MatSetType(A, MATSEQDENSE, ierr)) 51 PetscCall(MatSetUp(A, ierr)) 52 53! Access array storage 54 PetscCall(MatDenseGetArray(A, aa, ierr)) 55 56! Set matrix values directly 57 PetscCall(FillUpMatrix(m, n, aa)) 58 59 PetscCall(MatDenseRestoreArray(A, aa, ierr)) 60 61! Finalize matrix assembly 62 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) 63 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) 64 65! View matrix 66 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr)) 67 68! Clean up 69 PetscCall(MatDestroy(A, ierr)) 70end 71 72! ----------------------------------------------------------------- 73! 74! Demo2 - This subroutine demonstrates the use of user-allocated dense 75! matrix storage. 76! 77subroutine Demo2() 78 use petscmat 79 implicit none 80 81 PetscInt n, m 82 PetscErrorCode ierr 83 parameter(m=5, n=4) 84 Mat A 85 PetscScalar aa(m, n) 86 87! Create matrix 88 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, n, aa, A, ierr)) 89 90! Set matrix values directly 91 PetscCall(FillUpMatrix(m, n, aa)) 92 93! Finalize matrix assembly 94 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) 95 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) 96 97! View matrix 98 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr)) 99 100! Clean up 101 PetscCall(MatDestroy(A, ierr)) 102end 103 104! ----------------------------------------------------------------- 105 106subroutine FillUpMatrix(m, n, X) 107 PetscInt m, n, i, j 108 PetscScalar X(m, n) 109 110 do j = 1, n 111 do i = 1, m 112 X(i, j) = 1.0/real(i + j - 1) 113 end do 114 end do 115end 116