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 MatDenseGetArrayF90() is used for direct access to the 27! array that stores the dense matrix. 28! 29! Note the use of PETSC_NULL_SCALAR 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(MatDenseGetArrayF90(A,aa,ierr)) 55 56! Set matrix values directly 57 PetscCall(FillUpMatrix(m,n,aa)) 58 59 PetscCall(MatDenseRestoreArrayF90(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 return 71 end 72 73! ----------------------------------------------------------------- 74! 75! Demo2 - This subroutine demonstrates the use of user-allocated dense 76! matrix storage. 77! 78 subroutine Demo2() 79#include <petsc/finclude/petscmat.h> 80 use petscmat 81 implicit none 82 83 PetscInt n,m 84 PetscErrorCode ierr 85 parameter (m=5,n=4) 86 Mat A 87 PetscScalar aa(m,n) 88 89! Create matrix 90 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,m,n,aa,A,ierr)) 91 92! Set matrix values directly 93 PetscCall(FillUpMatrix(m,n,aa)) 94 95! Finalize matrix assembly 96 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)) 97 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)) 98 99! View matrix 100 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)) 101 102! Clean up 103 PetscCall(MatDestroy(A,ierr)) 104 return 105 end 106 107! ----------------------------------------------------------------- 108 109 subroutine FillUpMatrix(m,n,X) 110 PetscInt m,n,i,j 111 PetscScalar X(m,n) 112 113 do 10, j=1,n 114 do 20, i=1,m 115 X(i,j) = 1.0/real(i+j-1) 116 20 continue 117 10 continue 118 return 119 end 120