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