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