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: /*T
  6:    Concepts: Mat^mat submatrix, parallel
  7:    Processors: n
  8: T*/

 10: #include <petscmat.h>

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

 26:   PetscInitialize(&argc,&args,(char*)0,help);
 27:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 28:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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


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

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

 76:   PetscViewerASCIIPrintf(viewer,"Initial matrix:\n");
 77:   MatView(A,viewer);

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

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

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

212: /*TEST

214:    test:
215:       nsize: 2
216:       args: -total_subdomains 1
217:       output_file: output/ex183_2_1.out

219:    test:
220:       suffix: 2
221:       nsize: 3
222:       args: -total_subdomains 2
223:       output_file: output/ex183_3_2.out

225:    test:
226:       suffix: 3
227:       nsize: 4
228:       args: -total_subdomains 2
229:       output_file: output/ex183_4_2.out

231:    test:
232:       suffix: 4
233:       nsize: 6
234:       args: -total_subdomains 2
235:       output_file: output/ex183_6_2.out

237: TEST*/