1 2 /* 3 Defines the basic matrix operations for the ADJ adjacency list matrix data-structure. 4 */ 5 #include <../src/mat/impls/adj/mpi/mpiadj.h> /*I "petscmat.h" I*/ 6 #include <petscsf.h> 7 8 /* 9 * The interface should be easy to use for both MatCreateSubMatrix (parallel sub-matrix) and MatCreateSubMatrices (sequential sub-matrices) 10 * */ 11 static PetscErrorCode MatCreateSubMatrix_MPIAdj_data(Mat adj,IS irows, IS icols, PetscInt **sadj_xadj,PetscInt **sadj_adjncy,PetscInt **sadj_values) 12 { 13 PetscInt nlrows_is,icols_n,i,j,nroots,nleaves,owner,rlocalindex,*ncols_send,*ncols_recv; 14 PetscInt nlrows_mat,*adjncy_recv,Ncols_recv,Ncols_send,*xadj_recv,*values_recv; 15 PetscInt *ncols_recv_offsets,loc,rnclos,*sadjncy,*sxadj,*svalues,isvalue; 16 const PetscInt *irows_indices,*icols_indices,*xadj, *adjncy; 17 Mat_MPIAdj *a = (Mat_MPIAdj*)adj->data; 18 PetscLayout rmap; 19 MPI_Comm comm; 20 PetscSF sf; 21 PetscSFNode *iremote; 22 PetscBool done; 23 PetscErrorCode ierr; 24 25 PetscFunctionBegin; 26 ierr = PetscObjectGetComm((PetscObject)adj,&comm);CHKERRQ(ierr); 27 ierr = MatGetLayouts(adj,&rmap,NULL);CHKERRQ(ierr); 28 ierr = ISGetLocalSize(irows,&nlrows_is);CHKERRQ(ierr); 29 ierr = ISGetIndices(irows,&irows_indices);CHKERRQ(ierr); 30 ierr = PetscCalloc1(nlrows_is,&iremote);CHKERRQ(ierr); 31 /* construct sf graph*/ 32 nleaves = nlrows_is; 33 for (i=0; i<nlrows_is; i++){ 34 owner = -1; 35 rlocalindex = -1; 36 ierr = PetscLayoutFindOwnerIndex(rmap,irows_indices[i],&owner,&rlocalindex);CHKERRQ(ierr); 37 iremote[i].rank = owner; 38 iremote[i].index = rlocalindex; 39 } 40 ierr = MatGetRowIJ(adj,0,PETSC_FALSE,PETSC_FALSE,&nlrows_mat,&xadj,&adjncy,&done);CHKERRQ(ierr); 41 ierr = PetscCalloc4(nlrows_mat,&ncols_send,nlrows_is,&xadj_recv,nlrows_is+1,&ncols_recv_offsets,nlrows_is,&ncols_recv);CHKERRQ(ierr); 42 nroots = nlrows_mat; 43 for (i=0; i<nlrows_mat; i++){ 44 ncols_send[i] = xadj[i+1]-xadj[i]; 45 } 46 ierr = PetscSFCreate(comm,&sf);CHKERRQ(ierr); 47 ierr = PetscSFSetGraph(sf,nroots,nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 48 ierr = PetscSFSetType(sf,PETSCSFBASIC);CHKERRQ(ierr); 49 ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr); 50 ierr = PetscSFBcastBegin(sf,MPIU_INT,ncols_send,ncols_recv);CHKERRQ(ierr); 51 ierr = PetscSFBcastEnd(sf,MPIU_INT,ncols_send,ncols_recv);CHKERRQ(ierr); 52 ierr = PetscSFBcastBegin(sf,MPIU_INT,xadj,xadj_recv);CHKERRQ(ierr); 53 ierr = PetscSFBcastEnd(sf,MPIU_INT,xadj,xadj_recv);CHKERRQ(ierr); 54 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 55 Ncols_recv =0; 56 for (i=0; i<nlrows_is; i++){ 57 Ncols_recv += ncols_recv[i]; 58 ncols_recv_offsets[i+1] = ncols_recv[i]+ncols_recv_offsets[i]; 59 } 60 Ncols_send = 0; 61 for(i=0; i<nlrows_mat; i++){ 62 Ncols_send += ncols_send[i]; 63 } 64 ierr = PetscCalloc1(Ncols_recv,&iremote);CHKERRQ(ierr); 65 ierr = PetscCalloc1(Ncols_recv,&adjncy_recv);CHKERRQ(ierr); 66 nleaves = Ncols_recv; 67 Ncols_recv = 0; 68 for (i=0; i<nlrows_is; i++){ 69 ierr = PetscLayoutFindOwner(rmap,irows_indices[i],&owner);CHKERRQ(ierr); 70 for (j=0; j<ncols_recv[i]; j++){ 71 iremote[Ncols_recv].rank = owner; 72 iremote[Ncols_recv++].index = xadj_recv[i]+j; 73 } 74 } 75 ierr = ISRestoreIndices(irows,&irows_indices);CHKERRQ(ierr); 76 /*if we need to deal with edge weights ???*/ 77 if (a->values){isvalue=1;} else {isvalue=0;} 78 /*involve a global communication */ 79 /*ierr = MPIU_Allreduce(&isvalue,&isvalue,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr);*/ 80 if (isvalue){ierr = PetscCalloc1(Ncols_recv,&values_recv);CHKERRQ(ierr);} 81 nroots = Ncols_send; 82 ierr = PetscSFCreate(comm,&sf);CHKERRQ(ierr); 83 ierr = PetscSFSetGraph(sf,nroots,nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 84 ierr = PetscSFSetType(sf,PETSCSFBASIC);CHKERRQ(ierr); 85 ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr); 86 ierr = PetscSFBcastBegin(sf,MPIU_INT,adjncy,adjncy_recv);CHKERRQ(ierr); 87 ierr = PetscSFBcastEnd(sf,MPIU_INT,adjncy,adjncy_recv);CHKERRQ(ierr); 88 if (isvalue){ 89 ierr = PetscSFBcastBegin(sf,MPIU_INT,a->values,values_recv);CHKERRQ(ierr); 90 ierr = PetscSFBcastEnd(sf,MPIU_INT,a->values,values_recv);CHKERRQ(ierr); 91 } 92 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 93 ierr = MatRestoreRowIJ(adj,0,PETSC_FALSE,PETSC_FALSE,&nlrows_mat,&xadj,&adjncy,&done);CHKERRQ(ierr); 94 ierr = ISGetLocalSize(icols,&icols_n);CHKERRQ(ierr); 95 ierr = ISGetIndices(icols,&icols_indices);CHKERRQ(ierr); 96 rnclos = 0; 97 for (i=0; i<nlrows_is; i++){ 98 for (j=ncols_recv_offsets[i]; j<ncols_recv_offsets[i+1]; j++){ 99 ierr = PetscFindInt(adjncy_recv[j], icols_n, icols_indices, &loc);CHKERRQ(ierr); 100 if (loc<0){ 101 adjncy_recv[j] = -1; 102 if (isvalue) values_recv[j] = -1; 103 ncols_recv[i]--; 104 } else { 105 rnclos++; 106 } 107 } 108 } 109 ierr = ISRestoreIndices(icols,&icols_indices);CHKERRQ(ierr); 110 ierr = PetscCalloc1(rnclos,&sadjncy);CHKERRQ(ierr); 111 if (isvalue) {ierr = PetscCalloc1(rnclos,&svalues);CHKERRQ(ierr);} 112 ierr = PetscCalloc1(nlrows_is+1,&sxadj);CHKERRQ(ierr); 113 rnclos = 0; 114 for(i=0; i<nlrows_is; i++){ 115 for(j=ncols_recv_offsets[i]; j<ncols_recv_offsets[i+1]; j++){ 116 if (adjncy_recv[j]<0) continue; 117 sadjncy[rnclos] = adjncy_recv[j]; 118 if (isvalue) svalues[rnclos] = values_recv[j]; 119 rnclos++; 120 } 121 } 122 for(i=0; i<nlrows_is; i++){ 123 sxadj[i+1] = sxadj[i]+ncols_recv[i]; 124 } 125 if (sadj_xadj) { *sadj_xadj = sxadj;} else { ierr = PetscFree(sxadj);CHKERRQ(ierr);} 126 if (sadj_adjncy){ *sadj_adjncy = sadjncy;} else { ierr = PetscFree(sadjncy);CHKERRQ(ierr);} 127 if (sadj_values){ 128 if (isvalue) *sadj_values = svalues; else *sadj_values=0; 129 } else { 130 if (isvalue) {ierr = PetscFree(svalues);CHKERRQ(ierr);} 131 } 132 ierr = PetscFree4(ncols_send,xadj_recv,ncols_recv_offsets,ncols_recv);CHKERRQ(ierr); 133 ierr = PetscFree(adjncy_recv);CHKERRQ(ierr); 134 if (isvalue) {ierr = PetscFree(values_recv);CHKERRQ(ierr);} 135 PetscFunctionReturn(0); 136 } 137 138 static PetscErrorCode MatCreateSubMatrices_MPIAdj_Private(Mat mat,PetscInt n,const IS irow[],const IS icol[],PetscBool subcomm,MatReuse scall,Mat *submat[]) 139 { 140 PetscInt i,irow_n,icol_n,*sxadj,*sadjncy,*svalues; 141 PetscInt *indices,nindx,j,k,loc; 142 PetscMPIInt issame; 143 const PetscInt *irow_indices,*icol_indices; 144 MPI_Comm scomm_row,scomm_col,scomm_mat; 145 PetscErrorCode ierr; 146 147 PetscFunctionBegin; 148 nindx = 0; 149 /* 150 * Estimate a maximum number for allocating memory 151 */ 152 for(i=0; i<n; i++){ 153 ierr = ISGetLocalSize(irow[i],&irow_n);CHKERRQ(ierr); 154 ierr = ISGetLocalSize(icol[i],&icol_n);CHKERRQ(ierr); 155 nindx = nindx>(irow_n+icol_n)? nindx:(irow_n+icol_n); 156 } 157 ierr = PetscCalloc1(nindx,&indices);CHKERRQ(ierr); 158 /* construct a submat */ 159 for(i=0; i<n; i++){ 160 if (subcomm){ 161 ierr = PetscObjectGetComm((PetscObject)irow[i],&scomm_row);CHKERRQ(ierr); 162 ierr = PetscObjectGetComm((PetscObject)icol[i],&scomm_col);CHKERRQ(ierr); 163 ierr = MPI_Comm_compare(scomm_row,scomm_col,&issame);CHKERRQ(ierr); 164 if (issame != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row index set must have the same comm as the col index set\n"); 165 ierr = MPI_Comm_compare(scomm_row,PETSC_COMM_SELF,&issame);CHKERRQ(ierr); 166 if (issame == MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP," can not use PETSC_COMM_SELF as comm when extracting a parallel submatrix\n"); 167 } else { 168 scomm_row = PETSC_COMM_SELF; 169 } 170 /*get sub-matrix data*/ 171 sxadj=0; sadjncy=0; svalues=0; 172 ierr = MatCreateSubMatrix_MPIAdj_data(mat,irow[i],icol[i],&sxadj,&sadjncy,&svalues);CHKERRQ(ierr); 173 ierr = ISGetLocalSize(irow[i],&irow_n);CHKERRQ(ierr); 174 ierr = ISGetLocalSize(icol[i],&icol_n);CHKERRQ(ierr); 175 ierr = ISGetIndices(irow[i],&irow_indices);CHKERRQ(ierr); 176 ierr = PetscMemcpy(indices,irow_indices,sizeof(PetscInt)*irow_n);CHKERRQ(ierr); 177 ierr = ISRestoreIndices(irow[i],&irow_indices);CHKERRQ(ierr); 178 ierr = ISGetIndices(icol[i],&icol_indices);CHKERRQ(ierr); 179 ierr = PetscMemcpy(indices+irow_n,icol_indices,sizeof(PetscInt)*icol_n);CHKERRQ(ierr); 180 ierr = ISRestoreIndices(icol[i],&icol_indices);CHKERRQ(ierr); 181 nindx = irow_n+icol_n; 182 ierr = PetscSortRemoveDupsInt(&nindx,indices);CHKERRQ(ierr); 183 /* renumber columns */ 184 for(j=0; j<irow_n; j++){ 185 for(k=sxadj[j]; k<sxadj[j+1]; k++){ 186 ierr = PetscFindInt(sadjncy[k],nindx,indices,&loc);CHKERRQ(ierr); 187 #if PETSC_USE_DEBUG 188 if (loc<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"can not find col %d \n",sadjncy[k]); 189 #endif 190 sadjncy[k] = loc; 191 } 192 } 193 if (scall==MAT_INITIAL_MATRIX){ 194 ierr = MatCreateMPIAdj(scomm_row,irow_n,icol_n,sxadj,sadjncy,svalues,submat[i]);CHKERRQ(ierr); 195 } else { 196 Mat sadj = *(submat[i]); 197 Mat_MPIAdj *sa = (Mat_MPIAdj*)((sadj)->data); 198 ierr = PetscObjectGetComm((PetscObject)sadj,&scomm_mat);CHKERRQ(ierr); 199 ierr = MPI_Comm_compare(scomm_row,scomm_mat,&issame);CHKERRQ(ierr); 200 if (issame != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"submatrix must have the same comm as the col index set\n"); 201 ierr = PetscMemcpy(sa->i,sxadj,sizeof(PetscInt)*(irow_n+1));CHKERRQ(ierr); 202 ierr = PetscMemcpy(sa->j,sadjncy,sizeof(PetscInt)*sxadj[irow_n]);CHKERRQ(ierr); 203 if (svalues){ierr = PetscMemcpy(sa->values,svalues,sizeof(PetscInt)*sxadj[irow_n]);CHKERRQ(ierr);} 204 ierr = PetscFree(sxadj);CHKERRQ(ierr); 205 ierr = PetscFree(sadjncy);CHKERRQ(ierr); 206 if (svalues) {ierr = PetscFree(svalues);CHKERRQ(ierr);} 207 } 208 } 209 ierr = PetscFree(indices);CHKERRQ(ierr); 210 PetscFunctionReturn(0); 211 } 212 213 static PetscErrorCode MatCreateSubMatricesMPI_MPIAdj(Mat mat,PetscInt n, const IS irow[],const IS icol[],MatReuse scall,Mat *submat[]) 214 { 215 PetscErrorCode ierr; 216 /*get sub-matrices across a sub communicator */ 217 PetscFunctionBegin; 218 ierr = MatCreateSubMatrices_MPIAdj_Private(mat,n,irow,icol,PETSC_TRUE,scall,submat);CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 } 221 222 static PetscErrorCode MatCreateSubMatrices_MPIAdj(Mat mat,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[]) 223 { 224 PetscErrorCode ierr; 225 226 PetscFunctionBegin; 227 /*get sub-matrices based on PETSC_COMM_SELF */ 228 ierr = MatCreateSubMatrices_MPIAdj_Private(mat,n,irow,icol,PETSC_FALSE,scall,submat);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 static PetscErrorCode MatView_MPIAdj_ASCII(Mat A,PetscViewer viewer) 233 { 234 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 235 PetscErrorCode ierr; 236 PetscInt i,j,m = A->rmap->n; 237 const char *name; 238 PetscViewerFormat format; 239 240 PetscFunctionBegin; 241 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 242 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 243 if (format == PETSC_VIEWER_ASCII_INFO) { 244 PetscFunctionReturn(0); 245 } else if (format == PETSC_VIEWER_ASCII_MATLAB) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MATLAB format not supported"); 246 else { 247 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 248 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 249 for (i=0; i<m; i++) { 250 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"row %D:",i+A->rmap->rstart);CHKERRQ(ierr); 251 for (j=a->i[i]; j<a->i[i+1]; j++) { 252 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %D ",a->j[j]);CHKERRQ(ierr); 253 } 254 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr); 255 } 256 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 257 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 259 } 260 PetscFunctionReturn(0); 261 } 262 263 static PetscErrorCode MatView_MPIAdj(Mat A,PetscViewer viewer) 264 { 265 PetscErrorCode ierr; 266 PetscBool iascii; 267 268 PetscFunctionBegin; 269 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 270 if (iascii) { 271 ierr = MatView_MPIAdj_ASCII(A,viewer);CHKERRQ(ierr); 272 } 273 PetscFunctionReturn(0); 274 } 275 276 static PetscErrorCode MatDestroy_MPIAdj(Mat mat) 277 { 278 Mat_MPIAdj *a = (Mat_MPIAdj*)mat->data; 279 PetscErrorCode ierr; 280 281 PetscFunctionBegin; 282 #if defined(PETSC_USE_LOG) 283 PetscLogObjectState((PetscObject)mat,"Rows=%D, Cols=%D, NZ=%D",mat->rmap->n,mat->cmap->n,a->nz); 284 #endif 285 ierr = PetscFree(a->diag);CHKERRQ(ierr); 286 if (a->freeaij) { 287 if (a->freeaijwithfree) { 288 if (a->i) free(a->i); 289 if (a->j) free(a->j); 290 } else { 291 ierr = PetscFree(a->i);CHKERRQ(ierr); 292 ierr = PetscFree(a->j);CHKERRQ(ierr); 293 ierr = PetscFree(a->values);CHKERRQ(ierr); 294 } 295 } 296 ierr = PetscFree(a->rowvalues);CHKERRQ(ierr); 297 ierr = PetscFree(mat->data);CHKERRQ(ierr); 298 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 299 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMPIAdjSetPreallocation_C",NULL);CHKERRQ(ierr); 300 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMPIAdjCreateNonemptySubcommMat_C",NULL);CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 304 static PetscErrorCode MatSetOption_MPIAdj(Mat A,MatOption op,PetscBool flg) 305 { 306 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 switch (op) { 311 case MAT_SYMMETRIC: 312 case MAT_STRUCTURALLY_SYMMETRIC: 313 case MAT_HERMITIAN: 314 a->symmetric = flg; 315 break; 316 case MAT_SYMMETRY_ETERNAL: 317 break; 318 default: 319 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 320 break; 321 } 322 PetscFunctionReturn(0); 323 } 324 325 static PetscErrorCode MatGetRow_MPIAdj(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 326 { 327 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 row -= A->rmap->rstart; 332 if (row < 0 || row >= A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row out of range"); 333 *nz = a->i[row+1] - a->i[row]; 334 if (v) { 335 PetscInt j; 336 if (a->rowvalues_alloc < *nz) { 337 ierr = PetscFree(a->rowvalues);CHKERRQ(ierr); 338 a->rowvalues_alloc = PetscMax(a->rowvalues_alloc*2, *nz); 339 ierr = PetscMalloc1(a->rowvalues_alloc,&a->rowvalues);CHKERRQ(ierr); 340 } 341 for (j=0; j<*nz; j++){ 342 a->rowvalues[j] = a->values ? a->values[a->i[row]+j]:1.0; 343 } 344 *v = (*nz) ? a->rowvalues : NULL; 345 } 346 if (idx) *idx = (*nz) ? a->j + a->i[row] : NULL; 347 PetscFunctionReturn(0); 348 } 349 350 static PetscErrorCode MatRestoreRow_MPIAdj(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 351 { 352 PetscFunctionBegin; 353 PetscFunctionReturn(0); 354 } 355 356 static PetscErrorCode MatEqual_MPIAdj(Mat A,Mat B,PetscBool * flg) 357 { 358 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data,*b = (Mat_MPIAdj*)B->data; 359 PetscErrorCode ierr; 360 PetscBool flag; 361 362 PetscFunctionBegin; 363 /* If the matrix dimensions are not equal,or no of nonzeros */ 364 if ((A->rmap->n != B->rmap->n) ||(a->nz != b->nz)) { 365 flag = PETSC_FALSE; 366 } 367 368 /* if the a->i are the same */ 369 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),&flag);CHKERRQ(ierr); 370 371 /* if a->j are the same */ 372 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),&flag);CHKERRQ(ierr); 373 374 ierr = MPIU_Allreduce(&flag,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 375 PetscFunctionReturn(0); 376 } 377 378 static PetscErrorCode MatGetRowIJ_MPIAdj(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *m,const PetscInt *inia[],const PetscInt *inja[],PetscBool *done) 379 { 380 PetscInt i; 381 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 382 PetscInt **ia = (PetscInt**)inia,**ja = (PetscInt**)inja; 383 384 PetscFunctionBegin; 385 *m = A->rmap->n; 386 *ia = a->i; 387 *ja = a->j; 388 *done = PETSC_TRUE; 389 if (oshift) { 390 for (i=0; i<(*ia)[*m]; i++) { 391 (*ja)[i]++; 392 } 393 for (i=0; i<=(*m); i++) (*ia)[i]++; 394 } 395 PetscFunctionReturn(0); 396 } 397 398 static PetscErrorCode MatRestoreRowIJ_MPIAdj(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *m,const PetscInt *inia[],const PetscInt *inja[],PetscBool *done) 399 { 400 PetscInt i; 401 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 402 PetscInt **ia = (PetscInt**)inia,**ja = (PetscInt**)inja; 403 404 PetscFunctionBegin; 405 if (ia && a->i != *ia) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"ia passed back is not one obtained with MatGetRowIJ()"); 406 if (ja && a->j != *ja) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"ja passed back is not one obtained with MatGetRowIJ()"); 407 if (oshift) { 408 if (!ia) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"If oshift then you must passed in inia[] argument"); 409 if (!ja) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"If oshift then you must passed in inja[] argument"); 410 for (i=0; i<=(*m); i++) (*ia)[i]--; 411 for (i=0; i<(*ia)[*m]; i++) { 412 (*ja)[i]--; 413 } 414 } 415 PetscFunctionReturn(0); 416 } 417 418 PetscErrorCode MatConvertFrom_MPIAdj(Mat A,MatType type,MatReuse reuse,Mat *newmat) 419 { 420 Mat B; 421 PetscErrorCode ierr; 422 PetscInt i,m,N,nzeros = 0,*ia,*ja,len,rstart,cnt,j,*a; 423 const PetscInt *rj; 424 const PetscScalar *ra; 425 MPI_Comm comm; 426 427 PetscFunctionBegin; 428 ierr = MatGetSize(A,NULL,&N);CHKERRQ(ierr); 429 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 430 ierr = MatGetOwnershipRange(A,&rstart,NULL);CHKERRQ(ierr); 431 432 /* count the number of nonzeros per row */ 433 for (i=0; i<m; i++) { 434 ierr = MatGetRow(A,i+rstart,&len,&rj,NULL);CHKERRQ(ierr); 435 for (j=0; j<len; j++) { 436 if (rj[j] == i+rstart) {len--; break;} /* don't count diagonal */ 437 } 438 nzeros += len; 439 ierr = MatRestoreRow(A,i+rstart,&len,&rj,NULL);CHKERRQ(ierr); 440 } 441 442 /* malloc space for nonzeros */ 443 ierr = PetscMalloc1(nzeros+1,&a);CHKERRQ(ierr); 444 ierr = PetscMalloc1(N+1,&ia);CHKERRQ(ierr); 445 ierr = PetscMalloc1(nzeros+1,&ja);CHKERRQ(ierr); 446 447 nzeros = 0; 448 ia[0] = 0; 449 for (i=0; i<m; i++) { 450 ierr = MatGetRow(A,i+rstart,&len,&rj,&ra);CHKERRQ(ierr); 451 cnt = 0; 452 for (j=0; j<len; j++) { 453 if (rj[j] != i+rstart) { /* if not diagonal */ 454 a[nzeros+cnt] = (PetscInt) PetscAbsScalar(ra[j]); 455 ja[nzeros+cnt++] = rj[j]; 456 } 457 } 458 ierr = MatRestoreRow(A,i+rstart,&len,&rj,&ra);CHKERRQ(ierr); 459 nzeros += cnt; 460 ia[i+1] = nzeros; 461 } 462 463 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 464 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 465 ierr = MatSetSizes(B,m,PETSC_DETERMINE,PETSC_DETERMINE,N);CHKERRQ(ierr); 466 ierr = MatSetType(B,type);CHKERRQ(ierr); 467 ierr = MatMPIAdjSetPreallocation(B,ia,ja,a);CHKERRQ(ierr); 468 469 if (reuse == MAT_INPLACE_MATRIX) { 470 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 471 } else { 472 *newmat = B; 473 } 474 PetscFunctionReturn(0); 475 } 476 477 /* -------------------------------------------------------------------*/ 478 static struct _MatOps MatOps_Values = {0, 479 MatGetRow_MPIAdj, 480 MatRestoreRow_MPIAdj, 481 0, 482 /* 4*/ 0, 483 0, 484 0, 485 0, 486 0, 487 0, 488 /*10*/ 0, 489 0, 490 0, 491 0, 492 0, 493 /*15*/ 0, 494 MatEqual_MPIAdj, 495 0, 496 0, 497 0, 498 /*20*/ 0, 499 0, 500 MatSetOption_MPIAdj, 501 0, 502 /*24*/ 0, 503 0, 504 0, 505 0, 506 0, 507 /*29*/ 0, 508 0, 509 0, 510 0, 511 0, 512 /*34*/ 0, 513 0, 514 0, 515 0, 516 0, 517 /*39*/ 0, 518 MatCreateSubMatrices_MPIAdj, 519 0, 520 0, 521 0, 522 /*44*/ 0, 523 0, 524 MatShift_Basic, 525 0, 526 0, 527 /*49*/ 0, 528 MatGetRowIJ_MPIAdj, 529 MatRestoreRowIJ_MPIAdj, 530 0, 531 0, 532 /*54*/ 0, 533 0, 534 0, 535 0, 536 0, 537 /*59*/ 0, 538 MatDestroy_MPIAdj, 539 MatView_MPIAdj, 540 MatConvertFrom_MPIAdj, 541 0, 542 /*64*/ 0, 543 0, 544 0, 545 0, 546 0, 547 /*69*/ 0, 548 0, 549 0, 550 0, 551 0, 552 /*74*/ 0, 553 0, 554 0, 555 0, 556 0, 557 /*79*/ 0, 558 0, 559 0, 560 0, 561 0, 562 /*84*/ 0, 563 0, 564 0, 565 0, 566 0, 567 /*89*/ 0, 568 0, 569 0, 570 0, 571 0, 572 /*94*/ 0, 573 0, 574 0, 575 0, 576 0, 577 /*99*/ 0, 578 0, 579 0, 580 0, 581 0, 582 /*104*/ 0, 583 0, 584 0, 585 0, 586 0, 587 /*109*/ 0, 588 0, 589 0, 590 0, 591 0, 592 /*114*/ 0, 593 0, 594 0, 595 0, 596 0, 597 /*119*/ 0, 598 0, 599 0, 600 0, 601 0, 602 /*124*/ 0, 603 0, 604 0, 605 0, 606 MatCreateSubMatricesMPI_MPIAdj, 607 /*129*/ 0, 608 0, 609 0, 610 0, 611 0, 612 /*134*/ 0, 613 0, 614 0, 615 0, 616 0, 617 /*139*/ 0, 618 0, 619 0 620 }; 621 622 static PetscErrorCode MatMPIAdjSetPreallocation_MPIAdj(Mat B,PetscInt *i,PetscInt *j,PetscInt *values) 623 { 624 Mat_MPIAdj *b = (Mat_MPIAdj*)B->data; 625 PetscErrorCode ierr; 626 #if defined(PETSC_USE_DEBUG) 627 PetscInt ii; 628 #endif 629 630 PetscFunctionBegin; 631 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 632 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 633 634 #if defined(PETSC_USE_DEBUG) 635 if (i[0] != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"First i[] index must be zero, instead it is %D\n",i[0]); 636 for (ii=1; ii<B->rmap->n; ii++) { 637 if (i[ii] < 0 || i[ii] < i[ii-1]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i[%D]=%D index is out of range: i[%D]=%D",ii,i[ii],ii-1,i[ii-1]); 638 } 639 for (ii=0; ii<i[B->rmap->n]; ii++) { 640 if (j[ii] < 0 || j[ii] >= B->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index %D out of range %D\n",ii,j[ii]); 641 } 642 #endif 643 B->preallocated = PETSC_TRUE; 644 645 b->j = j; 646 b->i = i; 647 b->values = values; 648 649 b->nz = i[B->rmap->n]; 650 b->diag = 0; 651 b->symmetric = PETSC_FALSE; 652 b->freeaij = PETSC_TRUE; 653 654 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 655 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658 659 static PetscErrorCode MatMPIAdjCreateNonemptySubcommMat_MPIAdj(Mat A,Mat *B) 660 { 661 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 662 PetscErrorCode ierr; 663 const PetscInt *ranges; 664 MPI_Comm acomm,bcomm; 665 MPI_Group agroup,bgroup; 666 PetscMPIInt i,rank,size,nranks,*ranks; 667 668 PetscFunctionBegin; 669 *B = NULL; 670 ierr = PetscObjectGetComm((PetscObject)A,&acomm);CHKERRQ(ierr); 671 ierr = MPI_Comm_size(acomm,&size);CHKERRQ(ierr); 672 ierr = MPI_Comm_size(acomm,&rank);CHKERRQ(ierr); 673 ierr = MatGetOwnershipRanges(A,&ranges);CHKERRQ(ierr); 674 for (i=0,nranks=0; i<size; i++) { 675 if (ranges[i+1] - ranges[i] > 0) nranks++; 676 } 677 if (nranks == size) { /* All ranks have a positive number of rows, so we do not need to create a subcomm; */ 678 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 679 *B = A; 680 PetscFunctionReturn(0); 681 } 682 683 ierr = PetscMalloc1(nranks,&ranks);CHKERRQ(ierr); 684 for (i=0,nranks=0; i<size; i++) { 685 if (ranges[i+1] - ranges[i] > 0) ranks[nranks++] = i; 686 } 687 ierr = MPI_Comm_group(acomm,&agroup);CHKERRQ(ierr); 688 ierr = MPI_Group_incl(agroup,nranks,ranks,&bgroup);CHKERRQ(ierr); 689 ierr = PetscFree(ranks);CHKERRQ(ierr); 690 ierr = MPI_Comm_create(acomm,bgroup,&bcomm);CHKERRQ(ierr); 691 ierr = MPI_Group_free(&agroup);CHKERRQ(ierr); 692 ierr = MPI_Group_free(&bgroup);CHKERRQ(ierr); 693 if (bcomm != MPI_COMM_NULL) { 694 PetscInt m,N; 695 Mat_MPIAdj *b; 696 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 697 ierr = MatGetSize(A,NULL,&N);CHKERRQ(ierr); 698 ierr = MatCreateMPIAdj(bcomm,m,N,a->i,a->j,a->values,B);CHKERRQ(ierr); 699 b = (Mat_MPIAdj*)(*B)->data; 700 b->freeaij = PETSC_FALSE; 701 ierr = MPI_Comm_free(&bcomm);CHKERRQ(ierr); 702 } 703 PetscFunctionReturn(0); 704 } 705 706 PetscErrorCode MatMPIAdjToSeq_MPIAdj(Mat A,Mat *B) 707 { 708 PetscErrorCode ierr; 709 PetscInt M,N,*II,*J,NZ,nz,m,nzstart,i; 710 PetscInt *Values = NULL; 711 Mat_MPIAdj *adj = (Mat_MPIAdj*)A->data; 712 PetscMPIInt mnz,mm,*allnz,*allm,size,*dispnz,*dispm; 713 714 PetscFunctionBegin; 715 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 716 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 717 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 718 nz = adj->nz; 719 ierr = MPI_Allreduce(&nz,&NZ,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 720 ierr = PetscMPIIntCast(nz,&mnz);CHKERRQ(ierr); 721 ierr = PetscMalloc2(size,&allnz,size,&dispnz);CHKERRQ(ierr); 722 ierr = MPI_Allgather(&mnz,1,MPI_INT,allnz,size,MPI_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 723 dispnz[0] = 0; for (i=1; i<size; i++) dispnz[i] = dispnz[i-1]+ allnz[i-1]; 724 if (adj->values) { 725 ierr = PetscMalloc1(NZ,&Values);CHKERRQ(ierr); 726 ierr = MPI_Allgatherv(adj->values,mnz,MPIU_INT,Values,allnz,dispnz,MPIU_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 727 } 728 ierr = PetscMalloc1(NZ,&J);CHKERRQ(ierr); 729 ierr = MPI_Allgatherv(adj->j,mnz,MPIU_INT,J,allnz,dispnz,MPIU_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 730 ierr = PetscFree2(allnz,dispnz);CHKERRQ(ierr); 731 ierr = MPI_Scan(&nz,&nzstart,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 732 nzstart -= nz; 733 /* shift the i[] values so they will be correct after being received */ 734 for (i=0; i<m; i++) adj->i[i] += nzstart; 735 ierr = PetscMalloc1(M+1,&II);CHKERRQ(ierr); 736 ierr = PetscMPIIntCast(m,&mm);CHKERRQ(ierr); 737 ierr = PetscMalloc2(size,&allm,size,&dispm);CHKERRQ(ierr); 738 ierr = MPI_Allgather(&mm,1,MPI_INT,allm,size,MPI_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 739 dispm[0] = 0; for (i=1; i<size; i++) dispm[i] = dispm[i-1]+ allm[i-1]; 740 ierr = MPI_Allgatherv(adj->i,mm,MPIU_INT,II,allm,dispm,MPIU_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 741 ierr = PetscFree2(allm,dispm);CHKERRQ(ierr); 742 II[M] = NZ; 743 /* shift the i[] values back */ 744 for (i=0; i<m; i++) adj->i[i] -= nzstart; 745 ierr = MatCreateMPIAdj(PETSC_COMM_SELF,M,N,II,J,Values,B);CHKERRQ(ierr); 746 PetscFunctionReturn(0); 747 } 748 749 /*@ 750 MatMPIAdjCreateNonemptySubcommMat - create the same MPIAdj matrix on a subcommunicator containing only processes owning a positive number of rows 751 752 Collective 753 754 Input Arguments: 755 . A - original MPIAdj matrix 756 757 Output Arguments: 758 . B - matrix on subcommunicator, NULL on ranks that owned zero rows of A 759 760 Level: developer 761 762 Note: 763 This function is mostly useful for internal use by mesh partitioning packages that require that every process owns at least one row. 764 765 The matrix B should be destroyed with MatDestroy(). The arrays are not copied, so B should be destroyed before A is destroyed. 766 767 .seealso: MatCreateMPIAdj() 768 @*/ 769 PetscErrorCode MatMPIAdjCreateNonemptySubcommMat(Mat A,Mat *B) 770 { 771 PetscErrorCode ierr; 772 773 PetscFunctionBegin; 774 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 775 ierr = PetscUseMethod(A,"MatMPIAdjCreateNonemptySubcommMat_C",(Mat,Mat*),(A,B));CHKERRQ(ierr); 776 PetscFunctionReturn(0); 777 } 778 779 /*MC 780 MATMPIADJ - MATMPIADJ = "mpiadj" - A matrix type to be used for distributed adjacency matrices, 781 intended for use constructing orderings and partitionings. 782 783 Level: beginner 784 785 .seealso: MatCreateMPIAdj 786 M*/ 787 788 PETSC_EXTERN PetscErrorCode MatCreate_MPIAdj(Mat B) 789 { 790 Mat_MPIAdj *b; 791 PetscErrorCode ierr; 792 793 PetscFunctionBegin; 794 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 795 B->data = (void*)b; 796 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 797 B->assembled = PETSC_FALSE; 798 799 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPIAdjSetPreallocation_C",MatMPIAdjSetPreallocation_MPIAdj);CHKERRQ(ierr); 800 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPIAdjCreateNonemptySubcommMat_C",MatMPIAdjCreateNonemptySubcommMat_MPIAdj);CHKERRQ(ierr); 801 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPIAdjToSeq_C",MatMPIAdjToSeq_MPIAdj);CHKERRQ(ierr); 802 ierr = PetscObjectChangeTypeName((PetscObject)B,MATMPIADJ);CHKERRQ(ierr); 803 PetscFunctionReturn(0); 804 } 805 806 /*@C 807 MatMPIAdjToSeq- Converts an parallel MPIAdj matrix to complete MPIAdj on each process (needed by sequential preconditioners) 808 809 Logically Collective on MPI_Comm 810 811 Input Parameter: 812 . A - the matrix 813 814 Output Parameter: 815 . B - the same matrix on all processes 816 817 Level: intermediate 818 819 .seealso: MatCreate(), MatCreateMPIAdj(), MatSetValues() 820 @*/ 821 PetscErrorCode MatMPIAdjToSeq(Mat A,Mat *B) 822 { 823 PetscErrorCode ierr; 824 825 PetscFunctionBegin; 826 ierr = PetscUseMethod(A,"MatMPIAdjToSeq_C",(Mat,Mat*),(A,B));CHKERRQ(ierr); 827 PetscFunctionReturn(0); 828 } 829 830 /*@C 831 MatMPIAdjSetPreallocation - Sets the array used for storing the matrix elements 832 833 Logically Collective on MPI_Comm 834 835 Input Parameters: 836 + A - the matrix 837 . i - the indices into j for the start of each row 838 . j - the column indices for each row (sorted for each row). 839 The indices in i and j start with zero (NOT with one). 840 - values - [optional] edge weights 841 842 Level: intermediate 843 844 .seealso: MatCreate(), MatCreateMPIAdj(), MatSetValues() 845 @*/ 846 PetscErrorCode MatMPIAdjSetPreallocation(Mat B,PetscInt *i,PetscInt *j,PetscInt *values) 847 { 848 PetscErrorCode ierr; 849 850 PetscFunctionBegin; 851 ierr = PetscTryMethod(B,"MatMPIAdjSetPreallocation_C",(Mat,PetscInt*,PetscInt*,PetscInt*),(B,i,j,values));CHKERRQ(ierr); 852 PetscFunctionReturn(0); 853 } 854 855 /*@C 856 MatCreateMPIAdj - Creates a sparse matrix representing an adjacency list. 857 The matrix does not have numerical values associated with it, but is 858 intended for ordering (to reduce bandwidth etc) and partitioning. 859 860 Collective on MPI_Comm 861 862 Input Parameters: 863 + comm - MPI communicator 864 . m - number of local rows 865 . N - number of global columns 866 . i - the indices into j for the start of each row 867 . j - the column indices for each row (sorted for each row). 868 The indices in i and j start with zero (NOT with one). 869 - values -[optional] edge weights 870 871 Output Parameter: 872 . A - the matrix 873 874 Level: intermediate 875 876 Notes: This matrix object does not support most matrix operations, include 877 MatSetValues(). 878 You must NOT free the ii, values and jj arrays yourself. PETSc will free them 879 when the matrix is destroyed; you must allocate them with PetscMalloc(). If you 880 call from Fortran you need not create the arrays with PetscMalloc(). 881 Should not include the matrix diagonals. 882 883 If you already have a matrix, you can create its adjacency matrix by a call 884 to MatConvert, specifying a type of MATMPIADJ. 885 886 Possible values for MatSetOption() - MAT_STRUCTURALLY_SYMMETRIC 887 888 .seealso: MatCreate(), MatConvert(), MatGetOrdering() 889 @*/ 890 PetscErrorCode MatCreateMPIAdj(MPI_Comm comm,PetscInt m,PetscInt N,PetscInt *i,PetscInt *j,PetscInt *values,Mat *A) 891 { 892 PetscErrorCode ierr; 893 894 PetscFunctionBegin; 895 ierr = MatCreate(comm,A);CHKERRQ(ierr); 896 ierr = MatSetSizes(*A,m,PETSC_DETERMINE,PETSC_DETERMINE,N);CHKERRQ(ierr); 897 ierr = MatSetType(*A,MATMPIADJ);CHKERRQ(ierr); 898 ierr = MatMPIAdjSetPreallocation(*A,i,j,values);CHKERRQ(ierr); 899 PetscFunctionReturn(0); 900 } 901 902 903 904 905 906 907 908