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"; 4 5 /*T 6 Concepts: Mat^mat submatrix, parallel 7 Processors: n 8 T*/ 9 10 #include <petscmat.h> 11 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; 25 26 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 27 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 28 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 29 30 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex183","Mat");PetscCall(ierr); 31 m = 5; 32 PetscCall(PetscOptionsInt("-m","Local matrix size","MatSetSizes",m,&m,&flg)); 33 total_subdomains = size-1; 34 PetscCall(PetscOptionsInt("-total_subdomains","Number of submatrices where 0 < n < comm size","MatCreateSubMatricesMPI",total_subdomains,&total_subdomains,&flg)); 35 permute_indices = PETSC_FALSE; 36 PetscCall(PetscOptionsBool("-permute_indices","Whether to permute indices before breaking them into subdomains","ISCreateGeneral",permute_indices,&permute_indices,&flg)); 37 hash = 7; 38 PetscCall(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 PetscCall(PetscOptionsInt("-rep","Number of times to carry out submatrix extractions; currently only 1 & 2 are supported",NULL,rep,&rep,&flg)); 41 ierr = PetscOptionsEnd();PetscCall(ierr); 42 43 PetscCheckFalse(total_subdomains > size,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of subdomains %" PetscInt_FMT " must not exceed comm size %d",total_subdomains,size); 44 PetscCheckFalse(total_subdomains < 1 || total_subdomains > size,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of subdomains must be > 0 and <= %d (comm size), got total_subdomains = %" PetscInt_FMT,size,total_subdomains); 45 PetscCheckFalse(rep != 1 && rep != 2,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid number of test repetitions: %" PetscInt_FMT "; must be 1 or 2",rep); 46 47 viewer = PETSC_VIEWER_STDOUT_WORLD; 48 /* Create logically sparse, but effectively dense matrix for easy verification of submatrix extraction correctness. */ 49 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 50 PetscCall(MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE)); 51 PetscCall(MatSetFromOptions(A)); 52 PetscCall(MatSetUp(A)); 53 PetscCall(MatGetSize(A,NULL,&N)); 54 PetscCall(MatGetLocalSize(A,NULL,&n)); 55 PetscCall(MatGetBlockSize(A,&bs)); 56 PetscCall(MatSeqAIJSetPreallocation(A,n,NULL)); 57 PetscCall(MatMPIAIJSetPreallocation(A,n,NULL,N-n,NULL)); 58 PetscCall(MatSeqBAIJSetPreallocation(A,bs,n/bs,NULL)); 59 PetscCall(MatMPIBAIJSetPreallocation(A,bs,n/bs,NULL,(N-n)/bs,NULL)); 60 PetscCall(MatSeqSBAIJSetPreallocation(A,bs,n/bs,NULL)); 61 PetscCall(MatMPISBAIJSetPreallocation(A,bs,n/bs,NULL,(N-n)/bs,NULL)); 62 63 PetscCall(PetscMalloc2(N,&cols,N,&vals)); 64 PetscCall(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 PetscCall(MatSetValues(A,1,&i,N,cols,vals,INSERT_VALUES)); 71 } 72 PetscCall(PetscFree2(cols,vals)); 73 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 74 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 75 76 PetscCall(PetscViewerASCIIPrintf(viewer,"Initial matrix:\n")); 77 PetscCall(MatView(A,viewer)); 78 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 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD,color,rank,&subcomm)); 87 PetscCallMPI(MPI_Comm_size(subcomm,&subsize)); 88 PetscCallMPI(MPI_Comm_rank(subcomm,&subrank)); 89 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 90 nis = 1; 91 PetscCall(PetscMalloc2(rend-rstart,&rowindices,rend-rstart,&colindices)); 92 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 PetscCall(ISCreateGeneral(subcomm,rend-rstart,rowindices,PETSC_COPY_VALUES,&rowis[0])); 103 PetscCall(ISCreateGeneral(subcomm,rend-rstart,colindices,PETSC_COPY_VALUES,&colis[0])); 104 PetscCall(ISSort(rowis[0])); 105 PetscCall(ISSort(colis[0])); 106 PetscCall(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 PetscCall(PetscViewerASCIIPrintf(viewer,"Subdomains")); 113 if (permute_indices) { 114 PetscCall(PetscViewerASCIIPrintf(viewer," (hash=%" PetscInt_FMT ")",hash)); 115 } 116 PetscCall(PetscViewerASCIIPrintf(viewer,":\n")); 117 PetscCall(PetscViewerFlush(viewer)); 118 119 nsubdomains = 1; 120 for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s; 121 PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)rowis,&gnsubdomains,gsubdomainnums)); 122 PetscCall(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 PetscCall(PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)rowis[ss]),&subviewer)); 130 PetscCall(PetscViewerASCIIPrintf(subviewer,"Row IS %" PetscInt_FMT "\n",gs)); 131 PetscCall(ISView(rowis[ss],subviewer)); 132 PetscCall(PetscViewerFlush(subviewer)); 133 PetscCall(PetscViewerASCIIPrintf(subviewer,"Col IS %" PetscInt_FMT "\n",gs)); 134 PetscCall(ISView(colis[ss],subviewer)); 135 PetscCall(PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)rowis[ss]),&subviewer)); 136 ++s; 137 } 138 } 139 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 140 } 141 PetscCall(PetscViewerFlush(viewer)); 142 PetscCall(ISSort(rowis[0])); 143 PetscCall(ISSort(colis[0])); 144 nsubdomains = 1; 145 PetscCall(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 PetscCall(PetscViewerASCIIPrintf(viewer,"Submatrices (repetition 1):\n")); 152 for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s; 153 PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)submats,&gnsubdomains,gsubdomainnums)); 154 PetscCall(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 PetscCall(PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer)); 162 PetscCall(MatView(submats[ss],subviewer)); 163 PetscCall(PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer)); 164 ++s; 165 } 166 } 167 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 168 } 169 PetscCall(PetscViewerFlush(viewer)); 170 if (rep == 1) goto cleanup; 171 nsubdomains = 1; 172 PetscCall(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 PetscCall(PetscViewerASCIIPrintf(viewer,"Submatrices (repetition 2):\n")); 179 for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s; 180 PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)submats,&gnsubdomains,gsubdomainnums)); 181 PetscCall(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 PetscCall(PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer)); 189 PetscCall(MatView(submats[ss],subviewer)); 190 PetscCall(PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer)); 191 ++s; 192 } 193 } 194 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 195 } 196 PetscCall(PetscViewerFlush(viewer)); 197 cleanup: 198 for (k=0;k<nsubdomains;++k) { 199 PetscCall(MatDestroy(submats+k)); 200 } 201 PetscCall(PetscFree(submats)); 202 for (k=0;k<nis;++k) { 203 PetscCall(ISDestroy(rowis+k)); 204 PetscCall(ISDestroy(colis+k)); 205 } 206 PetscCall(MatDestroy(&A)); 207 PetscCallMPI(MPI_Comm_free(&subcomm)); 208 PetscCall(PetscFinalize()); 209 return 0; 210 } 211 212 /*TEST 213 214 test: 215 nsize: 2 216 args: -total_subdomains 1 217 output_file: output/ex183_2_1.out 218 219 test: 220 suffix: 2 221 nsize: 3 222 args: -total_subdomains 2 223 output_file: output/ex183_3_2.out 224 225 test: 226 suffix: 3 227 nsize: 4 228 args: -total_subdomains 2 229 output_file: output/ex183_4_2.out 230 231 test: 232 suffix: 4 233 nsize: 6 234 args: -total_subdomains 2 235 output_file: output/ex183_6_2.out 236 237 TEST*/ 238