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*/