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: !T*/
  9: !
 10: ! -----------------------------------------------------------------------

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

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

 34:       ptr => array
 35:       end subroutine

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

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

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

 57:       nullify(ptr)
 58:       end subroutine


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

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

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

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

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

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

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

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

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


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

144:       end

146: !/*TEST
147: !   build:
148: !     depends: ex21.c
149: !
150: !   test:
151: !
152: !TEST*/