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