Actual source code: ex183.c

  1: static char help[] = "Example of extracting an array of MPI submatrices from a given MPI matrix.\n"
  2:   "This test can only be run in parallel.\n"
  3:   "\n";

  5: #include <petscmat.h>

  7: int main(int argc, char **args)
  8: {
  9:   Mat             A,*submats;
 10:   MPI_Comm        subcomm;
 11:   PetscMPIInt     rank,size,subrank,subsize,color;
 12:   PetscInt        m,n,N,bs,rstart,rend,i,j,k,total_subdomains,hash,nsubdomains=1;
 13:   PetscInt        nis,*cols,gnsubdomains,gsubdomainnums[1],gsubdomainperm[1],s,gs;
 14:   PetscInt        *rowindices,*colindices,idx,rep;
 15:   PetscScalar     *vals;
 16:   IS              rowis[1],colis[1];
 17:   PetscViewer     viewer;
 18:   PetscBool       permute_indices,flg;
 19:   PetscErrorCode  ierr;

 21:   PetscInitialize(&argc,&args,(char*)0,help);
 22:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 23:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 25:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex183","Mat");
 26:   m = 5;
 27:   PetscOptionsInt("-m","Local matrix size","MatSetSizes",m,&m,&flg);
 28:   total_subdomains = size-1;
 29:   PetscOptionsInt("-total_subdomains","Number of submatrices where 0 < n < comm size","MatCreateSubMatricesMPI",total_subdomains,&total_subdomains,&flg);
 30:   permute_indices = PETSC_FALSE;
 31:   PetscOptionsBool("-permute_indices","Whether to permute indices before breaking them into subdomains","ISCreateGeneral",permute_indices,&permute_indices,&flg);
 32:   hash = 7;
 33:   PetscOptionsInt("-hash","Permutation factor, which has to be relatively prime to M = size*m (total matrix size)","ISCreateGeneral",hash,&hash,&flg);
 34:   rep = 2;
 35:   PetscOptionsInt("-rep","Number of times to carry out submatrix extractions; currently only 1 & 2 are supported",NULL,rep,&rep,&flg);
 36:   PetscOptionsEnd();


 42:   viewer = PETSC_VIEWER_STDOUT_WORLD;
 43:   /* Create logically sparse, but effectively dense matrix for easy verification of submatrix extraction correctness. */
 44:   MatCreate(PETSC_COMM_WORLD,&A);
 45:   MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE);
 46:   MatSetFromOptions(A);
 47:   MatSetUp(A);
 48:   MatGetSize(A,NULL,&N);
 49:   MatGetLocalSize(A,NULL,&n);
 50:   MatGetBlockSize(A,&bs);
 51:   MatSeqAIJSetPreallocation(A,n,NULL);
 52:   MatMPIAIJSetPreallocation(A,n,NULL,N-n,NULL);
 53:   MatSeqBAIJSetPreallocation(A,bs,n/bs,NULL);
 54:   MatMPIBAIJSetPreallocation(A,bs,n/bs,NULL,(N-n)/bs,NULL);
 55:   MatSeqSBAIJSetPreallocation(A,bs,n/bs,NULL);
 56:   MatMPISBAIJSetPreallocation(A,bs,n/bs,NULL,(N-n)/bs,NULL);

 58:   PetscMalloc2(N,&cols,N,&vals);
 59:   MatGetOwnershipRange(A,&rstart,&rend);
 60:   for (j = 0; j < N; ++j) cols[j] = j;
 61:   for (i=rstart; i<rend; i++) {
 62:     for (j=0;j<N;++j) {
 63:       vals[j] = i*10000+j;
 64:     }
 65:     MatSetValues(A,1,&i,N,cols,vals,INSERT_VALUES);
 66:   }
 67:   PetscFree2(cols,vals);
 68:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 71:   PetscViewerASCIIPrintf(viewer,"Initial matrix:\n");
 72:   MatView(A,viewer);

 74:   /*
 75:      Create subcomms and ISs so that each rank participates in one IS.
 76:      The IS either coalesces adjacent rank indices (contiguous),
 77:      or selects indices by scrambling them using a hash.
 78:   */
 79:   k = size/total_subdomains + (size%total_subdomains>0); /* There are up to k ranks to a color */
 80:   color = rank/k;
 81:   MPI_Comm_split(PETSC_COMM_WORLD,color,rank,&subcomm);
 82:   MPI_Comm_size(subcomm,&subsize);
 83:   MPI_Comm_rank(subcomm,&subrank);
 84:   MatGetOwnershipRange(A,&rstart,&rend);
 85:   nis = 1;
 86:   PetscMalloc2(rend-rstart,&rowindices,rend-rstart,&colindices);

 88:   for (j = rstart; j < rend; ++j) {
 89:     if (permute_indices) {
 90:       idx = (j*hash);
 91:     } else {
 92:       idx = j;
 93:     }
 94:     rowindices[j-rstart] = idx%N;
 95:     colindices[j-rstart] = (idx+m)%N;
 96:   }
 97:   ISCreateGeneral(subcomm,rend-rstart,rowindices,PETSC_COPY_VALUES,&rowis[0]);
 98:   ISCreateGeneral(subcomm,rend-rstart,colindices,PETSC_COPY_VALUES,&colis[0]);
 99:   ISSort(rowis[0]);
100:   ISSort(colis[0]);
101:   PetscFree2(rowindices,colindices);
102:   /*
103:     Now view the ISs.  To avoid deadlock when viewing a list of objects on different subcomms,
104:     we need to obtain the global numbers of our local objects and wait for the corresponding global
105:     number to be viewed.
106:   */
107:   PetscViewerASCIIPrintf(viewer,"Subdomains");
108:   if (permute_indices) {
109:     PetscViewerASCIIPrintf(viewer," (hash=%" PetscInt_FMT ")",hash);
110:   }
111:   PetscViewerASCIIPrintf(viewer,":\n");
112:   PetscViewerFlush(viewer);

114:   nsubdomains = 1;
115:   for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s;
116:   PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)rowis,&gnsubdomains,gsubdomainnums);
117:   PetscSortIntWithPermutation(nsubdomains,gsubdomainnums,gsubdomainperm);
118:   for (gs=0,s=0; gs < gnsubdomains;++gs) {
119:     if (s < nsubdomains) {
120:       PetscInt ss;
121:       ss = gsubdomainperm[s];
122:       if (gs == gsubdomainnums[ss]) { /* Global subdomain gs being viewed is my subdomain with local number ss. */
123:         PetscViewer subviewer = NULL;
124:         PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)rowis[ss]),&subviewer);
125:         PetscViewerASCIIPrintf(subviewer,"Row IS %" PetscInt_FMT "\n",gs);
126:         ISView(rowis[ss],subviewer);
127:         PetscViewerFlush(subviewer);
128:         PetscViewerASCIIPrintf(subviewer,"Col IS %" PetscInt_FMT "\n",gs);
129:         ISView(colis[ss],subviewer);
130:         PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)rowis[ss]),&subviewer);
131:         ++s;
132:       }
133:     }
134:     MPI_Barrier(PETSC_COMM_WORLD);
135:   }
136:   PetscViewerFlush(viewer);
137:   ISSort(rowis[0]);
138:   ISSort(colis[0]);
139:   nsubdomains = 1;
140:   MatCreateSubMatricesMPI(A,nsubdomains,rowis,colis,MAT_INITIAL_MATRIX,&submats);
141:   /*
142:     Now view the matrices.  To avoid deadlock when viewing a list of objects on different subcomms,
143:     we need to obtain the global numbers of our local objects and wait for the corresponding global
144:     number to be viewed.
145:   */
146:   PetscViewerASCIIPrintf(viewer,"Submatrices (repetition 1):\n");
147:   for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s;
148:   PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)submats,&gnsubdomains,gsubdomainnums);
149:   PetscSortIntWithPermutation(nsubdomains,gsubdomainnums,gsubdomainperm);
150:   for (gs=0,s=0; gs < gnsubdomains;++gs) {
151:     if (s < nsubdomains) {
152:       PetscInt ss;
153:       ss = gsubdomainperm[s];
154:       if (gs == gsubdomainnums[ss]) { /* Global subdomain gs being viewed is my subdomain with local number ss. */
155:         PetscViewer subviewer = NULL;
156:         PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer);
157:         MatView(submats[ss],subviewer);
158:         PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer);
159:         ++s;
160:       }
161:     }
162:     MPI_Barrier(PETSC_COMM_WORLD);
163:   }
164:   PetscViewerFlush(viewer);
165:   if (rep == 1) goto cleanup;
166:   nsubdomains = 1;
167:   MatCreateSubMatricesMPI(A,nsubdomains,rowis,colis,MAT_REUSE_MATRIX,&submats);
168:   /*
169:     Now view the matrices.  To avoid deadlock when viewing a list of objects on different subcomms,
170:     we need to obtain the global numbers of our local objects and wait for the corresponding global
171:     number to be viewed.
172:   */
173:   PetscViewerASCIIPrintf(viewer,"Submatrices (repetition 2):\n");
174:   for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s;
175:   PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)submats,&gnsubdomains,gsubdomainnums);
176:   PetscSortIntWithPermutation(nsubdomains,gsubdomainnums,gsubdomainperm);
177:   for (gs=0,s=0; gs < gnsubdomains;++gs) {
178:     if (s < nsubdomains) {
179:       PetscInt ss;
180:       ss = gsubdomainperm[s];
181:       if (gs == gsubdomainnums[ss]) { /* Global subdomain gs being viewed is my subdomain with local number ss. */
182:         PetscViewer subviewer = NULL;
183:         PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer);
184:         MatView(submats[ss],subviewer);
185:         PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer);
186:         ++s;
187:       }
188:     }
189:     MPI_Barrier(PETSC_COMM_WORLD);
190:   }
191:   PetscViewerFlush(viewer);
192:   cleanup:
193:   for (k=0;k<nsubdomains;++k) {
194:     MatDestroy(submats+k);
195:   }
196:   PetscFree(submats);
197:   for (k=0;k<nis;++k) {
198:     ISDestroy(rowis+k);
199:     ISDestroy(colis+k);
200:   }
201:   MatDestroy(&A);
202:   MPI_Comm_free(&subcomm);
203:   PetscFinalize();
204:   return 0;
205: }

207: /*TEST

209:    test:
210:       nsize: 2
211:       args: -total_subdomains 1
212:       output_file: output/ex183_2_1.out

214:    test:
215:       suffix: 2
216:       nsize: 3
217:       args: -total_subdomains 2
218:       output_file: output/ex183_3_2.out

220:    test:
221:       suffix: 3
222:       nsize: 4
223:       args: -total_subdomains 2
224:       output_file: output/ex183_4_2.out

226:    test:
227:       suffix: 4
228:       nsize: 6
229:       args: -total_subdomains 2
230:       output_file: output/ex183_6_2.out

232: TEST*/