1: !
  2: !
  3: !   This program demonstrates use of PETSc dense matrices.
  4: !
  5: #include <petsc/finclude/petscmat.h>
  6: module ex36fmodule
  7:   use petscmat
  8:   implicit none
  9: ! -----------------------------------------------------------------
 10: !
 11: !  Demo1 -  This subroutine demonstrates the use of PETSc-allocated dense
 12: !  matrix storage.  Here MatDenseGetArray() is used for direct access to the
 13: !  array that stores the dense matrix.
 14: !
 15: !  Note the use of PETSC_NULL_SCALAR_ARRAY in MatCreateSeqDense() to indicate that no
 16: !  storage is being provided by the user. (Do NOT pass a zero in that
 17: !  location.)
 18: !
 19: contains
 20:   subroutine Demo1(m, n)

 22:     PetscInt, intent(in) :: m, n
 23:     Mat A
 24:     PetscScalar, pointer :: aa(:, :)
 25:     PetscErrorCode ierr

 27:     ! Create matrix
 28:     PetscCall(MatCreate(PETSC_COMM_SELF, A, ierr))
 29:     PetscCall(MatSetSizes(A, m, n, m, n, ierr))
 30:     PetscCall(MatSetType(A, MATSEQDENSE, ierr))
 31:     PetscCall(MatSetUp(A, ierr))

 33:     ! Access array storage
 34:     PetscCall(MatDenseGetArray(A, aa, ierr))

 36:     ! Set matrix values directly
 37:     PetscCall(FillUpMatrix(m, n, aa))

 39:     PetscCall(MatDenseRestoreArray(A, aa, ierr))

 41:     ! View matrix
 42:     PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))

 44:     ! Clean up
 45:     PetscCall(MatDestroy(A, ierr))
 46:   end subroutine Demo1

 48: ! -----------------------------------------------------------------
 49: !
 50: !  Demo2 -  This subroutine demonstrates the use of user-provided dense
 51: !  matrix storage. Using allocate (typically heap memory)
 52: !
 53:   subroutine Demo2(m, n)

 55:     PetscInt, intent(in) :: m, n
 56:     Mat A
 57:     PetscScalar, pointer :: aa(:, :)
 58:     PetscErrorCode ierr

 60:     allocate (aa(m, n))

 62:     ! Create matrix
 63:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, n, aa, A, ierr))

 65:     ! Set matrix values directly
 66:     PetscCall(FillUpMatrix(m, n, aa))

 68:     ! View matrix
 69:     PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))

 71:     ! Clean up
 72:     PetscCall(MatDestroy(A, ierr))
 73:     deallocate (aa)
 74:   end subroutine Demo2

 76: ! -----------------------------------------------------------------
 77: !
 78: !  Demo3 -  This subroutine demonstrates the use of user-provided dense
 79: !  matrix storage. Using fixed dimensions (typically stack memory)
 80: !
 81:   subroutine Demo3(m, n)

 83:     PetscInt, intent(in) :: m, n
 84:     Mat A
 85:     PetscScalar :: aa(m, n)
 86:     PetscErrorCode ierr

 88:     ! Create matrix
 89:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, n, aa, A, ierr))

 91:     ! Set matrix values directly
 92:     PetscCall(FillUpMatrix(m, n, aa))

 94:     ! View matrix
 95:     PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))

 97:     ! Clean up
 98:     PetscCall(MatDestroy(A, ierr))
 99:     if (aa(1, 1) == 2.0) then
100:       print *, 'Error in a(1,1)'
101:     end if
102:   end subroutine Demo3

104: ! -----------------------------------------------------------------

106:   subroutine FillUpMatrix(m, n, X)
107:     PetscInt, intent(in) :: m, n
108:     PetscScalar, intent(out) :: X(m, n)
109:     PetscInt i, j

111:     do j = 1, n
112:       do i = 1, m
113:         X(i, j) = 1.0/real(i + j - 1)
114:       end do
115:     end do
116:   end subroutine FillUpMatrix
117: end module ex36fmodule

119: program main
120:   use ex36fmodule
121:   implicit none

123:   PetscErrorCode ierr
124:   PetscInt, parameter :: m = 5, n = 4

126:   PetscCallA(PetscInitialize(ierr))

128:   ! Demo of PETSc-allocated dense matrix storage
129:   call Demo1(m, n)

131:   ! Demo of user-provided dense matrix storage (heap)
132:   call Demo2(m, n)

134:   ! Demo of user-provided dense matrix storage (stack)
135:   call Demo3(m, n)

137:   PetscCallA(PetscFinalize(ierr))
138: end program main

140: !/*TEST
141: !
142: !   test:
143: !      nsize: 1
144: !      output_file: output/ex36f.out
145: !
146: !TEST*/