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 #include <petscmat.h> 6 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; 20 21 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 23 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 24 25 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex183","Mat");PetscCall(ierr); 26 m = 5; 27 PetscCall(PetscOptionsInt("-m","Local matrix size","MatSetSizes",m,&m,&flg)); 28 total_subdomains = size-1; 29 PetscCall(PetscOptionsInt("-total_subdomains","Number of submatrices where 0 < n < comm size","MatCreateSubMatricesMPI",total_subdomains,&total_subdomains,&flg)); 30 permute_indices = PETSC_FALSE; 31 PetscCall(PetscOptionsBool("-permute_indices","Whether to permute indices before breaking them into subdomains","ISCreateGeneral",permute_indices,&permute_indices,&flg)); 32 hash = 7; 33 PetscCall(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 PetscCall(PetscOptionsInt("-rep","Number of times to carry out submatrix extractions; currently only 1 & 2 are supported",NULL,rep,&rep,&flg)); 36 ierr = PetscOptionsEnd();PetscCall(ierr); 37 38 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); 39 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); 40 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); 41 42 viewer = PETSC_VIEWER_STDOUT_WORLD; 43 /* Create logically sparse, but effectively dense matrix for easy verification of submatrix extraction correctness. */ 44 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 45 PetscCall(MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE)); 46 PetscCall(MatSetFromOptions(A)); 47 PetscCall(MatSetUp(A)); 48 PetscCall(MatGetSize(A,NULL,&N)); 49 PetscCall(MatGetLocalSize(A,NULL,&n)); 50 PetscCall(MatGetBlockSize(A,&bs)); 51 PetscCall(MatSeqAIJSetPreallocation(A,n,NULL)); 52 PetscCall(MatMPIAIJSetPreallocation(A,n,NULL,N-n,NULL)); 53 PetscCall(MatSeqBAIJSetPreallocation(A,bs,n/bs,NULL)); 54 PetscCall(MatMPIBAIJSetPreallocation(A,bs,n/bs,NULL,(N-n)/bs,NULL)); 55 PetscCall(MatSeqSBAIJSetPreallocation(A,bs,n/bs,NULL)); 56 PetscCall(MatMPISBAIJSetPreallocation(A,bs,n/bs,NULL,(N-n)/bs,NULL)); 57 58 PetscCall(PetscMalloc2(N,&cols,N,&vals)); 59 PetscCall(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 PetscCall(MatSetValues(A,1,&i,N,cols,vals,INSERT_VALUES)); 66 } 67 PetscCall(PetscFree2(cols,vals)); 68 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 69 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 70 71 PetscCall(PetscViewerASCIIPrintf(viewer,"Initial matrix:\n")); 72 PetscCall(MatView(A,viewer)); 73 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 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD,color,rank,&subcomm)); 82 PetscCallMPI(MPI_Comm_size(subcomm,&subsize)); 83 PetscCallMPI(MPI_Comm_rank(subcomm,&subrank)); 84 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 85 nis = 1; 86 PetscCall(PetscMalloc2(rend-rstart,&rowindices,rend-rstart,&colindices)); 87 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 PetscCall(ISCreateGeneral(subcomm,rend-rstart,rowindices,PETSC_COPY_VALUES,&rowis[0])); 98 PetscCall(ISCreateGeneral(subcomm,rend-rstart,colindices,PETSC_COPY_VALUES,&colis[0])); 99 PetscCall(ISSort(rowis[0])); 100 PetscCall(ISSort(colis[0])); 101 PetscCall(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 PetscCall(PetscViewerASCIIPrintf(viewer,"Subdomains")); 108 if (permute_indices) { 109 PetscCall(PetscViewerASCIIPrintf(viewer," (hash=%" PetscInt_FMT ")",hash)); 110 } 111 PetscCall(PetscViewerASCIIPrintf(viewer,":\n")); 112 PetscCall(PetscViewerFlush(viewer)); 113 114 nsubdomains = 1; 115 for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s; 116 PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)rowis,&gnsubdomains,gsubdomainnums)); 117 PetscCall(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 PetscCall(PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)rowis[ss]),&subviewer)); 125 PetscCall(PetscViewerASCIIPrintf(subviewer,"Row IS %" PetscInt_FMT "\n",gs)); 126 PetscCall(ISView(rowis[ss],subviewer)); 127 PetscCall(PetscViewerFlush(subviewer)); 128 PetscCall(PetscViewerASCIIPrintf(subviewer,"Col IS %" PetscInt_FMT "\n",gs)); 129 PetscCall(ISView(colis[ss],subviewer)); 130 PetscCall(PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)rowis[ss]),&subviewer)); 131 ++s; 132 } 133 } 134 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 135 } 136 PetscCall(PetscViewerFlush(viewer)); 137 PetscCall(ISSort(rowis[0])); 138 PetscCall(ISSort(colis[0])); 139 nsubdomains = 1; 140 PetscCall(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 PetscCall(PetscViewerASCIIPrintf(viewer,"Submatrices (repetition 1):\n")); 147 for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s; 148 PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)submats,&gnsubdomains,gsubdomainnums)); 149 PetscCall(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 PetscCall(PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer)); 157 PetscCall(MatView(submats[ss],subviewer)); 158 PetscCall(PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer)); 159 ++s; 160 } 161 } 162 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 163 } 164 PetscCall(PetscViewerFlush(viewer)); 165 if (rep == 1) goto cleanup; 166 nsubdomains = 1; 167 PetscCall(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 PetscCall(PetscViewerASCIIPrintf(viewer,"Submatrices (repetition 2):\n")); 174 for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s; 175 PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD,1,(PetscObject*)submats,&gnsubdomains,gsubdomainnums)); 176 PetscCall(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 PetscCall(PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer)); 184 PetscCall(MatView(submats[ss],subviewer)); 185 PetscCall(PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)submats[ss]),&subviewer)); 186 ++s; 187 } 188 } 189 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 190 } 191 PetscCall(PetscViewerFlush(viewer)); 192 cleanup: 193 for (k=0;k<nsubdomains;++k) { 194 PetscCall(MatDestroy(submats+k)); 195 } 196 PetscCall(PetscFree(submats)); 197 for (k=0;k<nis;++k) { 198 PetscCall(ISDestroy(rowis+k)); 199 PetscCall(ISDestroy(colis+k)); 200 } 201 PetscCall(MatDestroy(&A)); 202 PetscCallMPI(MPI_Comm_free(&subcomm)); 203 PetscCall(PetscFinalize()); 204 return 0; 205 } 206 207 /*TEST 208 209 test: 210 nsize: 2 211 args: -total_subdomains 1 212 output_file: output/ex183_2_1.out 213 214 test: 215 suffix: 2 216 nsize: 3 217 args: -total_subdomains 2 218 output_file: output/ex183_3_2.out 219 220 test: 221 suffix: 3 222 nsize: 4 223 args: -total_subdomains 2 224 output_file: output/ex183_4_2.out 225 226 test: 227 suffix: 4 228 nsize: 6 229 args: -total_subdomains 2 230 output_file: output/ex183_6_2.out 231 232 TEST*/ 233