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: !/*T
  6: !   Concepts: vectors^basic routines;
  7: !   Processors: n
  8: !   depends: ex21.c
  9: !T*/
 10: !
 11: ! -----------------------------------------------------------------------

 13:       module mymodule
 14:  #include <petsc/finclude/petscsys.h>
 15:       type MyStruct
 16:         sequence
 17:         PetscScalar :: a,b,c
 18:       end type MyStruct
 19:       end module

 21: !
 22: !  These routines are used internally by the C functions VecGetArrayMyStruct() and VecRestoreArrayMyStruct()
 23: !  Because Fortran requires "knowing" exactly what derived types the pointers to point too, these have to be
 24: !  customized for exactly the derived type in question
 25: !
 26:       subroutine F90Array1dCreateMyStruct(array,start,len,ptr)
 27:  #include <petsc/finclude/petscsys.h>
 28:       use petscsys
 29:       use mymodule
 30:       implicit none
 31:       PetscInt start,len
 32:       type(MyStruct), target :: array(start:start+len-1)
 33:       type(MyStruct), pointer :: ptr(:)

 35:       ptr => array
 36:       end subroutine

 38:       subroutine F90Array1dAccessMyStruct(ptr,address)
 39:  #include <petsc/finclude/petscsys.h>
 40:       use petscsys
 41:       use mymodule
 42:       implicit none
 43:       type(MyStruct), pointer :: ptr(:)
 44:       PetscFortranAddr address
 45:       PetscInt start

 47:       start = lbound(ptr,1)
 48:       call F90Array1dGetAddrMyStruct(ptr(start),address)
 49:       end subroutine

 51:       subroutine F90Array1dDestroyMyStruct(ptr)
 52:  #include <petsc/finclude/petscsys.h>
 53:       use petscsys
 54:       use mymodule
 55:       implicit none
 56:       type(MyStruct), pointer :: ptr(:)

 58:       nullify(ptr)
 59:       end subroutine


 62:       program main
 63:  #include <petsc/finclude/petscvec.h>
 64:       use petscvec
 65:       use mymodule
 66:       implicit none

 68: !
 69: !
 70: !   These two routines are defined in ex21.c they create the Fortran pointer to the derived type
 71: !
 72:       Interface
 73:         Subroutine VecGetArrayMyStruct(v,array,ierr)
 74:           use petscvec
 75:           use mymodule
 76:           type(MyStruct), pointer :: array(:)
 77:           PetscErrorCode ierr
 78:           Vec     v
 79:         End Subroutine
 80:       End Interface

 82:       Interface
 83:         Subroutine VecRestoreArrayMyStruct(v,array,ierr)
 84:           use petscvec
 85:           use mymodule
 86:           type(MyStruct), pointer :: array(:)
 87:           PetscErrorCode ierr
 88:           Vec     v
 89:         End Subroutine
 90:       End Interface

 92: !
 93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94: !                   Variable declarations
 95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96: !
 97: !  Variables:
 98: !     x, y, w - vectors
 99: !     z       - array of vectors
100: !
101:       Vec              x,y
102:       type(MyStruct),  pointer :: xarray(:)
103:       PetscInt         n
104:       PetscErrorCode   ierr
105:       PetscBool        flg
106:       integer          i

108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: !                 Beginning of program
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

112:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
113:       if (ierr .ne. 0) then
114:         print*,'PetscInitialize failed'
115:         stop
116:       endif
117:       n     = 30

119:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr);CHKERRA(ierr)
120:       call VecCreate(PETSC_COMM_WORLD,x,ierr);CHKERRA(ierr)
121:       call VecSetSizes(x,PETSC_DECIDE,n,ierr);CHKERRA(ierr)
122:       call VecSetFromOptions(x,ierr);CHKERRA(ierr)
123:       call VecDuplicate(x,y,ierr);CHKERRA(ierr)

125:       call VecGetArrayMyStruct(x,xarray,ierr);CHKERRA(ierr)
126:       do i=1,10
127:       xarray(i)%a = i
128:       xarray(i)%b = 100*i
129:       xarray(i)%c = 10000*i
130:       enddo

132:       call VecRestoreArrayMyStruct(x,xarray,ierr);CHKERRA(ierr)
133:       call VecView(x,PETSC_VIEWER_STDOUT_SELF,ierr);CHKERRA(ierr)
134:       call VecGetArrayMyStruct(x,xarray,ierr);CHKERRA(ierr)
135:       do i = 1 , 10
136:         write(*,*) abs(xarray(i)%a),abs(xarray(i)%b),abs(xarray(i)%c)
137:       end do
138:       call VecRestoreArrayMyStruct(x,xarray,ierr);CHKERRA(ierr)


141:       call VecDestroy(x,ierr);CHKERRA(ierr)
142:       call VecDestroy(y,ierr);CHKERRA(ierr)
143:       call PetscFinalize(ierr)

145:       end