1! 2! 3! This program demonstrates use of PETSc dense matrices. 4! 5 program main 6#include <petsc/finclude/petscsys.h> 7 use petscsys 8 implicit none 9 10 PetscErrorCode ierr 11 12 PetscCallA(PetscInitialize(ierr)) 13 14! Demo of PETSc-allocated dense matrix storage 15 call Demo1() 16 17! Demo of user-allocated dense matrix storage 18 call Demo2() 19 20 PetscCallA(PetscFinalize(ierr)) 21 end 22 23! ----------------------------------------------------------------- 24! 25! Demo1 - This subroutine demonstrates the use of PETSc-allocated dense 26! matrix storage. Here MatDenseGetArray() is used for direct access to the 27! array that stores the dense matrix. 28! 29! Note the use of PETSC_NULL_SCALAR_ARRAY in MatCreateSeqDense() to indicate that no 30! storage is being provided by the user. (Do NOT pass a zero in that 31! location.) 32! 33 subroutine Demo1() 34#include <petsc/finclude/petscmat.h> 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)) 70 end 71 72! ----------------------------------------------------------------- 73! 74! Demo2 - This subroutine demonstrates the use of user-allocated dense 75! matrix storage. 76! 77 subroutine Demo2() 78#include <petsc/finclude/petscmat.h> 79 use petscmat 80 implicit none 81 82 PetscInt n,m 83 PetscErrorCode ierr 84 parameter (m=5,n=4) 85 Mat A 86 PetscScalar aa(m,n) 87 88! Create matrix 89 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,m,n,aa,A,ierr)) 90 91! Set matrix values directly 92 PetscCall(FillUpMatrix(m,n,aa)) 93 94! Finalize matrix assembly 95 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)) 96 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)) 97 98! View matrix 99 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)) 100 101! Clean up 102 PetscCall(MatDestroy(A,ierr)) 103 end 104 105! ----------------------------------------------------------------- 106 107 subroutine FillUpMatrix(m,n,X) 108 PetscInt m,n,i,j 109 PetscScalar X(m,n) 110 111 do 10, j=1,n 112 do 20, i=1,m 113 X(i,j) = 1.0/real(i+j-1) 114 20 continue 115 10 continue 116 end 117