Actual source code: ex11f90.F90
petsc-3.8.4 2018-03-24
1: program main
2: !-----------------------------------------------------------------------
3: !
4: ! Tests DMDAGetVecGetArray()
5: !-----------------------------------------------------------------------
6: !
8: #include <petsc/finclude/petscdm.h>
9: use petsc
10: implicit none
12: Type(tVec) g
13: Type(tDM) ada
15: PetscScalar,pointer :: x1(:),x2(:,:)
16: PetscScalar,pointer :: x3(:,:,:),x4(:,:,:,:)
17: PetscErrorCode ierr
18: PetscInt m,n,p,dof,s,i,j,k,xs,xl
19: PetscInt ys,yl
20: PetscInt zs,zl
22: m = 5
23: n = 6
24: p = 4;
25: s = 1
26: dof = 1
27: CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)
28: if (ierr .ne. 0) then
29: print*,'Unable to initialize PETSc'
30: stop
31: endif
32: call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m,dof,1,PETSC_NULL_INTEGER,ada,ierr);CHKERRA(ierr)
33: call DMSetUp(ada,ierr);CHKERRA(ierr)
34: call DMGetGlobalVector(ada,g,ierr);CHKERRA(ierr)
35: call DMDAGetCorners(ada,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xl,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
36: call DMDAVecGetArrayF90(ada,g,x1,ierr);CHKERRA(ierr)
37: do i=xs,xs+xl-1
38: ! CHKMEMQ
39: x1(i) = i
40: ! CHKMEMQ
41: enddo
42: call DMDAVecRestoreArrayF90(ada,g,x1,ierr);CHKERRA(ierr)
43: call VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
44: call DMRestoreGlobalVector(ada,g,ierr);CHKERRA(ierr)
45: call DMDestroy(ada,ierr);CHKERRA(ierr)
47: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,PETSC_DECIDE,PETSC_DECIDE,dof,s, &
48: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ada,ierr);CHKERRA(ierr)
49: call DMSetUp(ada,ierr);CHKERRA(ierr)
50: call DMGetGlobalVector(ada,g,ierr);CHKERRA(ierr)
51: call DMDAGetCorners(ada,xs,ys,PETSC_NULL_INTEGER,xl,yl,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
52: call DMDAVecGetArrayF90(ada,g,x2,ierr);CHKERRA(ierr)
53: do i=xs,xs+xl-1
54: do j=ys,ys+yl-1
55: ! CHKMEMQ
56: x2(i,j) = i + j
57: ! CHKMEMQ
58: enddo
59: enddo
60: call DMDAVecRestoreArrayF90(ada,g,x2,ierr);CHKERRA(ierr)
61: call VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
62: call DMRestoreGlobalVector(ada,g,ierr);CHKERRA(ierr)
63: call DMDestroy(ada,ierr);CHKERRA(ierr)
65: call DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, m,n,p,PETSC_DECIDE,PETSC_DECIDE, &
66: & PETSC_DECIDE,dof,s,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ada,ierr);CHKERRA(ierr)
67: call DMSetUp(ada,ierr);CHKERRA(ierr)
68: call DMGetGlobalVector(ada,g,ierr);CHKERRA(ierr)
69: call DMDAGetCorners(ada,xs,ys,zs,xl,yl,zl,ierr);CHKERRA(ierr)
70: call DMDAVecGetArrayF90(ada,g,x3,ierr);CHKERRA(ierr)
71: do i=xs,xs+xl-1
72: do j=ys,ys+yl-1
73: do k=zs,zs+zl-1
74: ! CHKMEMQ
75: x3(i,j,k) = i + j + k
76: ! CHKMEMQ
77: enddo
78: enddo
79: enddo
80: call DMDAVecRestoreArrayF90(ada,g,x3,ierr);CHKERRA(ierr)
81: call VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
82: call DMRestoreGlobalVector(ada,g,ierr);CHKERRA(ierr)
83: call DMDestroy(ada,ierr);CHKERRA(ierr)
85: !
86: ! Same tests but now with DOF > 1, so dimensions of array are one higher
87: !
88: dof = 2
89: call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m,dof,1,PETSC_NULL_INTEGER,ada,ierr);CHKERRA(ierr)
90: call DMSetUp(ada,ierr);CHKERRA(ierr)
91: call DMGetGlobalVector(ada,g,ierr);CHKERRA(ierr)
92: call DMDAGetCorners(ada,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xl,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
93: call DMDAVecGetArrayF90(ada,g,x2,ierr);CHKERRA(ierr)
94: do i=xs,xs+xl-1
95: ! CHKMEMQ
96: x2(0,i) = i
97: x2(1,i) = -i
98: ! CHKMEMQ
99: enddo
100: call DMDAVecRestoreArrayF90(ada,g,x1,ierr);CHKERRA(ierr)
101: call VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
102: call DMRestoreGlobalVector(ada,g,ierr);CHKERRA(ierr)
103: call DMDestroy(ada,ierr);CHKERRA(ierr)
105: dof = 2
106: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,PETSC_DECIDE,PETSC_DECIDE,dof,s, &
107: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ada,ierr);CHKERRA(ierr)
108: call DMSetUp(ada,ierr);CHKERRA(ierr)
109: call DMGetGlobalVector(ada,g,ierr);CHKERRA(ierr)
110: call DMDAGetCorners(ada,xs,ys,PETSC_NULL_INTEGER,xl,yl,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
111: call DMDAVecGetArrayF90(ada,g,x3,ierr);CHKERRA(ierr)
112: do i=xs,xs+xl-1
113: do j=ys,ys+yl-1
114: ! CHKMEMQ
115: x3(0,i,j) = i + j
116: x3(1,i,j) = -(i + j)
117: ! CHKMEMQ
118: enddo
119: enddo
120: call DMDAVecRestoreArrayF90(ada,g,x3,ierr);CHKERRA(ierr)
121: call VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
122: call DMRestoreGlobalVector(ada,g,ierr);CHKERRA(ierr)
123: call DMDestroy(ada,ierr);CHKERRA(ierr)
125: dof = 3
126: call DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,p,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,s, &
127: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ada,ierr);CHKERRA(ierr)
128: call DMSetUp(ada,ierr);CHKERRA(ierr)
129: call DMGetGlobalVector(ada,g,ierr);CHKERRA(ierr)
130: call DMDAGetCorners(ada,xs,ys,zs,xl,yl,zl,ierr);CHKERRA(ierr)
131: call DMDAVecGetArrayF90(ada,g,x4,ierr);CHKERRA(ierr)
132: do i=xs,xs+xl-1
133: do j=ys,ys+yl-1
134: do k=zs,zs+zl-1
135: ! CHKMEMQ
136: x4(0,i,j,k) = i + j + k
137: x4(1,i,j,k) = -(i + j + k)
138: x4(2,i,j,k) = i + j + k
139: ! CHKMEMQ
140: enddo
141: enddo
142: enddo
143: call DMDAVecRestoreArrayF90(ada,g,x4,ierr);CHKERRA(ierr)
144: call VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
145: call DMRestoreGlobalVector(ada,g,ierr);CHKERRA(ierr)
146: call DMDestroy(ada,ierr);CHKERRA(ierr)
148: CALL PetscFinalize(ierr)
149: stop
150: END PROGRAM