Actual source code: ex21f.F

  1: !
  2: !
  3:       program main
 4:  #include finclude/petscsys.h
 5:  #include finclude/petscvec.h
 6:  #include finclude/petscda.h

  8: !         For testing purposes this example also creates a
  9: !   DA context. Actually codes using SDA routines will probably
 10: !   not also work with DA contexts.


 13:       integer        MM,ierr,dof,stencil_width,flg,i,start,end,PP
 14:       integer        flg2,flg3,NN,m,n,p
 15:       PetscOffset    in_idx,out_idx
 16:       DAPeriodicType periodic
 17:       DAStencilType  stencil_type
 18:       DA             da
 19:       integer        sda
 20:       Vec            local,global
 21:       Vec            local_copy
 22:       PetscScalar    value,mone,in(1),out(1)
 23:       PetscScalar    norm,work
 24: 
 25:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 27:       m             = PETSC_DECIDE
 28:       n             = PETSC_DECIDE
 29:       p             = PETSC_DECIDE
 30:       MM            = 8
 31:       NN            = 6
 32:       PP            = 5
 33:       dof           = 1
 34:       stencil_width = 1
 35:       periodic      = DA_NONPERIODIC
 36:       stencil_type  = DA_STENCIL_STAR


 39:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-M',MM,flg,ierr)
 40:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-N',NN,flg,ierr)
 41:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dof',dof,flg,ierr)
 42:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-stencil_width',         &
 43:      &     stencil_width,flg,ierr)
 44:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-periodic',periodic,     &
 45:      &     flg,ierr)
 46:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-stencil_type',          &
 47:      &     stencil_type,flg,ierr)

 49:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-2d',flg2,ierr)
 50:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-3d',flg3,ierr)
 51:       if (flg2 .ne. 0) then
 52:          call DACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,        &
 53:      &        MM,NN,m,n,dof,stencil_width,PETSC_NULL_INTEGER,           &
 54:      &        PETSC_NULL_INTEGER,da,ierr)
 55:          call SDACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,       &
 56:      &        MM,NN,m,n,dof,stencil_width,PETSC_NULL_INTEGER,           &
 57:      &        PETSC_NULL_INTEGER,sda,ierr)
 58:       else if (flg3 .ne. 0) then
 59:          call DACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,        &
 60:      &        MM,NN,PP,m,n,p,dof,stencil_width,PETSC_NULL_INTEGER,      &
 61:      &        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
 62:          call SDACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,       &
 63:      &        MM,NN,PP,m,n,p,dof,stencil_width,PETSC_NULL_INTEGER,      &
 64:      &        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,sda,ierr)
 65:       else
 66:          call DACreate1d(PETSC_COMM_WORLD,periodic,MM,dof,              &
 67:      &        stencil_width,PETSC_NULL,da,ierr)
 68:          call SDACreate1d(PETSC_COMM_WORLD,periodic,MM,dof,             &
 69:      &        stencil_width,PETSC_NULL_INTEGER,sda,ierr)
 70:       endif

 72:       call DACreateGlobalVector(da,global,ierr)
 73:       call DACreateLocalVector(da,local,ierr)
 74:       call VecDuplicate(local,local_copy,ierr)

 76: 
 77: !   zero out vectors so that ghostpoints are zero
 78:       value = 0.0
 79:       call VecSet(local,value,ierr)
 80:       call VecSet(local_copy,value,ierr)

 82:       call VecGetOwnershipRange(global,start,end,ierr)
 83:       do 10, i=start,end-1
 84:         value = i + 1
 85:         call VecSetValues(global,1,i,value,INSERT_VALUES,ierr)
 86:  10   continue
 87:       call VecAssemblyBegin(global,ierr)
 88:       call VecAssemblyEnd(global,ierr)

 90:       call DAGlobalToLocalBegin(da,global,INSERT_VALUES,local,          &
 91:      &                          ierr)
 92:       call DAGlobalToLocalEnd(da,global,INSERT_VALUES,local,ierr)


 95:       call VecGetArray(local,out,out_idx,ierr)
 96:       call VecGetArray(local_copy,in,in_idx,ierr)
 97:       call SDALocalToLocalBegin(sda,out(out_idx+1),INSERT_VALUES,       &
 98:      &                          in(in_idx+1),ierr)
 99:       call SDALocalToLocalEnd(sda,out(out_idx+1),INSERT_VALUES,         &
100:      &                        in(in_idx+1),ierr)

102:       mone = -1.0
103:       call VecAXPY(local_copy,mone,local,ierr)
104:       call VecNorm(local_copy,NORM_MAX,work,ierr)
105:       call MPI_Allreduce(work, norm,1,MPIU_REAL_PRECISION,MPI_MAX,      &
106:      &                   PETSC_COMM_WORLD,ierr)
107:       print*,'Norm of difference ',norm,' should be zero'
108: 
109:       call DADestroy(da,ierr)
110:       call SDADestroy(sda,ierr)
111:       call VecDestroy(local_copy,ierr)

113:       call VecDestroy(local,ierr)
114:       call VecDestroy(global,ierr)
115:       call PetscFinalize(ierr)
116:       end