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 /*@ 707 MatMPIAdjCreateNonemptySubcommMat - create the same MPIAdj matrix on a subcommunicator containing only processes owning a positive number of rows 708 709 Collective 710 711 Input Arguments: 712 . A - original MPIAdj matrix 713 714 Output Arguments: 715 . B - matrix on subcommunicator, NULL on ranks that owned zero rows of A 716 717 Level: developer 718 719 Note: 720 This function is mostly useful for internal use by mesh partitioning packages that require that every process owns at least one row. 721 722 The matrix B should be destroyed with MatDestroy(). The arrays are not copied, so B should be destroyed before A is destroyed. 723 724 .seealso: MatCreateMPIAdj() 725 @*/ 726 PetscErrorCode MatMPIAdjCreateNonemptySubcommMat(Mat A,Mat *B) 727 { 728 PetscErrorCode ierr; 729 730 PetscFunctionBegin; 731 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 732 ierr = PetscUseMethod(A,"MatMPIAdjCreateNonemptySubcommMat_C",(Mat,Mat*),(A,B));CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 /*MC 737 MATMPIADJ - MATMPIADJ = "mpiadj" - A matrix type to be used for distributed adjacency matrices, 738 intended for use constructing orderings and partitionings. 739 740 Level: beginner 741 742 .seealso: MatCreateMPIAdj 743 M*/ 744 745 PETSC_EXTERN PetscErrorCode MatCreate_MPIAdj(Mat B) 746 { 747 Mat_MPIAdj *b; 748 PetscErrorCode ierr; 749 750 PetscFunctionBegin; 751 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 752 B->data = (void*)b; 753 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 754 B->assembled = PETSC_FALSE; 755 756 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPIAdjSetPreallocation_C",MatMPIAdjSetPreallocation_MPIAdj);CHKERRQ(ierr); 757 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPIAdjCreateNonemptySubcommMat_C",MatMPIAdjCreateNonemptySubcommMat_MPIAdj);CHKERRQ(ierr); 758 ierr = PetscObjectChangeTypeName((PetscObject)B,MATMPIADJ);CHKERRQ(ierr); 759 PetscFunctionReturn(0); 760 } 761 762 /*@C 763 MatMPIAdjSetPreallocation - Sets the array used for storing the matrix elements 764 765 Logically Collective on MPI_Comm 766 767 Input Parameters: 768 + A - the matrix 769 . i - the indices into j for the start of each row 770 . j - the column indices for each row (sorted for each row). 771 The indices in i and j start with zero (NOT with one). 772 - values - [optional] edge weights 773 774 Level: intermediate 775 776 .seealso: MatCreate(), MatCreateMPIAdj(), MatSetValues() 777 @*/ 778 PetscErrorCode MatMPIAdjSetPreallocation(Mat B,PetscInt *i,PetscInt *j,PetscInt *values) 779 { 780 PetscErrorCode ierr; 781 782 PetscFunctionBegin; 783 ierr = PetscTryMethod(B,"MatMPIAdjSetPreallocation_C",(Mat,PetscInt*,PetscInt*,PetscInt*),(B,i,j,values));CHKERRQ(ierr); 784 PetscFunctionReturn(0); 785 } 786 787 /*@C 788 MatCreateMPIAdj - Creates a sparse matrix representing an adjacency list. 789 The matrix does not have numerical values associated with it, but is 790 intended for ordering (to reduce bandwidth etc) and partitioning. 791 792 Collective on MPI_Comm 793 794 Input Parameters: 795 + comm - MPI communicator 796 . m - number of local rows 797 . N - number of global columns 798 . i - the indices into j for the start of each row 799 . j - the column indices for each row (sorted for each row). 800 The indices in i and j start with zero (NOT with one). 801 - values -[optional] edge weights 802 803 Output Parameter: 804 . A - the matrix 805 806 Level: intermediate 807 808 Notes: This matrix object does not support most matrix operations, include 809 MatSetValues(). 810 You must NOT free the ii, values and jj arrays yourself. PETSc will free them 811 when the matrix is destroyed; you must allocate them with PetscMalloc(). If you 812 call from Fortran you need not create the arrays with PetscMalloc(). 813 Should not include the matrix diagonals. 814 815 If you already have a matrix, you can create its adjacency matrix by a call 816 to MatConvert, specifying a type of MATMPIADJ. 817 818 Possible values for MatSetOption() - MAT_STRUCTURALLY_SYMMETRIC 819 820 .seealso: MatCreate(), MatConvert(), MatGetOrdering() 821 @*/ 822 PetscErrorCode MatCreateMPIAdj(MPI_Comm comm,PetscInt m,PetscInt N,PetscInt *i,PetscInt *j,PetscInt *values,Mat *A) 823 { 824 PetscErrorCode ierr; 825 826 PetscFunctionBegin; 827 ierr = MatCreate(comm,A);CHKERRQ(ierr); 828 ierr = MatSetSizes(*A,m,PETSC_DETERMINE,PETSC_DETERMINE,N);CHKERRQ(ierr); 829 ierr = MatSetType(*A,MATMPIADJ);CHKERRQ(ierr); 830 ierr = MatMPIAdjSetPreallocation(*A,i,j,values);CHKERRQ(ierr); 831 PetscFunctionReturn(0); 832 } 833 834 835 836 837 838 839 840