xref: /petsc/src/vec/is/sf/tests/ex15.c (revision 8dbb0df66754dee4fb72dff2ad56e76db1e6b7c7)
1 static char help[]= "  Test VecScatterRemap() on various vecscatter. \n\
2 We may do optimization based on index patterns. After index remapping by VecScatterRemap(), we need to \n\
3 make sure the vecscatter works as expected with the optimizaiton. \n\
4   VecScatterRemap() does not support all kinds of vecscatters. In addition, it only supports remapping \n\
5 entries where we read the data (i.e., todata in paralle scatter, fromdata in sequential scatter). This test \n\
6 tests VecScatterRemap on parallel to paralle (PtoP) vecscatter, sequential general to sequential \n\
7 general (SGToSG) vecscatter and sequential general to sequential stride 1 (SGToSS_Stride1) vecscatter.\n\n";
8 
9 #include <petscvec.h>
10 
11 int main(int argc,char **argv)
12 {
13   PetscErrorCode     ierr;
14   PetscInt           i,n,*ix,*iy,*tomap,start;
15   Vec                x,y;
16   PetscMPIInt        nproc,rank;
17   IS                 isx,isy;
18   const PetscInt     *ranges;
19   VecScatter         vscat;
20 
21   PetscFunctionBegin;
22   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
23   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&nproc);CHKERRMPI(ierr);
24   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
25 
26   if (nproc != 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This test must run with exactly two MPI ranks\n");
27 
28   /* ====================================================================
29      (1) test VecScatterRemap on a parallel to parallel (PtoP) vecscatter
30      ====================================================================
31    */
32 
33   n = 64;  /* long enough to trigger memcpy optimizations both in local scatter and remote scatter */
34 
35   /* create two MPI vectors x, y of length n=64, N=128 */
36   ierr = VecCreateMPI(PETSC_COMM_WORLD,n,PETSC_DECIDE,&x);CHKERRQ(ierr);
37   ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
38 
39   /* Initialize x as {0~127} */
40   ierr = VecGetOwnershipRanges(x,&ranges);CHKERRQ(ierr);
41   for (i=ranges[rank]; i<ranges[rank+1]; i++) { ierr = VecSetValue(x,i,(PetscScalar)i,INSERT_VALUES);CHKERRQ(ierr); }
42   ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
43   ierr = VecAssemblyEnd(x);CHKERRQ(ierr);
44 
45   /* create two general index sets isx = {0~127} and isy = {32~63,64~95,96~127,0~31}. isx is sequential, but we use
46      it as general and let PETSc detect the pattern and optimize it. indices in isy are set to make the vecscatter
47      have both local scatter and remote scatter (i.e., MPI communication)
48    */
49   ierr = PetscMalloc2(n,&ix,n,&iy);CHKERRQ(ierr);
50   start = ranges[rank];
51   for (i=ranges[rank]; i<ranges[rank+1]; i++) ix[i-start] = i;
52   ierr = ISCreateGeneral(PETSC_COMM_WORLD,n,ix,PETSC_COPY_VALUES,&isx);CHKERRQ(ierr);
53 
54   if (rank == 0) { for (i=0; i<n; i++) iy[i] = i+32; }
55   else for (i=0; i<n/2; i++) { iy[i] = i+96; iy[i+n/2] = i; }
56 
57   ierr = ISCreateGeneral(PETSC_COMM_WORLD,n,iy,PETSC_COPY_VALUES,&isy);CHKERRQ(ierr);
58 
59   /* create a vecscatter that shifts x to the tail by quater periodically and puts the results in y */
60   ierr = VecScatterCreate(x,isx,y,isy,&vscat);CHKERRQ(ierr);
61   ierr = VecScatterBegin(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
62   ierr = VecScatterEnd(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
63 
64   /* view y to check the result. y should be {Q3,Q0,Q1,Q2} of x, that is {96~127,0~31,32~63,64~95} */
65   ierr = PetscPrintf(PETSC_COMM_WORLD,"Before VecScatterRemap on PtoP, MPI vector y is:\n");CHKERRQ(ierr);
66   ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
67 
68   /* now call the weird subroutine VecScatterRemap to slightly change the vecscatter. It changes where we read vector
69      x entries to send out, but does not change the communication pattern (i.e., send/recv pairs and msg lengths).
70 
71      We create tomap as {32~63,0~31}. Originaly, we read from indices {0~64} of the local x to send out. The remap
72      does indices[i] = tomap[indices[i]]. Therefore, after the remap, we read from indices {32~63,0~31} of the local x.
73      isy is unchanged. So, we will shift x to {Q2,Q1,Q0,Q3}, that is {64~95,32~63,0~31,96~127}
74   */
75   ierr = PetscMalloc1(n,&tomap);CHKERRQ(ierr);
76   for (i=0; i<n/2; i++) { tomap[i] = i+n/2; tomap[i+n/2] = i; };
77   ierr = VecScatterRemap(vscat,tomap,NULL);CHKERRQ(ierr);
78   ierr = VecScatterBegin(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
79   ierr = VecScatterEnd(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
80 
81   /* view y to check the result. y should be {64~95,32~63,0~31,96~127} */
82   ierr = PetscPrintf(PETSC_COMM_WORLD,"After VecScatterRemap on PtoP, MPI vector y is:\n");CHKERRQ(ierr);
83   ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
84 
85   /* destroy everything before we recreate them in different types */
86   ierr = PetscFree2(ix,iy);CHKERRQ(ierr);
87   ierr = VecDestroy(&x);CHKERRQ(ierr);
88   ierr = VecDestroy(&y);CHKERRQ(ierr);
89   ierr = ISDestroy(&isx);CHKERRQ(ierr);
90   ierr = ISDestroy(&isy);CHKERRQ(ierr);
91   ierr = PetscFree(tomap);CHKERRQ(ierr);
92   ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr);
93 
94   /* ==========================================================================================
95      (2) test VecScatterRemap on a sequential general to sequential general (SGToSG) vecscatter
96      ==========================================================================================
97    */
98   n = 64; /* long enough to trigger memcpy optimizations in local scatter */
99 
100   /* create two seq vectors x, y of length n */
101   ierr = VecCreateSeq(PETSC_COMM_SELF,n,&x);CHKERRQ(ierr);
102   ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
103 
104   /* Initialize x as {0~63} */
105   for (i=0; i<n; i++) { ierr = VecSetValue(x,i,(PetscScalar)i,INSERT_VALUES);CHKERRQ(ierr); }
106   ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
107   ierr = VecAssemblyEnd(x);CHKERRQ(ierr);
108 
109   /* create two general index sets isx = isy = {0~63}, which are sequential, but we use them as
110      general and let PETSc detect the pattern and optimize it */
111   ierr = PetscMalloc2(n,&ix,n,&iy);CHKERRQ(ierr);
112   for (i=0; i<n; i++) ix[i] = i;
113   ierr = ISCreateGeneral(PETSC_COMM_SELF,n,ix,PETSC_COPY_VALUES,&isx);CHKERRQ(ierr);
114   ierr = ISDuplicate(isx,&isy);CHKERRQ(ierr);
115 
116   /* create a vecscatter that just copies x to y */
117   ierr = VecScatterCreate(x,isx,y,isy,&vscat);CHKERRQ(ierr);
118   ierr = VecScatterBegin(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
119   ierr = VecScatterEnd(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
120 
121   /* view y to check the result. y should be {0~63} */
122   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nBefore VecScatterRemap on SGToSG, SEQ vector y is:\n");CHKERRQ(ierr);
123   ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
124 
125   /* now call the weird subroutine VecScatterRemap to slightly change the vecscatter.
126 
127      Create tomap as {32~63,0~31}. Originaly, we read from indices {0~64} of seq x to write to y. The remap
128      does indices[i] = tomap[indices[i]]. Therefore, after the remap, we read from indices{32~63,0~31} of seq x.
129    */
130   ierr = PetscMalloc1(n,&tomap);CHKERRQ(ierr);
131   for (i=0; i<n/2; i++) { tomap[i] = i+n/2; tomap[i+n/2] = i; };
132   ierr = VecScatterRemap(vscat,tomap,NULL);CHKERRQ(ierr);
133   ierr = VecScatterBegin(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
134   ierr = VecScatterEnd(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
135 
136   /* view y to check the result. y should be {32~63,0~31} */
137   ierr = PetscPrintf(PETSC_COMM_WORLD,"After VecScatterRemap on SGToSG, SEQ vector y is:\n");CHKERRQ(ierr);
138   ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
139 
140   /* destroy everything before we recreate them in different types */
141   ierr = PetscFree2(ix,iy);CHKERRQ(ierr);
142   ierr = VecDestroy(&x);CHKERRQ(ierr);
143   ierr = VecDestroy(&y);CHKERRQ(ierr);
144   ierr = ISDestroy(&isx);CHKERRQ(ierr);
145   ierr = ISDestroy(&isy);CHKERRQ(ierr);
146   ierr = PetscFree(tomap);CHKERRQ(ierr);
147   ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr);
148 
149   /* ===================================================================================================
150      (3) test VecScatterRemap on a sequential general to sequential stride 1 (SGToSS_Stride1) vecscatter
151      ===================================================================================================
152    */
153   n = 64; /* long enough to trigger memcpy optimizations in local scatter */
154 
155   /* create two seq vectors x of length n, and y of length n/2 */
156   ierr = VecCreateSeq(PETSC_COMM_SELF,n,&x);CHKERRQ(ierr);
157   ierr = VecCreateSeq(PETSC_COMM_SELF,n/2,&y);CHKERRQ(ierr);
158 
159   /* Initialize x as {0~63} */
160   for (i=0; i<n; i++) { ierr = VecSetValue(x,i,(PetscScalar)i,INSERT_VALUES);CHKERRQ(ierr); }
161   ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
162   ierr = VecAssemblyEnd(x);CHKERRQ(ierr);
163 
164   /* create a general index set isx = {0:63:2}, which actually is a stride IS with first=0, n=32, step=2,
165      but we use it as general and let PETSc detect the pattern and optimize it. */
166   ierr = PetscMalloc2(n/2,&ix,n/2,&iy);CHKERRQ(ierr);
167   for (i=0; i<n/2; i++) ix[i] = i*2;
168   ierr = ISCreateGeneral(PETSC_COMM_SELF,n/2,ix,PETSC_COPY_VALUES,&isx);CHKERRQ(ierr);
169 
170   /* create a stride1 index set isy = {0~31}. We intentionally set the step to 1 to trigger optimizations */
171   ierr = ISCreateStride(PETSC_COMM_SELF,32,0,1,&isy);CHKERRQ(ierr);
172 
173   /* create a vecscatter that just copies even entries of x to y */
174   ierr = VecScatterCreate(x,isx,y,isy,&vscat);CHKERRQ(ierr);
175   ierr = VecScatterBegin(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
176   ierr = VecScatterEnd(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
177 
178   /* view y to check the result. y should be {0:63:2} */
179   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nBefore VecScatterRemap on SGToSS_Stride1, SEQ vector y is:\n");CHKERRQ(ierr);
180   ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
181 
182   /* now call the weird subroutine VecScatterRemap to slightly change the vecscatter.
183 
184      Create tomap as {32~63,0~31}. Originaly, we read from indices{0:63:2} of seq x to write to y. The remap
185      does indices[i] = tomap[indices[i]]. Therefore, after the remap, we read from indices{32:63:2,0:31:2} of seq x.
186    */
187   ierr = PetscMalloc1(n,&tomap);CHKERRQ(ierr);
188   for (i=0; i<n/2; i++) { tomap[i] = i+n/2; tomap[i+n/2] = i; };
189   ierr = VecScatterRemap(vscat,tomap,NULL);CHKERRQ(ierr);
190   ierr = VecScatterBegin(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
191   ierr = VecScatterEnd(vscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
192 
193   /* view y to check the result. y should be {32:63:2,0:31:2} */
194   ierr = PetscPrintf(PETSC_COMM_WORLD,"After VecScatterRemap on SGToSS_Stride1, SEQ vector y is:\n");CHKERRQ(ierr);
195   ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
196 
197   /* destroy everything before PetscFinalize */
198   ierr = PetscFree2(ix,iy);CHKERRQ(ierr);
199   ierr = VecDestroy(&x);CHKERRQ(ierr);
200   ierr = VecDestroy(&y);CHKERRQ(ierr);
201   ierr = ISDestroy(&isx);CHKERRQ(ierr);
202   ierr = ISDestroy(&isy);CHKERRQ(ierr);
203   ierr = PetscFree(tomap);CHKERRQ(ierr);
204   ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr);
205 
206   ierr = PetscFinalize();
207   return ierr;
208 }
209 
210 /*TEST
211 
212    test:
213       suffix: 1
214       nsize: 2
215       diff_args: -j
216       requires: double
217 TEST*/
218 
219