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);if (ierr) return ierr;
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();
43: if (total_subdomains > size) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of subdomains %D must not exceed comm size %D",total_subdomains,size);
44: if (total_subdomains < 1 || total_subdomains > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of subdomains must be > 0 and <= %D (comm size), got total_subdomains = %D",size,total_subdomains);
45: if (rep != 1 && rep != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid number of test repetitions: %D; must be 1 or 2",rep);
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=%D)",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 %D\n",gs);
131: ISView(rowis[ss],subviewer);
132: PetscViewerFlush(subviewer);
133: PetscViewerASCIIPrintf(subviewer,"Col IS %D\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 ierr;
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*/