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. The user declares an array (aa(1)) 28! and index variable (ia), which are then used together to manipulate 29! the array contents. 30! 31! Note the use of PETSC_NULL_SCALAR in MatCreateSeqDense() to indicate that no 32! storage is being provided by the user. (Do NOT pass a zero in that 33! location.) 34! 35 subroutine Demo1() 36#include <petsc/finclude/petscmat.h> 37 use petscmat 38 implicit none 39 40 Mat A 41 PetscInt n,m 42 PetscErrorCode ierr 43 PetscScalar aa(1) 44 PetscOffset ia 45 46 n = 4 47 m = 5 48 49! Create matrix 50 51 PetscCall(MatCreate(PETSC_COMM_SELF,A,ierr)) 52 PetscCall(MatSetSizes(A,m,n,m,n,ierr)) 53 PetscCall(MatSetType(A,MATSEQDENSE,ierr)) 54 PetscCall(MatSetUp(A,ierr)) 55 56! Access array storage 57 PetscCall(MatDenseGetArray(A,aa,ia,ierr)) 58 59! Set matrix values directly 60 PetscCall(FillUpMatrix(m,n,aa(ia+1)) 61 62 PetscCall(MatDenseRestoreArray(A,aa,ia,ierr)) 63 64! Finalize matrix assembly 65 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)) 66 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)) 67 68! View matrix 69 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)) 70 71! Clean up 72 PetscCall(MatDestroy(A,ierr)) 73 return 74 end 75 76! ----------------------------------------------------------------- 77! 78! Demo2 - This subroutine demonstrates the use of user-allocated dense 79! matrix storage. 80! 81 subroutine Demo2() 82#include <petsc/finclude/petscmat.h> 83 use petscmat 84 implicit none 85 86 PetscInt n,m 87 PetscErrorCode ierr 88 parameter (m=5,n=4) 89 Mat A 90 PetscScalar aa(m,n) 91 92! Create matrix 93 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,m,n,aa,A,ierr)) 94 95! Set matrix values directly 96 PetscCall(FillUpMatrix(m,n,aa) 97 98! Finalize matrix assembly 99 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)) 100 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)) 101 102! View matrix 103 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)) 104 105! Clean up 106 PetscCall(MatDestroy(A,ierr)) 107 return 108 end 109 110! ----------------------------------------------------------------- 111 112 subroutine FillUpMatrix(m,n,X) 113 PetscInt m,n,i,j 114 PetscScalar X(m,n) 115 116 do 10, j=1,n 117 do 20, i=1,m 118 X(i,j) = 1.0/real(i+j-1) 119 20 continue 120 10 continue 121 return 122 end 123