Actual source code: mpiadj.c
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, 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;
25: PetscObjectGetComm((PetscObject)adj, &comm);
26: MatGetLayouts(adj, &rmap, NULL);
27: ISGetLocalSize(irows, &nlrows_is);
28: ISGetIndices(irows, &irows_indices);
29: PetscMalloc1(nlrows_is, &iremote);
30: /* construct sf graph*/
31: nleaves = nlrows_is;
32: for (i = 0; i < nlrows_is; i++) {
33: owner = -1;
34: rlocalindex = -1;
35: PetscLayoutFindOwnerIndex(rmap, irows_indices[i], &owner, &rlocalindex);
36: iremote[i].rank = owner;
37: iremote[i].index = rlocalindex;
38: }
39: MatGetRowIJ(adj, 0, PETSC_FALSE, PETSC_FALSE, &nlrows_mat, &xadj, &adjncy, &done);
40: PetscCalloc4(nlrows_mat, &ncols_send, nlrows_is, &xadj_recv, nlrows_is + 1, &ncols_recv_offsets, nlrows_is, &ncols_recv);
41: nroots = nlrows_mat;
42: for (i = 0; i < nlrows_mat; i++) ncols_send[i] = xadj[i + 1] - xadj[i];
43: PetscSFCreate(comm, &sf);
44: PetscSFSetGraph(sf, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
45: PetscSFSetType(sf, PETSCSFBASIC);
46: PetscSFSetFromOptions(sf);
47: PetscSFBcastBegin(sf, MPIU_INT, ncols_send, ncols_recv, MPI_REPLACE);
48: PetscSFBcastEnd(sf, MPIU_INT, ncols_send, ncols_recv, MPI_REPLACE);
49: PetscSFBcastBegin(sf, MPIU_INT, xadj, xadj_recv, MPI_REPLACE);
50: PetscSFBcastEnd(sf, MPIU_INT, xadj, xadj_recv, MPI_REPLACE);
51: PetscSFDestroy(&sf);
52: Ncols_recv = 0;
53: for (i = 0; i < nlrows_is; i++) {
54: Ncols_recv += ncols_recv[i];
55: ncols_recv_offsets[i + 1] = ncols_recv[i] + ncols_recv_offsets[i];
56: }
57: Ncols_send = 0;
58: for (i = 0; i < nlrows_mat; i++) Ncols_send += ncols_send[i];
59: PetscCalloc1(Ncols_recv, &iremote);
60: PetscCalloc1(Ncols_recv, &adjncy_recv);
61: nleaves = Ncols_recv;
62: Ncols_recv = 0;
63: for (i = 0; i < nlrows_is; i++) {
64: PetscLayoutFindOwner(rmap, irows_indices[i], &owner);
65: for (j = 0; j < ncols_recv[i]; j++) {
66: iremote[Ncols_recv].rank = owner;
67: iremote[Ncols_recv++].index = xadj_recv[i] + j;
68: }
69: }
70: ISRestoreIndices(irows, &irows_indices);
71: /*if we need to deal with edge weights ???*/
72: if (a->useedgeweights) PetscCalloc1(Ncols_recv, &values_recv);
73: nroots = Ncols_send;
74: PetscSFCreate(comm, &sf);
75: PetscSFSetGraph(sf, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
76: PetscSFSetType(sf, PETSCSFBASIC);
77: PetscSFSetFromOptions(sf);
78: PetscSFBcastBegin(sf, MPIU_INT, adjncy, adjncy_recv, MPI_REPLACE);
79: PetscSFBcastEnd(sf, MPIU_INT, adjncy, adjncy_recv, MPI_REPLACE);
80: if (a->useedgeweights) {
81: PetscSFBcastBegin(sf, MPIU_INT, a->values, values_recv, MPI_REPLACE);
82: PetscSFBcastEnd(sf, MPIU_INT, a->values, values_recv, MPI_REPLACE);
83: }
84: PetscSFDestroy(&sf);
85: MatRestoreRowIJ(adj, 0, PETSC_FALSE, PETSC_FALSE, &nlrows_mat, &xadj, &adjncy, &done);
86: ISGetLocalSize(icols, &icols_n);
87: ISGetIndices(icols, &icols_indices);
88: rnclos = 0;
89: for (i = 0; i < nlrows_is; i++) {
90: for (j = ncols_recv_offsets[i]; j < ncols_recv_offsets[i + 1]; j++) {
91: PetscFindInt(adjncy_recv[j], icols_n, icols_indices, &loc);
92: if (loc < 0) {
93: adjncy_recv[j] = -1;
94: if (a->useedgeweights) values_recv[j] = -1;
95: ncols_recv[i]--;
96: } else {
97: rnclos++;
98: }
99: }
100: }
101: ISRestoreIndices(icols, &icols_indices);
102: PetscCalloc1(rnclos, &sadjncy);
103: if (a->useedgeweights) PetscCalloc1(rnclos, &svalues);
104: PetscCalloc1(nlrows_is + 1, &sxadj);
105: rnclos = 0;
106: for (i = 0; i < nlrows_is; i++) {
107: for (j = ncols_recv_offsets[i]; j < ncols_recv_offsets[i + 1]; j++) {
108: if (adjncy_recv[j] < 0) continue;
109: sadjncy[rnclos] = adjncy_recv[j];
110: if (a->useedgeweights) svalues[rnclos] = values_recv[j];
111: rnclos++;
112: }
113: }
114: for (i = 0; i < nlrows_is; i++) sxadj[i + 1] = sxadj[i] + ncols_recv[i];
115: if (sadj_xadj) {
116: *sadj_xadj = sxadj;
117: } else PetscFree(sxadj);
118: if (sadj_adjncy) {
119: *sadj_adjncy = sadjncy;
120: } else PetscFree(sadjncy);
121: if (sadj_values) {
122: if (a->useedgeweights) *sadj_values = svalues;
123: else *sadj_values = NULL;
124: } else {
125: if (a->useedgeweights) PetscFree(svalues);
126: }
127: PetscFree4(ncols_send, xadj_recv, ncols_recv_offsets, ncols_recv);
128: PetscFree(adjncy_recv);
129: if (a->useedgeweights) PetscFree(values_recv);
130: return 0;
131: }
133: static PetscErrorCode MatCreateSubMatrices_MPIAdj_Private(Mat mat, PetscInt n, const IS irow[], const IS icol[], PetscBool subcomm, MatReuse scall, Mat *submat[])
134: {
135: PetscInt i, irow_n, icol_n, *sxadj, *sadjncy, *svalues;
136: PetscInt *indices, nindx, j, k, loc;
137: PetscMPIInt issame;
138: const PetscInt *irow_indices, *icol_indices;
139: MPI_Comm scomm_row, scomm_col, scomm_mat;
141: nindx = 0;
142: /*
143: * Estimate a maximum number for allocating memory
144: */
145: for (i = 0; i < n; i++) {
146: ISGetLocalSize(irow[i], &irow_n);
147: ISGetLocalSize(icol[i], &icol_n);
148: nindx = nindx > (irow_n + icol_n) ? nindx : (irow_n + icol_n);
149: }
150: PetscMalloc1(nindx, &indices);
151: /* construct a submat */
152: // if (scall == MAT_INITIAL_MATRIX) PetscMalloc1(n,submat);
154: for (i = 0; i < n; i++) {
155: if (subcomm) {
156: PetscObjectGetComm((PetscObject)irow[i], &scomm_row);
157: PetscObjectGetComm((PetscObject)icol[i], &scomm_col);
158: MPI_Comm_compare(scomm_row, scomm_col, &issame);
160: MPI_Comm_compare(scomm_row, PETSC_COMM_SELF, &issame);
162: } else {
163: scomm_row = PETSC_COMM_SELF;
164: }
165: /*get sub-matrix data*/
166: sxadj = NULL;
167: sadjncy = NULL;
168: svalues = NULL;
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: PetscArraycpy(indices, irow_indices, irow_n);
174: ISRestoreIndices(irow[i], &irow_indices);
175: ISGetIndices(icol[i], &icol_indices);
176: PetscArraycpy(indices + irow_n, icol_indices, 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);
185: sadjncy[k] = loc;
186: }
187: }
188: if (scall == MAT_INITIAL_MATRIX) {
189: MatCreateMPIAdj(scomm_row, irow_n, icol_n, sxadj, sadjncy, svalues, submat[i]);
190: } else {
191: Mat sadj = *(submat[i]);
192: Mat_MPIAdj *sa = (Mat_MPIAdj *)((sadj)->data);
193: PetscObjectGetComm((PetscObject)sadj, &scomm_mat);
194: MPI_Comm_compare(scomm_row, scomm_mat, &issame);
196: PetscArraycpy(sa->i, sxadj, irow_n + 1);
197: PetscArraycpy(sa->j, sadjncy, sxadj[irow_n]);
198: if (svalues) PetscArraycpy(sa->values, svalues, sxadj[irow_n]);
199: PetscFree(sxadj);
200: PetscFree(sadjncy);
201: PetscFree(svalues);
202: }
203: }
204: PetscFree(indices);
205: return 0;
206: }
208: static PetscErrorCode MatCreateSubMatricesMPI_MPIAdj(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
209: {
210: /*get sub-matrices across a sub communicator */
211: MatCreateSubMatrices_MPIAdj_Private(mat, n, irow, icol, PETSC_TRUE, scall, submat);
212: return 0;
213: }
215: static PetscErrorCode MatCreateSubMatrices_MPIAdj(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
216: {
217: /*get sub-matrices based on PETSC_COMM_SELF */
218: MatCreateSubMatrices_MPIAdj_Private(mat, n, irow, icol, PETSC_FALSE, scall, submat);
219: return 0;
220: }
222: static PetscErrorCode MatView_MPIAdj_ASCII(Mat A, PetscViewer viewer)
223: {
224: Mat_MPIAdj *a = (Mat_MPIAdj *)A->data;
225: PetscInt i, j, m = A->rmap->n;
226: const char *name;
227: PetscViewerFormat format;
229: PetscObjectGetName((PetscObject)A, &name);
230: PetscViewerGetFormat(viewer, &format);
231: if (format == PETSC_VIEWER_ASCII_INFO) {
232: return 0;
233: } else {
235: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
236: PetscViewerASCIIPushSynchronized(viewer);
237: for (i = 0; i < m; i++) {
238: PetscViewerASCIISynchronizedPrintf(viewer, "row %" PetscInt_FMT ":", i + A->rmap->rstart);
239: for (j = a->i[i]; j < a->i[i + 1]; j++) {
240: if (a->values) {
241: PetscViewerASCIISynchronizedPrintf(viewer, " (%" PetscInt_FMT ", %" PetscInt_FMT ") ", a->j[j], a->values[j]);
242: } else {
243: PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT " ", a->j[j]);
244: }
245: }
246: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
247: }
248: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
249: PetscViewerFlush(viewer);
250: PetscViewerASCIIPopSynchronized(viewer);
251: }
252: return 0;
253: }
255: static PetscErrorCode MatView_MPIAdj(Mat A, PetscViewer viewer)
256: {
257: PetscBool iascii;
259: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
260: if (iascii) MatView_MPIAdj_ASCII(A, viewer);
261: return 0;
262: }
264: static PetscErrorCode MatDestroy_MPIAdj(Mat mat)
265: {
266: Mat_MPIAdj *a = (Mat_MPIAdj *)mat->data;
268: #if defined(PETSC_USE_LOG)
269: PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, mat->rmap->n, mat->cmap->n, a->nz);
270: #endif
271: PetscFree(a->diag);
272: if (a->freeaij) {
273: if (a->freeaijwithfree) {
274: if (a->i) free(a->i);
275: if (a->j) free(a->j);
276: } else {
277: PetscFree(a->i);
278: PetscFree(a->j);
279: PetscFree(a->values);
280: }
281: }
282: PetscFree(a->rowvalues);
283: PetscFree(mat->data);
284: PetscObjectChangeTypeName((PetscObject)mat, NULL);
285: PetscObjectComposeFunction((PetscObject)mat, "MatMPIAdjSetPreallocation_C", NULL);
286: PetscObjectComposeFunction((PetscObject)mat, "MatMPIAdjCreateNonemptySubcommMat_C", NULL);
287: PetscObjectComposeFunction((PetscObject)mat, "MatMPIAdjToSeq_C", NULL);
288: PetscObjectComposeFunction((PetscObject)mat, "MatMPIAdjToSeqRankZero_C", NULL);
289: return 0;
290: }
292: static PetscErrorCode MatSetOption_MPIAdj(Mat A, MatOption op, PetscBool flg)
293: {
294: Mat_MPIAdj *a = (Mat_MPIAdj *)A->data;
296: switch (op) {
297: case MAT_SYMMETRIC:
298: case MAT_STRUCTURALLY_SYMMETRIC:
299: case MAT_HERMITIAN:
300: case MAT_SPD:
301: a->symmetric = flg;
302: break;
303: case MAT_SYMMETRY_ETERNAL:
304: case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
305: case MAT_SPD_ETERNAL:
306: break;
307: default:
308: PetscInfo(A, "Option %s ignored\n", MatOptions[op]);
309: break;
310: }
311: return 0;
312: }
314: static PetscErrorCode MatGetRow_MPIAdj(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
315: {
316: Mat_MPIAdj *a = (Mat_MPIAdj *)A->data;
318: row -= A->rmap->rstart;
320: *nz = a->i[row + 1] - a->i[row];
321: if (v) {
322: PetscInt j;
323: if (a->rowvalues_alloc < *nz) {
324: PetscFree(a->rowvalues);
325: a->rowvalues_alloc = PetscMax(a->rowvalues_alloc * 2, *nz);
326: PetscMalloc1(a->rowvalues_alloc, &a->rowvalues);
327: }
328: for (j = 0; j < *nz; j++) a->rowvalues[j] = a->values ? a->values[a->i[row] + j] : 1.0;
329: *v = (*nz) ? a->rowvalues : NULL;
330: }
331: if (idx) *idx = (*nz) ? a->j + a->i[row] : NULL;
332: return 0;
333: }
335: static PetscErrorCode MatRestoreRow_MPIAdj(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
336: {
337: return 0;
338: }
340: static PetscErrorCode MatEqual_MPIAdj(Mat A, Mat B, PetscBool *flg)
341: {
342: Mat_MPIAdj *a = (Mat_MPIAdj *)A->data, *b = (Mat_MPIAdj *)B->data;
343: PetscBool flag;
345: /* If the matrix dimensions are not equal,or no of nonzeros */
346: if ((A->rmap->n != B->rmap->n) || (a->nz != b->nz)) flag = PETSC_FALSE;
348: /* if the a->i are the same */
349: PetscArraycmp(a->i, b->i, A->rmap->n + 1, &flag);
351: /* if a->j are the same */
352: PetscMemcmp(a->j, b->j, (a->nz) * sizeof(PetscInt), &flag);
354: MPIU_Allreduce(&flag, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A));
355: return 0;
356: }
358: static PetscErrorCode MatGetRowIJ_MPIAdj(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *m, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done)
359: {
360: PetscInt i;
361: Mat_MPIAdj *a = (Mat_MPIAdj *)A->data;
362: PetscInt **ia = (PetscInt **)inia, **ja = (PetscInt **)inja;
364: *m = A->rmap->n;
365: *ia = a->i;
366: *ja = a->j;
367: *done = PETSC_TRUE;
368: if (oshift) {
369: for (i = 0; i < (*ia)[*m]; i++) (*ja)[i]++;
370: for (i = 0; i <= (*m); i++) (*ia)[i]++;
371: }
372: return 0;
373: }
375: static PetscErrorCode MatRestoreRowIJ_MPIAdj(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *m, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done)
376: {
377: PetscInt i;
378: Mat_MPIAdj *a = (Mat_MPIAdj *)A->data;
379: PetscInt **ia = (PetscInt **)inia, **ja = (PetscInt **)inja;
383: if (oshift) {
386: for (i = 0; i <= (*m); i++) (*ia)[i]--;
387: for (i = 0; i < (*ia)[*m]; i++) (*ja)[i]--;
388: }
389: return 0;
390: }
392: PetscErrorCode MatConvertFrom_MPIAdj(Mat A, MatType type, MatReuse reuse, Mat *newmat)
393: {
394: Mat B;
395: PetscInt i, m, N, nzeros = 0, *ia, *ja, len, rstart, cnt, j, *a;
396: const PetscInt *rj;
397: const PetscScalar *ra;
398: MPI_Comm comm;
400: MatGetSize(A, NULL, &N);
401: MatGetLocalSize(A, &m, NULL);
402: MatGetOwnershipRange(A, &rstart, NULL);
404: /* count the number of nonzeros per row */
405: for (i = 0; i < m; i++) {
406: MatGetRow(A, i + rstart, &len, &rj, NULL);
407: for (j = 0; j < len; j++) {
408: if (rj[j] == i + rstart) {
409: len--;
410: break;
411: } /* don't count diagonal */
412: }
413: nzeros += len;
414: MatRestoreRow(A, i + rstart, &len, &rj, NULL);
415: }
417: /* malloc space for nonzeros */
418: PetscMalloc1(nzeros + 1, &a);
419: PetscMalloc1(N + 1, &ia);
420: PetscMalloc1(nzeros + 1, &ja);
422: nzeros = 0;
423: ia[0] = 0;
424: for (i = 0; i < m; i++) {
425: MatGetRow(A, i + rstart, &len, &rj, &ra);
426: cnt = 0;
427: for (j = 0; j < len; j++) {
428: if (rj[j] != i + rstart) { /* if not diagonal */
429: a[nzeros + cnt] = (PetscInt)PetscAbsScalar(ra[j]);
430: ja[nzeros + cnt++] = rj[j];
431: }
432: }
433: MatRestoreRow(A, i + rstart, &len, &rj, &ra);
434: nzeros += cnt;
435: ia[i + 1] = nzeros;
436: }
438: PetscObjectGetComm((PetscObject)A, &comm);
439: MatCreate(comm, &B);
440: MatSetSizes(B, m, PETSC_DETERMINE, PETSC_DETERMINE, N);
441: MatSetType(B, type);
442: MatMPIAdjSetPreallocation(B, ia, ja, a);
444: if (reuse == MAT_INPLACE_MATRIX) {
445: MatHeaderReplace(A, &B);
446: } else {
447: *newmat = B;
448: }
449: return 0;
450: }
452: PetscErrorCode MatSetValues_MPIAdj(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode im)
453: {
454: Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data;
455: PetscInt rStart, rEnd, cStart, cEnd;
458: MatGetOwnershipRange(A, &rStart, &rEnd);
459: MatGetOwnershipRangeColumn(A, &cStart, &cEnd);
460: if (!adj->ht) {
461: PetscHSetIJCreate(&adj->ht);
462: MatStashCreate_Private(PetscObjectComm((PetscObject)A), 1, &A->stash);
463: PetscLayoutSetUp(A->rmap);
464: PetscLayoutSetUp(A->cmap);
465: }
466: for (PetscInt r = 0; r < m; ++r) {
467: PetscHashIJKey key;
469: key.i = rows[r];
470: if (key.i < 0) continue;
471: if ((key.i < rStart) || (key.i >= rEnd)) {
472: MatStashValuesRow_Private(&A->stash, key.i, n, cols, values, PETSC_FALSE);
473: } else {
474: for (PetscInt c = 0; c < n; ++c) {
475: key.j = cols[c];
476: if (key.j < 0 || key.i == key.j) continue;
477: PetscHSetIJAdd(adj->ht, key);
478: }
479: }
480: }
481: return 0;
482: }
484: PetscErrorCode MatAssemblyBegin_MPIAdj(Mat A, MatAssemblyType type)
485: {
486: PetscInt nstash, reallocs;
487: Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data;
489: if (!adj->ht) {
490: PetscHSetIJCreate(&adj->ht);
491: PetscLayoutSetUp(A->rmap);
492: PetscLayoutSetUp(A->cmap);
493: }
494: MatStashScatterBegin_Private(A, &A->stash, A->rmap->range);
495: MatStashGetInfo_Private(&A->stash, &nstash, &reallocs);
496: PetscInfo(A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs);
497: return 0;
498: }
500: PetscErrorCode MatAssemblyEnd_MPIAdj(Mat A, MatAssemblyType type)
501: {
502: PetscScalar *val;
503: PetscInt *row, *col, m, rstart, *rowstarts;
504: PetscInt i, j, ncols, flg, nz;
505: PetscMPIInt n;
506: Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data;
507: PetscHashIter hi;
508: PetscHashIJKey key;
509: PetscHSetIJ ht = adj->ht;
511: while (1) {
512: MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg);
513: if (!flg) break;
515: for (i = 0; i < n;) {
516: /* Identify the consecutive vals belonging to the same row */
517: for (j = i, rstart = row[j]; j < n; j++) {
518: if (row[j] != rstart) break;
519: }
520: if (j < n) ncols = j - i;
521: else ncols = n - i;
522: /* Set all these values with a single function call */
523: MatSetValues_MPIAdj(A, 1, row + i, ncols, col + i, val + i, INSERT_VALUES);
524: i = j;
525: }
526: }
527: MatStashScatterEnd_Private(&A->stash);
528: MatStashDestroy_Private(&A->stash);
530: MatGetLocalSize(A, &m, NULL);
531: MatGetOwnershipRange(A, &rstart, NULL);
532: PetscCalloc1(m + 1, &rowstarts);
533: PetscHashIterBegin(ht, hi);
534: for (; !PetscHashIterAtEnd(ht, hi);) {
535: PetscHashIterGetKey(ht, hi, key);
536: rowstarts[key.i - rstart + 1]++;
537: PetscHashIterNext(ht, hi);
538: }
539: for (i = 1; i < m + 1; i++) rowstarts[i] = rowstarts[i - 1] + rowstarts[i];
541: PetscHSetIJGetSize(ht, &nz);
542: PetscMalloc1(nz, &col);
543: PetscHashIterBegin(ht, hi);
544: for (; !PetscHashIterAtEnd(ht, hi);) {
545: PetscHashIterGetKey(ht, hi, key);
546: col[rowstarts[key.i - rstart]++] = key.j;
547: PetscHashIterNext(ht, hi);
548: }
549: PetscHSetIJDestroy(&ht);
550: for (i = m; i > 0; i--) rowstarts[i] = rowstarts[i - 1];
551: rowstarts[0] = 0;
553: for (PetscInt i = 0; i < m; i++) PetscSortInt(rowstarts[i + 1] - rowstarts[i], &col[rowstarts[i]]);
555: adj->i = rowstarts;
556: adj->j = col;
557: adj->nz = rowstarts[m];
558: adj->freeaij = PETSC_TRUE;
559: return 0;
560: }
562: /* -------------------------------------------------------------------*/
563: static struct _MatOps MatOps_Values = {MatSetValues_MPIAdj,
564: MatGetRow_MPIAdj,
565: MatRestoreRow_MPIAdj,
566: NULL,
567: /* 4*/ NULL,
568: NULL,
569: NULL,
570: NULL,
571: NULL,
572: NULL,
573: /*10*/ NULL,
574: NULL,
575: NULL,
576: NULL,
577: NULL,
578: /*15*/ NULL,
579: MatEqual_MPIAdj,
580: NULL,
581: NULL,
582: NULL,
583: /*20*/ MatAssemblyBegin_MPIAdj,
584: MatAssemblyEnd_MPIAdj,
585: MatSetOption_MPIAdj,
586: NULL,
587: /*24*/ NULL,
588: NULL,
589: NULL,
590: NULL,
591: NULL,
592: /*29*/ NULL,
593: NULL,
594: NULL,
595: NULL,
596: NULL,
597: /*34*/ NULL,
598: NULL,
599: NULL,
600: NULL,
601: NULL,
602: /*39*/ NULL,
603: MatCreateSubMatrices_MPIAdj,
604: NULL,
605: NULL,
606: NULL,
607: /*44*/ NULL,
608: NULL,
609: MatShift_Basic,
610: NULL,
611: NULL,
612: /*49*/ NULL,
613: MatGetRowIJ_MPIAdj,
614: MatRestoreRowIJ_MPIAdj,
615: NULL,
616: NULL,
617: /*54*/ NULL,
618: NULL,
619: NULL,
620: NULL,
621: NULL,
622: /*59*/ NULL,
623: MatDestroy_MPIAdj,
624: MatView_MPIAdj,
625: MatConvertFrom_MPIAdj,
626: NULL,
627: /*64*/ NULL,
628: NULL,
629: NULL,
630: NULL,
631: NULL,
632: /*69*/ NULL,
633: NULL,
634: NULL,
635: NULL,
636: NULL,
637: /*74*/ NULL,
638: NULL,
639: NULL,
640: NULL,
641: NULL,
642: /*79*/ NULL,
643: NULL,
644: NULL,
645: NULL,
646: NULL,
647: /*84*/ NULL,
648: NULL,
649: NULL,
650: NULL,
651: NULL,
652: /*89*/ NULL,
653: NULL,
654: NULL,
655: NULL,
656: NULL,
657: /*94*/ NULL,
658: NULL,
659: NULL,
660: NULL,
661: NULL,
662: /*99*/ NULL,
663: NULL,
664: NULL,
665: NULL,
666: NULL,
667: /*104*/ NULL,
668: NULL,
669: NULL,
670: NULL,
671: NULL,
672: /*109*/ NULL,
673: NULL,
674: NULL,
675: NULL,
676: NULL,
677: /*114*/ NULL,
678: NULL,
679: NULL,
680: NULL,
681: NULL,
682: /*119*/ NULL,
683: NULL,
684: NULL,
685: NULL,
686: NULL,
687: /*124*/ NULL,
688: NULL,
689: NULL,
690: NULL,
691: MatCreateSubMatricesMPI_MPIAdj,
692: /*129*/ NULL,
693: NULL,
694: NULL,
695: NULL,
696: NULL,
697: /*134*/ NULL,
698: NULL,
699: NULL,
700: NULL,
701: NULL,
702: /*139*/ NULL,
703: NULL,
704: NULL,
705: NULL,
706: NULL,
707: /*144*/ NULL,
708: NULL,
709: NULL,
710: NULL,
711: NULL,
712: NULL,
713: /*150*/ NULL};
715: static PetscErrorCode MatMPIAdjSetPreallocation_MPIAdj(Mat B, PetscInt *i, PetscInt *j, PetscInt *values)
716: {
717: Mat_MPIAdj *b = (Mat_MPIAdj *)B->data;
718: PetscBool useedgeweights;
720: PetscLayoutSetUp(B->rmap);
721: PetscLayoutSetUp(B->cmap);
722: if (values) useedgeweights = PETSC_TRUE;
723: else useedgeweights = PETSC_FALSE;
724: /* Make everybody knows if they are using edge weights or not */
725: MPIU_Allreduce((int *)&useedgeweights, (int *)&b->useedgeweights, 1, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)B));
727: if (PetscDefined(USE_DEBUG)) {
728: PetscInt ii;
731: for (ii = 1; ii < B->rmap->n; ii++) {
733: }
735: }
736: b->j = j;
737: b->i = i;
738: b->values = values;
740: b->nz = i[B->rmap->n];
741: b->diag = NULL;
742: b->symmetric = PETSC_FALSE;
743: b->freeaij = PETSC_TRUE;
745: B->ops->assemblybegin = NULL;
746: B->ops->assemblyend = NULL;
747: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
748: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
749: MatStashDestroy_Private(&B->stash);
750: return 0;
751: }
753: static PetscErrorCode MatMPIAdjCreateNonemptySubcommMat_MPIAdj(Mat A, Mat *B)
754: {
755: Mat_MPIAdj *a = (Mat_MPIAdj *)A->data;
756: const PetscInt *ranges;
757: MPI_Comm acomm, bcomm;
758: MPI_Group agroup, bgroup;
759: PetscMPIInt i, rank, size, nranks, *ranks;
761: *B = NULL;
762: PetscObjectGetComm((PetscObject)A, &acomm);
763: MPI_Comm_size(acomm, &size);
764: MPI_Comm_size(acomm, &rank);
765: MatGetOwnershipRanges(A, &ranges);
766: for (i = 0, nranks = 0; i < size; i++) {
767: if (ranges[i + 1] - ranges[i] > 0) nranks++;
768: }
769: if (nranks == size) { /* All ranks have a positive number of rows, so we do not need to create a subcomm; */
770: PetscObjectReference((PetscObject)A);
771: *B = A;
772: return 0;
773: }
775: PetscMalloc1(nranks, &ranks);
776: for (i = 0, nranks = 0; i < size; i++) {
777: if (ranges[i + 1] - ranges[i] > 0) ranks[nranks++] = i;
778: }
779: MPI_Comm_group(acomm, &agroup);
780: MPI_Group_incl(agroup, nranks, ranks, &bgroup);
781: PetscFree(ranks);
782: MPI_Comm_create(acomm, bgroup, &bcomm);
783: MPI_Group_free(&agroup);
784: MPI_Group_free(&bgroup);
785: if (bcomm != MPI_COMM_NULL) {
786: PetscInt m, N;
787: Mat_MPIAdj *b;
788: MatGetLocalSize(A, &m, NULL);
789: MatGetSize(A, NULL, &N);
790: MatCreateMPIAdj(bcomm, m, N, a->i, a->j, a->values, B);
791: b = (Mat_MPIAdj *)(*B)->data;
792: b->freeaij = PETSC_FALSE;
793: MPI_Comm_free(&bcomm);
794: }
795: return 0;
796: }
798: PetscErrorCode MatMPIAdjToSeq_MPIAdj(Mat A, Mat *B)
799: {
800: PetscInt M, N, *II, *J, NZ, nz, m, nzstart, i;
801: PetscInt *Values = NULL;
802: Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data;
803: PetscMPIInt mnz, mm, *allnz, *allm, size, *dispnz, *dispm;
805: MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
806: MatGetSize(A, &M, &N);
807: MatGetLocalSize(A, &m, NULL);
808: nz = adj->nz;
810: MPI_Allreduce(&nz, &NZ, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A));
812: PetscMPIIntCast(nz, &mnz);
813: PetscMalloc2(size, &allnz, size, &dispnz);
814: MPI_Allgather(&mnz, 1, MPI_INT, allnz, 1, MPI_INT, PetscObjectComm((PetscObject)A));
815: dispnz[0] = 0;
816: for (i = 1; i < size; i++) dispnz[i] = dispnz[i - 1] + allnz[i - 1];
817: if (adj->values) {
818: PetscMalloc1(NZ, &Values);
819: MPI_Allgatherv(adj->values, mnz, MPIU_INT, Values, allnz, dispnz, MPIU_INT, PetscObjectComm((PetscObject)A));
820: }
821: PetscMalloc1(NZ, &J);
822: MPI_Allgatherv(adj->j, mnz, MPIU_INT, J, allnz, dispnz, MPIU_INT, PetscObjectComm((PetscObject)A));
823: PetscFree2(allnz, dispnz);
824: MPI_Scan(&nz, &nzstart, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A));
825: nzstart -= nz;
826: /* shift the i[] values so they will be correct after being received */
827: for (i = 0; i < m; i++) adj->i[i] += nzstart;
828: PetscMalloc1(M + 1, &II);
829: PetscMPIIntCast(m, &mm);
830: PetscMalloc2(size, &allm, size, &dispm);
831: MPI_Allgather(&mm, 1, MPI_INT, allm, 1, MPI_INT, PetscObjectComm((PetscObject)A));
832: dispm[0] = 0;
833: for (i = 1; i < size; i++) dispm[i] = dispm[i - 1] + allm[i - 1];
834: MPI_Allgatherv(adj->i, mm, MPIU_INT, II, allm, dispm, MPIU_INT, PetscObjectComm((PetscObject)A));
835: PetscFree2(allm, dispm);
836: II[M] = NZ;
837: /* shift the i[] values back */
838: for (i = 0; i < m; i++) adj->i[i] -= nzstart;
839: MatCreateMPIAdj(PETSC_COMM_SELF, M, N, II, J, Values, B);
840: return 0;
841: }
843: PetscErrorCode MatMPIAdjToSeqRankZero_MPIAdj(Mat A, Mat *B)
844: {
845: PetscInt M, N, *II, *J, NZ, nz, m, nzstart, i;
846: PetscInt *Values = NULL;
847: Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data;
848: PetscMPIInt mnz, mm, *allnz = NULL, *allm, size, *dispnz, *dispm, rank;
850: MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
851: MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
852: MatGetSize(A, &M, &N);
853: MatGetLocalSize(A, &m, NULL);
854: nz = adj->nz;
856: MPI_Allreduce(&nz, &NZ, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A));
858: PetscMPIIntCast(nz, &mnz);
859: if (!rank) PetscMalloc2(size, &allnz, size, &dispnz);
860: MPI_Gather(&mnz, 1, MPI_INT, allnz, 1, MPI_INT, 0, PetscObjectComm((PetscObject)A));
861: if (!rank) {
862: dispnz[0] = 0;
863: for (i = 1; i < size; i++) dispnz[i] = dispnz[i - 1] + allnz[i - 1];
864: if (adj->values) {
865: PetscMalloc1(NZ, &Values);
866: MPI_Gatherv(adj->values, mnz, MPIU_INT, Values, allnz, dispnz, MPIU_INT, 0, PetscObjectComm((PetscObject)A));
867: }
868: PetscMalloc1(NZ, &J);
869: MPI_Gatherv(adj->j, mnz, MPIU_INT, J, allnz, dispnz, MPIU_INT, 0, PetscObjectComm((PetscObject)A));
870: PetscFree2(allnz, dispnz);
871: } else {
872: if (adj->values) MPI_Gatherv(adj->values, mnz, MPIU_INT, NULL, NULL, NULL, MPIU_INT, 0, PetscObjectComm((PetscObject)A));
873: MPI_Gatherv(adj->j, mnz, MPIU_INT, NULL, NULL, NULL, MPIU_INT, 0, PetscObjectComm((PetscObject)A));
874: }
875: MPI_Scan(&nz, &nzstart, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A));
876: nzstart -= nz;
877: /* shift the i[] values so they will be correct after being received */
878: for (i = 0; i < m; i++) adj->i[i] += nzstart;
879: PetscMPIIntCast(m, &mm);
880: if (!rank) {
881: PetscMalloc1(M + 1, &II);
882: PetscMalloc2(size, &allm, size, &dispm);
883: MPI_Gather(&mm, 1, MPI_INT, allm, 1, MPI_INT, 0, PetscObjectComm((PetscObject)A));
884: dispm[0] = 0;
885: for (i = 1; i < size; i++) dispm[i] = dispm[i - 1] + allm[i - 1];
886: MPI_Gatherv(adj->i, mm, MPIU_INT, II, allm, dispm, MPIU_INT, 0, PetscObjectComm((PetscObject)A));
887: PetscFree2(allm, dispm);
888: II[M] = NZ;
889: } else {
890: MPI_Gather(&mm, 1, MPI_INT, NULL, 1, MPI_INT, 0, PetscObjectComm((PetscObject)A));
891: MPI_Gatherv(adj->i, mm, MPIU_INT, NULL, NULL, NULL, MPIU_INT, 0, PetscObjectComm((PetscObject)A));
892: }
893: /* shift the i[] values back */
894: for (i = 0; i < m; i++) adj->i[i] -= nzstart;
895: if (!rank) MatCreateMPIAdj(PETSC_COMM_SELF, M, N, II, J, Values, B);
896: return 0;
897: }
899: /*@
900: MatMPIAdjCreateNonemptySubcommMat - create the same `MATMPIADJ` matrix on a subcommunicator containing only processes owning a positive number of rows
902: Collective
904: Input Parameter:
905: . A - original MPIAdj matrix
907: Output Parameter:
908: . B - matrix on subcommunicator, NULL on ranks that owned zero rows of A
910: Level: developer
912: Note:
913: This function is mostly useful for internal use by mesh partitioning packages that require that every process owns at least one row.
915: The matrix B should be destroyed with `MatDestroy()`. The arrays are not copied, so B should be destroyed before A is destroyed.
917: .seealso: `MATMPIADJ`, `MatCreateMPIAdj()`
918: @*/
919: PetscErrorCode MatMPIAdjCreateNonemptySubcommMat(Mat A, Mat *B)
920: {
922: PetscUseMethod(A, "MatMPIAdjCreateNonemptySubcommMat_C", (Mat, Mat *), (A, B));
923: return 0;
924: }
926: /*MC
927: MATMPIADJ - MATMPIADJ = "mpiadj" - A matrix type to be used for distributed adjacency matrices,
928: intended for use constructing orderings and partitionings.
930: Level: beginner
932: Note:
933: You can provide values to the matrix using `MatMPIAdjSetPreallocation()`, `MatCreateMPIAdj()`, or
934: by calling `MatSetValues()` and `MatAssemblyBegin()` followed by `MatAssemblyEnd()`
936: .seealso: `MatCreateMPIAdj()`, `MatMPIAdjSetPreallocation()`, `MatSetValues()`
937: M*/
939: PETSC_EXTERN PetscErrorCode MatCreate_MPIAdj(Mat B)
940: {
941: Mat_MPIAdj *b;
943: PetscNew(&b);
944: B->data = (void *)b;
945: PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps));
946: B->assembled = PETSC_FALSE;
947: B->preallocated = PETSC_TRUE; /* so that MatSetValues() may be used */
949: PetscObjectComposeFunction((PetscObject)B, "MatMPIAdjSetPreallocation_C", MatMPIAdjSetPreallocation_MPIAdj);
950: PetscObjectComposeFunction((PetscObject)B, "MatMPIAdjCreateNonemptySubcommMat_C", MatMPIAdjCreateNonemptySubcommMat_MPIAdj);
951: PetscObjectComposeFunction((PetscObject)B, "MatMPIAdjToSeq_C", MatMPIAdjToSeq_MPIAdj);
952: PetscObjectComposeFunction((PetscObject)B, "MatMPIAdjToSeqRankZero_C", MatMPIAdjToSeqRankZero_MPIAdj);
953: PetscObjectChangeTypeName((PetscObject)B, MATMPIADJ);
954: return 0;
955: }
957: /*@C
958: MatMPIAdjToSeq - Converts an parallel `MATMPIADJ` matrix to complete `MATMPIADJ` on each process (needed by sequential partitioner)
960: Logically Collective
962: Input Parameter:
963: . A - the matrix
965: Output Parameter:
966: . B - the same matrix on all processes
968: Level: intermediate
970: .seealso: `MATMPIADJ`, `MatCreate()`, `MatCreateMPIAdj()`, `MatSetValues()`, `MatMPIAdjToSeqRankZero()`
971: @*/
972: PetscErrorCode MatMPIAdjToSeq(Mat A, Mat *B)
973: {
974: PetscUseMethod(A, "MatMPIAdjToSeq_C", (Mat, Mat *), (A, B));
975: return 0;
976: }
978: /*@C
979: MatMPIAdjToSeqRankZero - Converts an parallel `MATMPIADJ` matrix to complete `MATMPIADJ` on rank zero (needed by sequential partitioner)
981: Logically Collective
983: Input Parameter:
984: . A - the matrix
986: Output Parameter:
987: . B - the same matrix on rank zero, not set on other ranks
989: Level: intermediate
991: Note:
992: This routine has the advantage on systems with multiple ranks per node since only one copy of the matrix
993: is stored on the first node, instead of the number of ranks copies. This can allow partitioning much larger
994: parallel graph sequentially.
996: .seealso: `MATMPIADJ`, `MatCreate()`, `MatCreateMPIAdj()`, `MatSetValues()`, `MatMPIAdjToSeq()`
997: @*/
998: PetscErrorCode MatMPIAdjToSeqRankZero(Mat A, Mat *B)
999: {
1000: PetscUseMethod(A, "MatMPIAdjToSeqRankZero_C", (Mat, Mat *), (A, B));
1001: return 0;
1002: }
1004: /*@C
1005: MatMPIAdjSetPreallocation - Sets the array used for storing the matrix elements
1007: Logically Collective
1009: Input Parameters:
1010: + A - the matrix
1011: . i - the indices into j for the start of each row
1012: . j - the column indices for each row (sorted for each row).
1013: The indices in i and j start with zero (NOT with one).
1014: - values - [optional] edge weights
1016: Level: intermediate
1018: .seealso: `MatCreate()`, `MatCreateMPIAdj()`, `MatSetValues()`, `MATMPIADJ`, `MatCreateMPIAdj()`
1019: @*/
1020: PetscErrorCode MatMPIAdjSetPreallocation(Mat B, PetscInt *i, PetscInt *j, PetscInt *values)
1021: {
1022: PetscTryMethod(B, "MatMPIAdjSetPreallocation_C", (Mat, PetscInt *, PetscInt *, PetscInt *), (B, i, j, values));
1023: return 0;
1024: }
1026: /*@C
1027: MatCreateMPIAdj - Creates a sparse matrix representing an adjacency list.
1028: The matrix does not have numerical values associated with it, but is
1029: intended for ordering (to reduce bandwidth etc) and partitioning.
1031: Collective
1033: Input Parameters:
1034: + comm - MPI communicator
1035: . m - number of local rows
1036: . N - number of global columns
1037: . i - the indices into j for the start of each row
1038: . j - the column indices for each row (sorted for each row).
1039: The indices in i and j start with zero (NOT with one).
1040: - values -[optional] edge weights
1042: Output Parameter:
1043: . A - the matrix
1045: Level: intermediate
1047: Notes:
1048: You must NOT free the ii, values and jj arrays yourself. PETSc will free them
1049: when the matrix is destroyed; you must allocate them with `PetscMalloc()`. If you
1050: call from Fortran you need not create the arrays with `PetscMalloc()`.
1052: You should not include the matrix diagonals.
1054: If you already have a matrix, you can create its adjacency matrix by a call
1055: to `MatConvert()`, specifying a type of `MATMPIADJ`.
1057: Possible values for `MatSetOption()` - `MAT_STRUCTURALLY_SYMMETRIC`
1059: .seealso: `MatCreate()`, `MatConvert()`, `MatGetOrdering()`, `MATMPIADJ`, `MatMPIAdjSetPreallocation()`
1060: @*/
1061: PetscErrorCode MatCreateMPIAdj(MPI_Comm comm, PetscInt m, PetscInt N, PetscInt *i, PetscInt *j, PetscInt *values, Mat *A)
1062: {
1063: MatCreate(comm, A);
1064: MatSetSizes(*A, m, PETSC_DETERMINE, PETSC_DETERMINE, N);
1065: MatSetType(*A, MATMPIADJ);
1066: MatMPIAdjSetPreallocation(*A, i, j, values);
1067: return 0;
1068: }