1: !
  2: !
  3: !    Demonstrates how one may access entries of a PETSc Vec as if it was an array of Fortran derived types
  4: !
  5: !
  6: ! -----------------------------------------------------------------------
  7: #include <petsc/finclude/petscvec.h>
  8: module ex21module
  9:   use petscsys
 10:   implicit none
 11:   type MyStruct
 12:     sequence
 13:     PetscScalar :: a, b, c
 14:   end type MyStruct
 15: !
 16: !
 17: !   These two routines are defined in ex21.c they create the Fortran pointer to the derived type
 18: !
 19:   interface
 20:     subroutine VecGetArrayMyStruct(v, array, ierr)
 21:       use petscvec
 22:       import MyStruct
 23:       type(MyStruct), pointer :: array(:)
 24:       PetscErrorCode ierr
 25:       Vec v
 26:     end subroutine

 28:     subroutine VecRestoreArrayMyStruct(v, array, ierr)
 29:       use petscvec
 30:       import MyStruct
 31:       type(MyStruct), pointer :: array(:)
 32:       PetscErrorCode ierr
 33:       Vec v
 34:     end subroutine
 35:   end interface
 36: end module

 38: !  These routines are used internally by the C functions VecGetArrayMyStruct() and VecRestoreArrayMyStruct()
 39: !  Because Fortran requires "knowing" exactly what derived types the pointers to point too, these have to be
 40: !  customized for exactly the derived type in question
 41: !
 42: subroutine F90Array1dCreateMyStruct(array, start, len, ptr)
 43:   use petscsys
 44:   use ex21module
 45:   implicit none
 46:   PetscInt start, len
 47:   type(MyStruct), target :: array(start:start + len - 1)
 48:   type(MyStruct), pointer :: ptr(:)

 50:   ptr => array
 51: end subroutine

 53: subroutine F90Array1dAccessMyStruct(ptr, address)
 54:   use petscsys
 55:   use ex21module
 56:   implicit none
 57:   type(MyStruct), pointer :: ptr(:)
 58:   PetscFortranAddr address
 59:   PetscInt start

 61:   start = lbound(ptr, 1)
 62:   call F90Array1dGetAddrMyStruct(ptr(start), address)
 63: end subroutine

 65: subroutine F90Array1dDestroyMyStruct(ptr)
 66:   use ex21module
 67:   implicit none
 68:   type(MyStruct), pointer :: ptr(:)

 70:   nullify (ptr)
 71: end subroutine

 73: program main
 74:   use petscvec
 75:   use ex21module
 76:   implicit none

 78: !
 79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 80: !                   Variable declarations
 81: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 82: !
 83: !  Variables:
 84: !     x, y, w - vectors
 85: !     z       - array of vectors
 86: !
 87:   Vec x, y
 88:   type(MyStruct), pointer :: xarray(:)
 89:   PetscInt n
 90:   PetscErrorCode ierr
 91:   PetscBool flg
 92:   integer i

 94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95: !                 Beginning of program
 96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 98:   PetscCallA(PetscInitialize(ierr))

100:   n = 30
101:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
102:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
103:   PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr))
104:   PetscCallA(VecSetFromOptions(x, ierr))
105:   PetscCallA(VecDuplicate(x, y, ierr))

107:   PetscCallA(VecGetArrayMyStruct(x, xarray, ierr))
108:   do i = 1, 10
109:     xarray(i)%a = i
110:     xarray(i)%b = 100*i
111:     xarray(i)%c = 10000*i
112:   end do

114:   PetscCallA(VecRestoreArrayMyStruct(x, xarray, ierr))
115:   PetscCallA(VecView(x, PETSC_VIEWER_STDOUT_SELF, ierr))
116:   PetscCallA(VecGetArrayMyStruct(x, xarray, ierr))
117:   do i = 1, 10
118:     write (*, *) abs(xarray(i)%a), abs(xarray(i)%b), abs(xarray(i)%c)
119:   end do
120:   PetscCallA(VecRestoreArrayMyStruct(x, xarray, ierr))

122:   PetscCallA(VecDestroy(x, ierr))
123:   PetscCallA(VecDestroy(y, ierr))
124:   PetscCallA(PetscFinalize(ierr))

126: end

128: !/*TEST
129: !
130: !   build:
131: !     depends: ex21.c
132: !
133: !   test:
134: !
135: !TEST*/