Actual source code: ex36f.F
petsc-3.12.2 2019-11-22
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
12: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
13: if (ierr .ne. 0) then
14: print*,'Unable to initialize PETSc'
15: stop
16: endif
18: ! Demo of PETSc-allocated dense matrix storage
19: call Demo1()
21: ! Demo of user-allocated dense matrix storage
22: call Demo2()
24: call PetscFinalize(ierr)
25: end
27: ! -----------------------------------------------------------------
28: !
29: ! Demo1 - This subroutine demonstrates the use of PETSc-allocated dense
30: ! matrix storage. Here MatDenseGetArray() is used for direct access to the
31: ! array that stores the dense matrix. The user declares an array (aa(1))
32: ! and index variable (ia), which are then used together to manipulate
33: ! the array contents.
34: !
35: ! Note the use of PETSC_NULL_SCALAR in MatCreateSeqDense() to indicate that no
36: ! storage is being provided by the user. (Do NOT pass a zero in that
37: ! location.)
38: !
39: subroutine Demo1()
40: #include <petsc/finclude/petscmat.h>
41: use petscmat
42: implicit none
44: Mat A
45: PetscInt n,m
46: PetscErrorCode ierr
47: PetscScalar aa(1)
48: PetscOffset ia
50: n = 4
51: m = 5
53: ! Create matrix
55: ! call MatCreateSeqDense(PETSC_COMM_SELF,m,n,PETSC_NULL_SCALAR, &
56: ! & A,ierr)
58: ! Using MatCreate() instead of MatCreateSeqDense() as above to avoid Nag F90 errors
59: ! However both cases are equivalent
61: call MatCreate(PETSC_COMM_SELF,A,ierr)
62: call MatSetSizes(A,m,n,m,n,ierr)
63: call MatSetType(A,MATSEQDENSE,ierr)
64: call MatSetUp(A,ierr)
66: ! Access array storage
67: call MatDenseGetArray(A,aa,ia,ierr)
69: ! Set matrix values directly
70: call FillUpMatrix(m,n,aa(ia+1))
72: call MatDenseRestoreArray(A,aa,ia,ierr)
74: ! Finalize matrix assembly
75: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
76: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
78: ! View matrix
79: call MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)
81: ! Clean up
82: call MatDestroy(A,ierr)
83: return
84: end
86: ! -----------------------------------------------------------------
87: !
88: ! Demo2 - This subroutine demonstrates the use of user-allocated dense
89: ! matrix storage.
90: !
91: subroutine Demo2()
92: #include <petsc/finclude/petscmat.h>
93: use petscmat
94: implicit none
96: PetscInt n,m
97: PetscErrorCode ierr
98: parameter (m=5,n=4)
99: Mat A
100: PetscScalar aa(m,n)
102: ! Create matrix
103: call MatCreateSeqDense(PETSC_COMM_SELF,m,n,aa,A,ierr)
105: ! Set matrix values directly
106: call FillUpMatrix(m,n,aa)
108: ! Finalize matrix assembly
109: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
110: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
112: ! View matrix
113: call MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)
115: ! Clean up
116: call MatDestroy(A,ierr)
117: return
118: end
120: ! -----------------------------------------------------------------
122: subroutine FillUpMatrix(m,n,X)
123: PetscInt m,n,i,j
124: PetscScalar X(m,n)
126: do 10, j=1,n
127: do 20, i=1,m
128: X(i,j) = 1.0/real(i+j-1)
129: 20 continue
130: 10 continue
131: return
132: end