Actual source code: mpibaij.c
1: #include <../src/mat/impls/baij/mpi/mpibaij.h>
3: #include <petsc/private/hashseti.h>
4: #include <petscblaslapack.h>
5: #include <petscsf.h>
7: #if PetscDefined(HAVE_LIBXSMM)
8: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIBAIJLIBXSMM(Mat, MatType, MatReuse, Mat *);
9: #endif
11: static PetscErrorCode MatDestroy_MPIBAIJ(Mat mat)
12: {
13: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
15: PetscFunctionBegin;
16: PetscCall(PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ",Cols=%" PetscInt_FMT, mat->rmap->N, mat->cmap->N));
17: PetscCall(MatStashDestroy_Private(&mat->stash));
18: PetscCall(MatStashDestroy_Private(&mat->bstash));
19: PetscCall(MatDestroy(&baij->A));
20: PetscCall(MatDestroy(&baij->B));
21: #if PetscDefined(USE_CTABLE)
22: PetscCall(PetscHMapIDestroy(&baij->colmap));
23: #else
24: PetscCall(PetscFree(baij->colmap));
25: #endif
26: PetscCall(PetscFree(baij->garray));
27: PetscCall(VecDestroy(&baij->lvec));
28: PetscCall(VecScatterDestroy(&baij->Mvctx));
29: PetscCall(PetscFree2(baij->rowvalues, baij->rowindices));
30: PetscCall(PetscFree(baij->barray));
31: PetscCall(PetscFree2(baij->hd, baij->ht));
32: PetscCall(PetscFree(baij->rangebs));
33: PetscCall(PetscFree(mat->data));
35: PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
36: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatStoreValues_C", NULL));
37: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatRetrieveValues_C", NULL));
38: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatGetMultPetscSF_C", NULL));
39: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIBAIJSetPreallocation_C", NULL));
40: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIBAIJSetPreallocationCSR_C", NULL));
41: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalScaleLocal_C", NULL));
42: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSetHashTableFactor_C", NULL));
43: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpibaij_mpidense_C", NULL));
44: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_mpisbaij_C", NULL));
45: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_mpiadj_C", NULL));
46: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_mpiaij_C", NULL));
47: #if PetscDefined(HAVE_HYPRE)
48: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_hypre_C", NULL));
49: #endif
50: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_is_C", NULL));
51: #if PetscDefined(HAVE_LIBXSMM)
52: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_mpibaijlibxsmm_C", NULL));
53: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpibaijlibxsmm_mpidense_C", NULL));
54: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaijlibxsmm_mpibaij_C", NULL));
55: #endif
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: /* defines MatSetValues_MPI_Hash(), MatAssemblyBegin_MPI_Hash(), and MatAssemblyEnd_MPI_Hash() */
60: #define TYPE BAIJ
61: #include "../src/mat/impls/aij/mpi/mpihashmat.h"
62: #undef TYPE
64: #if PetscDefined(HAVE_HYPRE)
65: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat, MatType, MatReuse, Mat *);
66: #endif
68: static PetscErrorCode MatGetRowMaxAbs_MPIBAIJ(Mat A, Vec v, PetscInt idx[])
69: {
70: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
71: PetscInt i, *idxb = NULL, m = A->rmap->n, bs = A->cmap->bs;
72: PetscScalar *vv;
73: Vec vB, vA;
74: const PetscScalar *va, *vb;
76: PetscFunctionBegin;
77: PetscCall(MatCreateVecs(a->A, NULL, &vA));
78: PetscCall(MatGetRowMaxAbs(a->A, vA, idx));
80: PetscCall(VecGetArrayRead(vA, &va));
81: if (idx) {
82: for (i = 0; i < m; i++) {
83: if (PetscAbsScalar(va[i])) idx[i] += A->cmap->rstart;
84: }
85: }
87: PetscCall(MatCreateVecs(a->B, NULL, &vB));
88: PetscCall(PetscMalloc1(m, &idxb));
89: PetscCall(MatGetRowMaxAbs(a->B, vB, idxb));
91: PetscCall(VecGetArrayWrite(v, &vv));
92: PetscCall(VecGetArrayRead(vB, &vb));
93: for (i = 0; i < m; i++) {
94: if (PetscAbsScalar(va[i]) < PetscAbsScalar(vb[i])) {
95: vv[i] = vb[i];
96: if (idx) idx[i] = bs * a->garray[idxb[i] / bs] + (idxb[i] % bs);
97: } else {
98: vv[i] = va[i];
99: if (idx && PetscAbsScalar(va[i]) == PetscAbsScalar(vb[i]) && idxb[i] != -1 && idx[i] > bs * a->garray[idxb[i] / bs] + (idxb[i] % bs)) idx[i] = bs * a->garray[idxb[i] / bs] + (idxb[i] % bs);
100: }
101: }
102: PetscCall(VecRestoreArrayWrite(v, &vv));
103: PetscCall(VecRestoreArrayRead(vA, &va));
104: PetscCall(VecRestoreArrayRead(vB, &vb));
105: PetscCall(PetscFree(idxb));
106: PetscCall(VecDestroy(&vA));
107: PetscCall(VecDestroy(&vB));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: static PetscErrorCode MatGetRowSumAbs_MPIBAIJ(Mat A, Vec v)
112: {
113: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
114: Vec vB, vA;
116: PetscFunctionBegin;
117: PetscCall(MatCreateVecs(a->A, NULL, &vA));
118: PetscCall(MatGetRowSumAbs(a->A, vA));
119: PetscCall(MatCreateVecs(a->B, NULL, &vB));
120: PetscCall(MatGetRowSumAbs(a->B, vB));
121: PetscCall(VecAXPY(vA, 1.0, vB));
122: PetscCall(VecDestroy(&vB));
123: PetscCall(VecCopy(vA, v));
124: PetscCall(VecDestroy(&vA));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: static PetscErrorCode MatStoreValues_MPIBAIJ(Mat mat)
129: {
130: Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)mat->data;
132: PetscFunctionBegin;
133: PetscCall(MatStoreValues(aij->A));
134: PetscCall(MatStoreValues(aij->B));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode MatRetrieveValues_MPIBAIJ(Mat mat)
139: {
140: Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)mat->data;
142: PetscFunctionBegin;
143: PetscCall(MatRetrieveValues(aij->A));
144: PetscCall(MatRetrieveValues(aij->B));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*
149: Local utility routine that creates a mapping from the global column
150: number to the local number in the off-diagonal part of the local
151: storage of the matrix. This is done in a non scalable way since the
152: length of colmap equals the global matrix length.
153: */
154: PetscErrorCode MatCreateColmap_MPIBAIJ_Private(Mat mat)
155: {
156: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
157: Mat_SeqBAIJ *B = (Mat_SeqBAIJ *)baij->B->data;
158: PetscInt nbs = B->nbs, i, bs = mat->rmap->bs;
160: PetscFunctionBegin;
161: #if PetscDefined(USE_CTABLE)
162: PetscCall(PetscHMapICreateWithSize(baij->nbs, &baij->colmap));
163: for (i = 0; i < nbs; i++) PetscCall(PetscHMapISet(baij->colmap, baij->garray[i] + 1, i * bs + 1));
164: #else
165: PetscCall(PetscCalloc1(baij->Nbs + 1, &baij->colmap));
166: for (i = 0; i < nbs; i++) baij->colmap[baij->garray[i]] = i * bs + 1;
167: #endif
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: #define MatSetValues_SeqBAIJ_A_Private(row, col, value, addv, orow, ocol) \
172: do { \
173: brow = row / bs; \
174: rp = PetscSafePointerPlusOffset(aj, ai[brow]); \
175: ap = PetscSafePointerPlusOffset(aa, bs2 * ai[brow]); \
176: rmax = aimax[brow]; \
177: nrow = ailen[brow]; \
178: bcol = col / bs; \
179: ridx = row % bs; \
180: cidx = col % bs; \
181: low = 0; \
182: high = nrow; \
183: while (high - low > 3) { \
184: t = (low + high) / 2; \
185: if (rp[t] > bcol) high = t; \
186: else low = t; \
187: } \
188: for (_i = low; _i < high; _i++) { \
189: if (rp[_i] > bcol) break; \
190: if (rp[_i] == bcol) { \
191: bap = ap + bs2 * _i + bs * cidx + ridx; \
192: if (addv == ADD_VALUES) *bap += value; \
193: else *bap = value; \
194: goto a_noinsert; \
195: } \
196: } \
197: if (a->nonew == 1) goto a_noinsert; \
198: PetscCheck(a->nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \
199: MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, aimax, a->nonew, MatScalar); \
200: N = nrow++ - 1; \
201: /* shift up all the later entries in this row */ \
202: PetscCall(PetscArraymove(rp + _i + 1, rp + _i, N - _i + 1)); \
203: PetscCall(PetscArraymove(ap + bs2 * (_i + 1), ap + bs2 * _i, bs2 * (N - _i + 1))); \
204: PetscCall(PetscArrayzero(ap + bs2 * _i, bs2)); \
205: rp[_i] = bcol; \
206: ap[bs2 * _i + bs * cidx + ridx] = value; \
207: a_noinsert:; \
208: ailen[brow] = nrow; \
209: } while (0)
211: #define MatSetValues_SeqBAIJ_B_Private(row, col, value, addv, orow, ocol) \
212: do { \
213: brow = row / bs; \
214: rp = PetscSafePointerPlusOffset(bj, bi[brow]); \
215: ap = PetscSafePointerPlusOffset(ba, bs2 * bi[brow]); \
216: rmax = bimax[brow]; \
217: nrow = bilen[brow]; \
218: bcol = col / bs; \
219: ridx = row % bs; \
220: cidx = col % bs; \
221: low = 0; \
222: high = nrow; \
223: while (high - low > 3) { \
224: t = (low + high) / 2; \
225: if (rp[t] > bcol) high = t; \
226: else low = t; \
227: } \
228: for (_i = low; _i < high; _i++) { \
229: if (rp[_i] > bcol) break; \
230: if (rp[_i] == bcol) { \
231: bap = ap + bs2 * _i + bs * cidx + ridx; \
232: if (addv == ADD_VALUES) *bap += value; \
233: else *bap = value; \
234: goto b_noinsert; \
235: } \
236: } \
237: if (b->nonew == 1) goto b_noinsert; \
238: PetscCheck(b->nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \
239: MatSeqXAIJReallocateAIJ(B, b->mbs, bs2, nrow, brow, bcol, rmax, ba, bi, bj, rp, ap, bimax, b->nonew, MatScalar); \
240: N = nrow++ - 1; \
241: /* shift up all the later entries in this row */ \
242: PetscCall(PetscArraymove(rp + _i + 1, rp + _i, N - _i + 1)); \
243: PetscCall(PetscArraymove(ap + bs2 * (_i + 1), ap + bs2 * _i, bs2 * (N - _i + 1))); \
244: PetscCall(PetscArrayzero(ap + bs2 * _i, bs2)); \
245: rp[_i] = bcol; \
246: ap[bs2 * _i + bs * cidx + ridx] = value; \
247: b_noinsert:; \
248: bilen[brow] = nrow; \
249: } while (0)
251: PetscErrorCode MatSetValues_MPIBAIJ(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
252: {
253: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
254: MatScalar value;
255: PetscBool roworiented = baij->roworiented;
256: PetscInt i, j, row, col;
257: PetscInt rstart_orig = mat->rmap->rstart;
258: PetscInt rend_orig = mat->rmap->rend, cstart_orig = mat->cmap->rstart;
259: PetscInt cend_orig = mat->cmap->rend, bs = mat->rmap->bs;
261: /* Some Variables required in the macro */
262: Mat A = baij->A;
263: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
264: PetscInt *aimax = a->imax, *ai = a->i, *ailen = a->ilen, *aj = a->j;
265: MatScalar *aa = a->a;
267: Mat B = baij->B;
268: Mat_SeqBAIJ *b = (Mat_SeqBAIJ *)B->data;
269: PetscInt *bimax = b->imax, *bi = b->i, *bilen = b->ilen, *bj = b->j;
270: MatScalar *ba = b->a;
272: PetscInt *rp, ii, nrow, _i, rmax, N, brow, bcol;
273: PetscInt low, high, t, ridx, cidx, bs2 = a->bs2;
274: MatScalar *ap, *bap;
276: PetscFunctionBegin;
277: for (i = 0; i < m; i++) {
278: if (im[i] < 0) continue;
279: PetscCheck(im[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], mat->rmap->N - 1);
280: if (im[i] >= rstart_orig && im[i] < rend_orig) {
281: row = im[i] - rstart_orig;
282: for (j = 0; j < n; j++) {
283: if (in[j] >= cstart_orig && in[j] < cend_orig) {
284: col = in[j] - cstart_orig;
285: if (roworiented) value = v[i * n + j];
286: else value = v[i + j * m];
287: MatSetValues_SeqBAIJ_A_Private(row, col, value, addv, im[i], in[j]);
288: } else if (in[j] < 0) {
289: continue;
290: } else {
291: PetscCheck(in[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[j], mat->cmap->N - 1);
292: if (mat->was_assembled) {
293: if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
294: #if PetscDefined(USE_CTABLE)
295: PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] / bs + 1, 0, &col));
296: col = col - 1;
297: #else
298: col = baij->colmap[in[j] / bs] - 1;
299: #endif
300: if (col < 0 && !((Mat_SeqBAIJ *)baij->B->data)->nonew) {
301: PetscCall(MatDisAssemble_MPIBAIJ(mat));
302: col = in[j];
303: /* Reinitialize the variables required by MatSetValues_SeqBAIJ_B_Private() */
304: B = baij->B;
305: b = (Mat_SeqBAIJ *)B->data;
306: bimax = b->imax;
307: bi = b->i;
308: bilen = b->ilen;
309: bj = b->j;
310: ba = b->a;
311: } else {
312: PetscCheck(col >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", im[i], in[j]);
313: col += in[j] % bs;
314: }
315: } else col = in[j];
316: if (roworiented) value = v[i * n + j];
317: else value = v[i + j * m];
318: MatSetValues_SeqBAIJ_B_Private(row, col, value, addv, im[i], in[j]);
319: /* PetscCall(MatSetValues_SeqBAIJ(baij->B,1,&row,1,&col,&value,addv)); */
320: }
321: }
322: } else {
323: PetscCheck(!mat->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", im[i]);
324: if (!baij->donotstash) {
325: mat->assembled = PETSC_FALSE;
326: if (roworiented) {
327: PetscCall(MatStashValuesRow_Private(&mat->stash, im[i], n, in, v + i * n, PETSC_FALSE));
328: } else {
329: PetscCall(MatStashValuesCol_Private(&mat->stash, im[i], n, in, v + i, m, PETSC_FALSE));
330: }
331: }
332: }
333: }
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static inline PetscErrorCode MatSetValuesBlocked_SeqBAIJ_Inlined(Mat A, PetscInt row, PetscInt col, const PetscScalar v[], InsertMode is, PetscInt orow, PetscInt ocol)
338: {
339: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
340: PetscInt *rp, low, high, t, ii, jj, nrow, i, rmax, N;
341: PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen;
342: PetscInt *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs;
343: PetscBool roworiented = a->roworiented;
344: const PetscScalar *value = v;
345: MatScalar *ap, *aa = a->a, *bap;
347: PetscFunctionBegin;
348: rp = aj + ai[row];
349: ap = aa + bs2 * ai[row];
350: rmax = imax[row];
351: nrow = ailen[row];
352: value = v;
353: low = 0;
354: high = nrow;
355: while (high - low > 7) {
356: t = (low + high) / 2;
357: if (rp[t] > col) high = t;
358: else low = t;
359: }
360: for (i = low; i < high; i++) {
361: if (rp[i] > col) break;
362: if (rp[i] == col) {
363: bap = ap + bs2 * i;
364: if (roworiented) {
365: if (is == ADD_VALUES) {
366: for (ii = 0; ii < bs; ii++) {
367: for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
368: }
369: } else {
370: for (ii = 0; ii < bs; ii++) {
371: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
372: }
373: }
374: } else {
375: if (is == ADD_VALUES) {
376: for (ii = 0; ii < bs; ii++, value += bs) {
377: for (jj = 0; jj < bs; jj++) bap[jj] += value[jj];
378: bap += bs;
379: }
380: } else {
381: for (ii = 0; ii < bs; ii++, value += bs) {
382: for (jj = 0; jj < bs; jj++) bap[jj] = value[jj];
383: bap += bs;
384: }
385: }
386: }
387: goto noinsert2;
388: }
389: }
390: if (nonew == 1) goto noinsert2;
391: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new global block indexed nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", orow, ocol);
392: MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
393: N = nrow++ - 1;
394: high++;
395: /* shift up all the later entries in this row */
396: PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
397: PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
398: rp[i] = col;
399: bap = ap + bs2 * i;
400: if (roworiented) {
401: for (ii = 0; ii < bs; ii++) {
402: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
403: }
404: } else {
405: for (ii = 0; ii < bs; ii++) {
406: for (jj = 0; jj < bs; jj++) *bap++ = *value++;
407: }
408: }
409: noinsert2:;
410: ailen[row] = nrow;
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: /*
415: This routine should be optimized so that the block copy at ** Here a copy is required ** below is not needed
416: by passing additional stride information into the MatSetValuesBlocked_SeqBAIJ_Inlined() routine
417: */
418: static PetscErrorCode MatSetValuesBlocked_MPIBAIJ(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
419: {
420: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
421: const PetscScalar *value;
422: MatScalar *barray = baij->barray;
423: PetscBool roworiented = baij->roworiented;
424: PetscInt i, j, ii, jj, row, col, rstart = baij->rstartbs;
425: PetscInt rend = baij->rendbs, cstart = baij->cstartbs, stepval;
426: PetscInt cend = baij->cendbs, bs = mat->rmap->bs, bs2 = baij->bs2;
428: PetscFunctionBegin;
429: if (!barray) {
430: PetscCall(PetscMalloc1(bs2, &barray));
431: baij->barray = barray;
432: }
434: if (roworiented) stepval = (n - 1) * bs;
435: else stepval = (m - 1) * bs;
437: for (i = 0; i < m; i++) {
438: if (im[i] < 0) continue;
439: PetscCheck(im[i] < baij->Mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block indexed row too large %" PetscInt_FMT " max %" PetscInt_FMT, im[i], baij->Mbs - 1);
440: if (im[i] >= rstart && im[i] < rend) {
441: row = im[i] - rstart;
442: for (j = 0; j < n; j++) {
443: /* If NumCol = 1 then a copy is not required */
444: if ((roworiented) && (n == 1)) {
445: barray = (MatScalar *)v + i * bs2;
446: } else if ((!roworiented) && (m == 1)) {
447: barray = (MatScalar *)v + j * bs2;
448: } else { /* Here a copy is required */
449: if (roworiented) {
450: value = v + (i * (stepval + bs) + j) * bs;
451: } else {
452: value = v + (j * (stepval + bs) + i) * bs;
453: }
454: for (ii = 0; ii < bs; ii++, value += bs + stepval) {
455: for (jj = 0; jj < bs; jj++) barray[jj] = value[jj];
456: barray += bs;
457: }
458: barray -= bs2;
459: }
461: if (in[j] >= cstart && in[j] < cend) {
462: col = in[j] - cstart;
463: PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->A, row, col, barray, addv, im[i], in[j]));
464: } else if (in[j] < 0) {
465: continue;
466: } else {
467: PetscCheck(in[j] < baij->Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block indexed column too large %" PetscInt_FMT " max %" PetscInt_FMT, in[j], baij->Nbs - 1);
468: if (mat->was_assembled) {
469: if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
471: #if PetscDefined(USE_CTABLE)
472: PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] + 1, 0, &col));
473: col = col < 1 ? -1 : (col - 1) / bs;
474: #else
475: col = baij->colmap[in[j]] < 1 ? -1 : (baij->colmap[in[j]] - 1) / bs;
476: #endif
477: if (col < 0 && !((Mat_SeqBAIJ *)baij->B->data)->nonew) {
478: PetscCall(MatDisAssemble_MPIBAIJ(mat));
479: col = in[j];
480: } else PetscCheck(col >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new blocked indexed nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", im[i], in[j]);
481: } else col = in[j];
482: PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->B, row, col, barray, addv, im[i], in[j]));
483: }
484: }
485: } else {
486: PetscCheck(!mat->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process block indexed row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", im[i]);
487: if (!baij->donotstash) {
488: if (roworiented) {
489: PetscCall(MatStashValuesRowBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
490: } else {
491: PetscCall(MatStashValuesColBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
492: }
493: }
494: }
495: }
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: #define HASH_KEY 0.6180339887
500: #define HASH(size, key, tmp) (tmp = (key) * HASH_KEY, (PetscInt)((size) * (tmp - (PetscInt)tmp)))
501: /* #define HASH(size,key) ((PetscInt)((size)*fmod(((key)*HASH_KEY),1))) */
502: /* #define HASH(size,key,tmp) ((PetscInt)((size)*fmod(((key)*HASH_KEY),1))) */
503: static PetscErrorCode MatSetValues_MPIBAIJ_HT(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
504: {
505: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
506: PetscBool roworiented = baij->roworiented;
507: PetscInt i, j, row, col;
508: PetscInt rstart_orig = mat->rmap->rstart;
509: PetscInt rend_orig = mat->rmap->rend, Nbs = baij->Nbs;
510: PetscInt h1, key, size = baij->ht_size, bs = mat->rmap->bs, *HT = baij->ht, idx;
511: PetscReal tmp;
512: MatScalar **HD = baij->hd, value;
513: PetscInt total_ct = baij->ht_total_ct, insert_ct = baij->ht_insert_ct;
515: PetscFunctionBegin;
516: for (i = 0; i < m; i++) {
517: if (PetscDefined(USE_DEBUG)) {
518: PetscCheck(im[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row");
519: PetscCheck(im[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], mat->rmap->N - 1);
520: }
521: row = im[i];
522: if (row >= rstart_orig && row < rend_orig) {
523: for (j = 0; j < n; j++) {
524: col = in[j];
525: if (roworiented) value = v[i * n + j];
526: else value = v[i + j * m];
527: /* Look up PetscInto the Hash Table */
528: key = (row / bs) * Nbs + (col / bs) + 1;
529: h1 = HASH(size, key, tmp);
531: idx = h1;
532: if (PetscDefined(USE_DEBUG)) {
533: insert_ct++;
534: total_ct++;
535: if (HT[idx] != key) {
536: for (idx = h1; (idx < size) && (HT[idx] != key); idx++, total_ct++);
537: if (idx == size) {
538: for (idx = 0; (idx < h1) && (HT[idx] != key); idx++, total_ct++);
539: PetscCheck(idx != h1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "(%" PetscInt_FMT ",%" PetscInt_FMT ") has no entry in the hash table", row, col);
540: }
541: }
542: } else if (HT[idx] != key) {
543: for (idx = h1; (idx < size) && (HT[idx] != key); idx++);
544: if (idx == size) {
545: for (idx = 0; (idx < h1) && (HT[idx] != key); idx++);
546: PetscCheck(idx != h1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "(%" PetscInt_FMT ",%" PetscInt_FMT ") has no entry in the hash table", row, col);
547: }
548: }
549: /* A HASH table entry is found, so insert the values at the correct address */
550: if (addv == ADD_VALUES) *(HD[idx] + (col % bs) * bs + (row % bs)) += value;
551: else *(HD[idx] + (col % bs) * bs + (row % bs)) = value;
552: }
553: } else if (!baij->donotstash) {
554: if (roworiented) {
555: PetscCall(MatStashValuesRow_Private(&mat->stash, im[i], n, in, v + i * n, PETSC_FALSE));
556: } else {
557: PetscCall(MatStashValuesCol_Private(&mat->stash, im[i], n, in, v + i, m, PETSC_FALSE));
558: }
559: }
560: }
561: if (PetscDefined(USE_DEBUG)) {
562: baij->ht_total_ct += total_ct;
563: baij->ht_insert_ct += insert_ct;
564: }
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: static PetscErrorCode MatSetValuesBlocked_MPIBAIJ_HT(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
569: {
570: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
571: PetscBool roworiented = baij->roworiented;
572: PetscInt i, j, ii, jj, row, col;
573: PetscInt rstart = baij->rstartbs;
574: PetscInt rend = mat->rmap->rend, stepval, bs = mat->rmap->bs, bs2 = baij->bs2, nbs2 = n * bs2;
575: PetscInt h1, key, size = baij->ht_size, idx, *HT = baij->ht, Nbs = baij->Nbs;
576: PetscReal tmp;
577: MatScalar **HD = baij->hd, *baij_a;
578: const PetscScalar *v_t, *value;
579: PetscInt total_ct = baij->ht_total_ct, insert_ct = baij->ht_insert_ct;
581: PetscFunctionBegin;
582: if (roworiented) stepval = (n - 1) * bs;
583: else stepval = (m - 1) * bs;
585: for (i = 0; i < m; i++) {
586: if (PetscDefined(USE_DEBUG)) {
587: PetscCheck(im[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row: %" PetscInt_FMT, im[i]);
588: PetscCheck(im[i] < baij->Mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], baij->Mbs - 1);
589: }
590: row = im[i];
591: v_t = v + i * nbs2;
592: if (row >= rstart && row < rend) {
593: for (j = 0; j < n; j++) {
594: col = in[j];
596: /* Look up into the Hash Table */
597: key = row * Nbs + col + 1;
598: h1 = HASH(size, key, tmp);
600: idx = h1;
601: if (PetscDefined(USE_DEBUG)) {
602: total_ct++;
603: insert_ct++;
604: if (HT[idx] != key) {
605: for (idx = h1; (idx < size) && (HT[idx] != key); idx++, total_ct++);
606: if (idx == size) {
607: for (idx = 0; (idx < h1) && (HT[idx] != key); idx++, total_ct++);
608: PetscCheck(idx != h1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "(%" PetscInt_FMT ",%" PetscInt_FMT ") has no entry in the hash table", row, col);
609: }
610: }
611: } else if (HT[idx] != key) {
612: for (idx = h1; (idx < size) && (HT[idx] != key); idx++);
613: if (idx == size) {
614: for (idx = 0; (idx < h1) && (HT[idx] != key); idx++);
615: PetscCheck(idx != h1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "(%" PetscInt_FMT ",%" PetscInt_FMT ") has no entry in the hash table", row, col);
616: }
617: }
618: baij_a = HD[idx];
619: if (roworiented) {
620: /*value = v + i*(stepval+bs)*bs + j*bs;*/
621: /* value = v + (i*(stepval+bs)+j)*bs; */
622: value = v_t;
623: v_t += bs;
624: if (addv == ADD_VALUES) {
625: for (ii = 0; ii < bs; ii++, value += stepval) {
626: for (jj = ii; jj < bs2; jj += bs) baij_a[jj] += *value++;
627: }
628: } else {
629: for (ii = 0; ii < bs; ii++, value += stepval) {
630: for (jj = ii; jj < bs2; jj += bs) baij_a[jj] = *value++;
631: }
632: }
633: } else {
634: value = v + j * (stepval + bs) * bs + i * bs;
635: if (addv == ADD_VALUES) {
636: for (ii = 0; ii < bs; ii++, value += stepval, baij_a += bs) {
637: for (jj = 0; jj < bs; jj++) baij_a[jj] += *value++;
638: }
639: } else {
640: for (ii = 0; ii < bs; ii++, value += stepval, baij_a += bs) {
641: for (jj = 0; jj < bs; jj++) baij_a[jj] = *value++;
642: }
643: }
644: }
645: }
646: } else {
647: if (!baij->donotstash) {
648: if (roworiented) {
649: PetscCall(MatStashValuesRowBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
650: } else {
651: PetscCall(MatStashValuesColBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
652: }
653: }
654: }
655: }
656: if (PetscDefined(USE_DEBUG)) {
657: baij->ht_total_ct += total_ct;
658: baij->ht_insert_ct += insert_ct;
659: }
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: static PetscErrorCode MatGetValues_MPIBAIJ(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
664: {
665: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
666: PetscInt bs = mat->rmap->bs, i, j, bsrstart = mat->rmap->rstart, bsrend = mat->rmap->rend;
667: PetscInt bscstart = mat->cmap->rstart, bscend = mat->cmap->rend, row, col, data;
668: PetscBool roworiented = baij->roworiented;
669: PetscScalar *value;
671: PetscFunctionBegin;
672: for (i = 0; i < m; i++) {
673: if (idxm[i] < 0) continue; /* negative row */
674: PetscCheck(idxm[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm[i], mat->rmap->N - 1);
675: PetscCheck(idxm[i] >= bsrstart && idxm[i] < bsrend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local values currently supported");
676: row = idxm[i] - bsrstart;
677: for (j = 0; j < n; j++) {
678: if (idxn[j] < 0) continue; /* negative column */
679: PetscCheck(idxn[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, idxn[j], mat->cmap->N - 1);
680: value = roworiented ? &v[j + i * n] : &v[i + j * m];
681: if (idxn[j] >= bscstart && idxn[j] < bscend) {
682: col = idxn[j] - bscstart;
683: PetscCall(MatGetValues_SeqBAIJ(baij->A, 1, &row, 1, &col, value));
684: } else {
685: if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
686: #if PetscDefined(USE_CTABLE)
687: PetscCall(PetscHMapIGetWithDefault(baij->colmap, idxn[j] / bs + 1, 0, &data));
688: data--;
689: #else
690: data = baij->colmap[idxn[j] / bs] - 1;
691: #endif
692: if (data < 0 || baij->garray[data / bs] != idxn[j] / bs) *value = 0.0;
693: else {
694: col = data + idxn[j] % bs;
695: PetscCall(MatGetValues_SeqBAIJ(baij->B, 1, &row, 1, &col, value));
696: }
697: }
698: }
699: }
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode MatNorm_MPIBAIJ(Mat mat, NormType type, PetscReal *nrm)
704: {
705: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
706: Mat_SeqBAIJ *amat = (Mat_SeqBAIJ *)baij->A->data, *bmat = (Mat_SeqBAIJ *)baij->B->data;
707: PetscInt i, j, bs2 = baij->bs2, bs = baij->A->rmap->bs, nz, row, col;
708: PetscReal sum = 0.0;
709: MatScalar *v;
711: PetscFunctionBegin;
712: if (baij->size == 1) {
713: PetscCall(MatNorm(baij->A, type, nrm));
714: } else {
715: if (type == NORM_FROBENIUS) {
716: v = amat->a;
717: nz = amat->nz * bs2;
718: for (i = 0; i < nz; i++) {
719: sum += PetscRealPart(PetscConj(*v) * (*v));
720: v++;
721: }
722: v = bmat->a;
723: nz = bmat->nz * bs2;
724: for (i = 0; i < nz; i++) {
725: sum += PetscRealPart(PetscConj(*v) * (*v));
726: v++;
727: }
728: PetscCallMPI(MPIU_Allreduce(&sum, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)mat)));
729: *nrm = PetscSqrtReal(*nrm);
730: } else if (type == NORM_1) { /* max column sum */
731: Vec col, bcol;
732: PetscScalar *array;
733: PetscInt *jj, *garray = baij->garray;
735: PetscCall(MatCreateVecs(mat, &col, NULL));
736: PetscCall(VecGetArrayWrite(col, &array));
737: v = amat->a;
738: jj = amat->j;
739: for (i = 0; i < amat->nz; i++) {
740: for (j = 0; j < bs; j++) {
741: PetscInt col = bs * *jj + j; /* column index */
743: for (row = 0; row < bs; row++) array[col] += PetscAbsScalar(*v++);
744: }
745: jj++;
746: }
747: PetscCall(VecRestoreArrayWrite(col, &array));
748: PetscCall(MatCreateVecs(baij->B, &bcol, NULL));
749: PetscCall(VecGetArrayWrite(bcol, &array));
750: v = bmat->a;
751: jj = bmat->j;
752: for (i = 0; i < bmat->nz; i++) {
753: for (j = 0; j < bs; j++) {
754: PetscInt col = bs * *jj + j; /* column index */
756: for (row = 0; row < bs; row++) array[col] += PetscAbsScalar(*v++);
757: }
758: jj++;
759: }
760: PetscCall(VecSetValuesBlocked(col, bmat->nbs, garray, array, ADD_VALUES));
761: PetscCall(VecRestoreArrayWrite(bcol, &array));
762: PetscCall(VecDestroy(&bcol));
763: PetscCall(VecAssemblyBegin(col));
764: PetscCall(VecAssemblyEnd(col));
765: PetscCall(VecNorm(col, NORM_INFINITY, nrm));
766: PetscCall(VecDestroy(&col));
767: } else if (type == NORM_INFINITY) { /* max row sum */
768: PetscReal *sums;
769: PetscCall(PetscMalloc1(bs, &sums));
770: sum = 0.0;
771: for (j = 0; j < amat->mbs; j++) {
772: for (row = 0; row < bs; row++) sums[row] = 0.0;
773: v = amat->a + bs2 * amat->i[j];
774: nz = amat->i[j + 1] - amat->i[j];
775: for (i = 0; i < nz; i++) {
776: for (col = 0; col < bs; col++) {
777: for (row = 0; row < bs; row++) {
778: sums[row] += PetscAbsScalar(*v);
779: v++;
780: }
781: }
782: }
783: v = bmat->a + bs2 * bmat->i[j];
784: nz = bmat->i[j + 1] - bmat->i[j];
785: for (i = 0; i < nz; i++) {
786: for (col = 0; col < bs; col++) {
787: for (row = 0; row < bs; row++) {
788: sums[row] += PetscAbsScalar(*v);
789: v++;
790: }
791: }
792: }
793: for (row = 0; row < bs; row++) {
794: if (sums[row] > sum) sum = sums[row];
795: }
796: }
797: PetscCallMPI(MPIU_Allreduce(&sum, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)mat)));
798: PetscCall(PetscFree(sums));
799: } else SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support for this norm yet");
800: }
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: /*
805: Creates the hash table, and sets the table
806: This table is created only once.
807: If new entries need to be added to the matrix
808: then the hash table has to be destroyed and
809: recreated.
810: */
811: static PetscErrorCode MatCreateHashTable_MPIBAIJ_Private(Mat mat, PetscReal factor)
812: {
813: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
814: Mat A = baij->A, B = baij->B;
815: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)B->data;
816: PetscInt i, j, k, nz = a->nz + b->nz, h1, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
817: PetscInt ht_size, bs2 = baij->bs2, rstart = baij->rstartbs;
818: PetscInt cstart = baij->cstartbs, *garray = baij->garray, row, col, Nbs = baij->Nbs;
819: PetscInt *HT, key;
820: MatScalar **HD;
821: PetscReal tmp;
822: PetscInt ct = 0, max = 0;
824: PetscFunctionBegin;
825: if (baij->ht) PetscFunctionReturn(PETSC_SUCCESS);
827: baij->ht_size = (PetscInt)(factor * nz);
828: ht_size = baij->ht_size;
830: /* Allocate Memory for Hash Table */
831: PetscCall(PetscCalloc2(ht_size, &baij->hd, ht_size, &baij->ht));
832: HD = baij->hd;
833: HT = baij->ht;
835: /* Loop Over A */
836: for (i = 0; i < a->mbs; i++) {
837: for (j = ai[i]; j < ai[i + 1]; j++) {
838: row = i + rstart;
839: col = aj[j] + cstart;
841: key = row * Nbs + col + 1;
842: h1 = HASH(ht_size, key, tmp);
843: for (k = 0; k < ht_size; k++) {
844: if (!HT[(h1 + k) % ht_size]) {
845: HT[(h1 + k) % ht_size] = key;
846: HD[(h1 + k) % ht_size] = a->a + j * bs2;
847: break;
848: } else if (PetscDefined(USE_INFO)) ct++;
849: }
850: if (PetscDefined(USE_INFO) && k > max) max = k;
851: }
852: }
853: /* Loop Over B */
854: for (i = 0; i < b->mbs; i++) {
855: for (j = bi[i]; j < bi[i + 1]; j++) {
856: row = i + rstart;
857: col = garray[bj[j]];
858: key = row * Nbs + col + 1;
859: h1 = HASH(ht_size, key, tmp);
860: for (k = 0; k < ht_size; k++) {
861: if (!HT[(h1 + k) % ht_size]) {
862: HT[(h1 + k) % ht_size] = key;
863: HD[(h1 + k) % ht_size] = b->a + j * bs2;
864: break;
865: } else if (PetscDefined(USE_INFO)) ct++;
866: }
867: if (PetscDefined(USE_INFO) && k > max) max = k;
868: }
869: }
871: /* Print Summary */
872: if (PetscDefined(USE_INFO)) {
873: for (i = 0, j = 0; i < ht_size; i++) {
874: if (HT[i]) j++;
875: }
876: PetscCall(PetscInfo(mat, "Average Search = %5.2g,max search = %" PetscInt_FMT "\n", (!j) ? 0.0 : (double)(((PetscReal)(ct + j)) / j), max));
877: }
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: static PetscErrorCode MatAssemblyBegin_MPIBAIJ(Mat mat, MatAssemblyType mode)
882: {
883: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
884: PetscInt nstash, reallocs;
886: PetscFunctionBegin;
887: if (baij->donotstash || mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);
889: PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
890: PetscCall(MatStashScatterBegin_Private(mat, &mat->bstash, baij->rangebs));
891: PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs));
892: PetscCall(PetscInfo(mat, "Stash has %" PetscInt_FMT " entries,uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
893: PetscCall(MatStashGetInfo_Private(&mat->bstash, &nstash, &reallocs));
894: PetscCall(PetscInfo(mat, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: static PetscErrorCode MatAssemblyEnd_MPIBAIJ(Mat mat, MatAssemblyType mode)
899: {
900: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
901: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)baij->A->data;
902: PetscInt i, j, rstart, ncols, flg, bs2 = baij->bs2;
903: PetscInt *row, *col;
904: PetscBool r1, r2, r3, all_assembled;
905: MatScalar *val;
906: PetscMPIInt n;
908: PetscFunctionBegin;
909: /* do not use 'b=(Mat_SeqBAIJ*)baij->B->data' as B can be reset in disassembly */
910: if (!baij->donotstash && !mat->nooffprocentries) {
911: while (1) {
912: PetscCall(MatStashScatterGetMesg_Private(&mat->stash, &n, &row, &col, &val, &flg));
913: if (!flg) break;
915: for (i = 0; i < n;) {
916: /* Now identify the consecutive vals belonging to the same row */
917: for (j = i, rstart = row[j]; j < n; j++) {
918: if (row[j] != rstart) break;
919: }
920: if (j < n) ncols = j - i;
921: else ncols = n - i;
922: /* Now assemble all these values with a single function call */
923: PetscCall(MatSetValues_MPIBAIJ(mat, 1, row + i, ncols, col + i, val + i, mat->insertmode));
924: i = j;
925: }
926: }
927: PetscCall(MatStashScatterEnd_Private(&mat->stash));
928: /* Now process the block-stash. Since the values are stashed column-oriented,
929: set the row-oriented flag to column-oriented, and after MatSetValues()
930: restore the original flags */
931: r1 = baij->roworiented;
932: r2 = a->roworiented;
933: r3 = ((Mat_SeqBAIJ *)baij->B->data)->roworiented;
935: baij->roworiented = PETSC_FALSE;
936: a->roworiented = PETSC_FALSE;
937: ((Mat_SeqBAIJ *)baij->B->data)->roworiented = PETSC_FALSE;
938: while (1) {
939: PetscCall(MatStashScatterGetMesg_Private(&mat->bstash, &n, &row, &col, &val, &flg));
940: if (!flg) break;
942: for (i = 0; i < n;) {
943: /* Now identify the consecutive vals belonging to the same row */
944: for (j = i, rstart = row[j]; j < n; j++) {
945: if (row[j] != rstart) break;
946: }
947: if (j < n) ncols = j - i;
948: else ncols = n - i;
949: PetscCall(MatSetValuesBlocked_MPIBAIJ(mat, 1, row + i, ncols, col + i, val + i * bs2, mat->insertmode));
950: i = j;
951: }
952: }
953: PetscCall(MatStashScatterEnd_Private(&mat->bstash));
955: baij->roworiented = r1;
956: a->roworiented = r2;
957: ((Mat_SeqBAIJ *)baij->B->data)->roworiented = r3;
958: }
960: PetscCall(MatAssemblyBegin(baij->A, mode));
961: PetscCall(MatAssemblyEnd(baij->A, mode));
963: /* determine if any process has disassembled, if so we must
964: also disassemble ourselves, in order that we may reassemble. */
965: /*
966: if nonzero structure of submatrix B cannot change then we know that
967: no process disassembled thus we can skip this stuff
968: */
969: if (!((Mat_SeqBAIJ *)baij->B->data)->nonew) {
970: PetscCallMPI(MPIU_Allreduce(&mat->was_assembled, &all_assembled, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
971: if (mat->was_assembled && !all_assembled) PetscCall(MatDisAssemble_MPIBAIJ(mat));
972: }
974: if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) PetscCall(MatSetUpMultiply_MPIBAIJ(mat));
975: PetscCall(MatAssemblyBegin(baij->B, mode));
976: PetscCall(MatAssemblyEnd(baij->B, mode));
978: if (PetscDefined(USE_INFO) && baij->ht && mode == MAT_FINAL_ASSEMBLY) {
979: PetscCall(PetscInfo(mat, "Average Hash Table Search in MatSetValues = %5.2f\n", (double)((PetscReal)baij->ht_total_ct) / baij->ht_insert_ct));
981: baij->ht_total_ct = 0;
982: baij->ht_insert_ct = 0;
983: }
984: if (baij->ht_flag && !baij->ht && mode == MAT_FINAL_ASSEMBLY) {
985: PetscCall(MatCreateHashTable_MPIBAIJ_Private(mat, baij->ht_fact));
987: mat->ops->setvalues = MatSetValues_MPIBAIJ_HT;
988: mat->ops->setvaluesblocked = MatSetValuesBlocked_MPIBAIJ_HT;
989: }
991: PetscCall(PetscFree2(baij->rowvalues, baij->rowindices));
993: baij->rowvalues = NULL;
995: /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
996: if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqBAIJ *)baij->A->data)->nonew) {
997: mat->nonzerostate = baij->A->nonzerostate + baij->B->nonzerostate;
998: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &mat->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)mat)));
999: }
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: #include <petscdraw.h>
1004: static PetscErrorCode MatView_MPIBAIJ_ASCIIorDraworSocket(Mat mat, PetscViewer viewer)
1005: {
1006: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
1007: PetscMPIInt rank = baij->rank;
1008: PetscInt bs = mat->rmap->bs;
1009: PetscBool isascii, isdraw;
1010: PetscViewer sviewer;
1011: PetscViewerFormat format;
1013: PetscFunctionBegin;
1014: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1015: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1016: if (isascii) {
1017: PetscCall(PetscViewerGetFormat(viewer, &format));
1018: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1019: MatInfo info;
1020: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank));
1021: PetscCall(MatGetInfo(mat, MAT_LOCAL, &info));
1022: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1023: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " bs %" PetscInt_FMT " mem %g\n", rank, mat->rmap->n, (PetscInt)info.nz_used, (PetscInt)info.nz_allocated,
1024: mat->rmap->bs, info.memory));
1025: PetscCall(MatGetInfo(baij->A, MAT_LOCAL, &info));
1026: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] on-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
1027: PetscCall(MatGetInfo(baij->B, MAT_LOCAL, &info));
1028: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] off-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
1029: PetscCall(PetscViewerFlush(viewer));
1030: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1031: PetscCall(PetscViewerASCIIPrintf(viewer, "Information on VecScatter used in matrix-vector product: \n"));
1032: PetscCall(VecScatterView(baij->Mvctx, viewer));
1033: PetscFunctionReturn(PETSC_SUCCESS);
1034: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_FACTOR_INFO) PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: if (isdraw) {
1038: PetscDraw draw;
1039: PetscBool isnull;
1040: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1041: PetscCall(PetscDrawIsNull(draw, &isnull));
1042: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: {
1046: /* assemble the entire matrix onto first processor. */
1047: Mat A;
1048: Mat_SeqBAIJ *Aloc;
1049: PetscInt M = mat->rmap->N, N = mat->cmap->N, *ai, *aj, col, i, j, k, *rvals, mbs = baij->mbs;
1050: MatScalar *a;
1051: const char *matname;
1053: /* Here we are creating a temporary matrix, so will assume MPIBAIJ is acceptable */
1054: /* Perhaps this should be the type of mat? */
1055: PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &A));
1056: if (rank == 0) {
1057: PetscCall(MatSetSizes(A, M, N, M, N));
1058: } else {
1059: PetscCall(MatSetSizes(A, 0, 0, M, N));
1060: }
1061: PetscCall(MatSetType(A, MATMPIBAIJ));
1062: PetscCall(MatMPIBAIJSetPreallocation(A, mat->rmap->bs, 0, NULL, 0, NULL));
1063: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
1065: /* copy over the A part */
1066: Aloc = (Mat_SeqBAIJ *)baij->A->data;
1067: ai = Aloc->i;
1068: aj = Aloc->j;
1069: a = Aloc->a;
1070: PetscCall(PetscMalloc1(bs, &rvals));
1072: for (i = 0; i < mbs; i++) {
1073: rvals[0] = bs * (baij->rstartbs + i);
1074: for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
1075: for (j = ai[i]; j < ai[i + 1]; j++) {
1076: col = (baij->cstartbs + aj[j]) * bs;
1077: for (k = 0; k < bs; k++) {
1078: PetscCall(MatSetValues_MPIBAIJ(A, bs, rvals, 1, &col, a, INSERT_VALUES));
1079: col++;
1080: a += bs;
1081: }
1082: }
1083: }
1084: /* copy over the B part */
1085: Aloc = (Mat_SeqBAIJ *)baij->B->data;
1086: ai = Aloc->i;
1087: aj = Aloc->j;
1088: a = Aloc->a;
1089: for (i = 0; i < mbs; i++) {
1090: rvals[0] = bs * (baij->rstartbs + i);
1091: for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
1092: for (j = ai[i]; j < ai[i + 1]; j++) {
1093: col = baij->garray[aj[j]] * bs;
1094: for (k = 0; k < bs; k++) {
1095: PetscCall(MatSetValues_MPIBAIJ(A, bs, rvals, 1, &col, a, INSERT_VALUES));
1096: col++;
1097: a += bs;
1098: }
1099: }
1100: }
1101: PetscCall(PetscFree(rvals));
1102: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1103: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1104: /*
1105: Everyone has to call to draw the matrix since the graphics waits are
1106: synchronized across all processors that share the PetscDraw object
1107: */
1108: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
1109: if (((PetscObject)mat)->name) PetscCall(PetscObjectGetName((PetscObject)mat, &matname));
1110: if (rank == 0) {
1111: if (((PetscObject)mat)->name) PetscCall(PetscObjectSetName((PetscObject)((Mat_MPIBAIJ *)A->data)->A, matname));
1112: PetscCall(MatView_SeqBAIJ(((Mat_MPIBAIJ *)A->data)->A, sviewer));
1113: }
1114: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
1115: PetscCall(MatDestroy(&A));
1116: }
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: /* Used for both MPIBAIJ and MPISBAIJ matrices */
1121: PetscErrorCode MatView_MPIBAIJ_Binary(Mat mat, PetscViewer viewer)
1122: {
1123: Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)mat->data;
1124: Mat_SeqBAIJ *A = (Mat_SeqBAIJ *)aij->A->data;
1125: Mat_SeqBAIJ *B = (Mat_SeqBAIJ *)aij->B->data;
1126: const PetscInt *garray = aij->garray;
1127: PetscInt header[4], M, N, m, rs, cs, bs, cnt, i, j, ja, jb, k, l;
1128: PetscCount nz, hnz;
1129: PetscInt *rowlens, *colidxs;
1130: PetscScalar *matvals;
1131: PetscMPIInt rank;
1133: PetscFunctionBegin;
1134: PetscCall(PetscViewerSetUp(viewer));
1136: M = mat->rmap->N;
1137: N = mat->cmap->N;
1138: m = mat->rmap->n;
1139: rs = mat->rmap->rstart;
1140: cs = mat->cmap->rstart;
1141: bs = mat->rmap->bs;
1142: nz = bs * bs * (A->nz + B->nz);
1144: /* write matrix header */
1145: header[0] = MAT_FILE_CLASSID;
1146: header[1] = M;
1147: header[2] = N;
1148: PetscCallMPI(MPI_Reduce(&nz, &hnz, 1, MPIU_COUNT, MPI_SUM, 0, PetscObjectComm((PetscObject)mat)));
1149: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank));
1150: if (rank == 0) PetscCall(PetscIntCast(hnz, &header[3]));
1151: PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));
1153: /* fill in and store row lengths */
1154: PetscCall(PetscMalloc1(m, &rowlens));
1155: for (cnt = 0, i = 0; i < A->mbs; i++)
1156: for (j = 0; j < bs; j++) rowlens[cnt++] = bs * (A->i[i + 1] - A->i[i] + B->i[i + 1] - B->i[i]);
1157: PetscCall(PetscViewerBinaryWriteAll(viewer, rowlens, m, rs, M, PETSC_INT));
1158: PetscCall(PetscFree(rowlens));
1160: /* fill in and store column indices */
1161: PetscCall(PetscMalloc1(nz, &colidxs));
1162: for (cnt = 0, i = 0; i < A->mbs; i++) {
1163: for (k = 0; k < bs; k++) {
1164: for (jb = B->i[i]; jb < B->i[i + 1]; jb++) {
1165: if (garray[B->j[jb]] > cs / bs) break;
1166: for (l = 0; l < bs; l++) colidxs[cnt++] = bs * garray[B->j[jb]] + l;
1167: }
1168: for (ja = A->i[i]; ja < A->i[i + 1]; ja++)
1169: for (l = 0; l < bs; l++) colidxs[cnt++] = bs * A->j[ja] + l + cs;
1170: for (; jb < B->i[i + 1]; jb++)
1171: for (l = 0; l < bs; l++) colidxs[cnt++] = bs * garray[B->j[jb]] + l;
1172: }
1173: }
1174: PetscCheck(cnt == nz, PETSC_COMM_SELF, PETSC_ERR_LIB, "Internal PETSc error: cnt = %" PetscInt_FMT " nz = %" PetscCount_FMT, cnt, nz);
1175: PetscCall(PetscViewerBinaryWriteAll(viewer, colidxs, nz, PETSC_DECIDE, PETSC_DECIDE, PETSC_INT));
1176: PetscCall(PetscFree(colidxs));
1178: /* fill in and store nonzero values */
1179: PetscCall(PetscMalloc1(nz, &matvals));
1180: for (cnt = 0, i = 0; i < A->mbs; i++) {
1181: for (k = 0; k < bs; k++) {
1182: for (jb = B->i[i]; jb < B->i[i + 1]; jb++) {
1183: if (garray[B->j[jb]] > cs / bs) break;
1184: for (l = 0; l < bs; l++) matvals[cnt++] = B->a[bs * (bs * jb + l) + k];
1185: }
1186: for (ja = A->i[i]; ja < A->i[i + 1]; ja++)
1187: for (l = 0; l < bs; l++) matvals[cnt++] = A->a[bs * (bs * ja + l) + k];
1188: for (; jb < B->i[i + 1]; jb++)
1189: for (l = 0; l < bs; l++) matvals[cnt++] = B->a[bs * (bs * jb + l) + k];
1190: }
1191: }
1192: PetscCall(PetscViewerBinaryWriteAll(viewer, matvals, nz, PETSC_DECIDE, PETSC_DECIDE, PETSC_SCALAR));
1193: PetscCall(PetscFree(matvals));
1195: /* write block size option to the viewer's .info file */
1196: PetscCall(MatView_Binary_BlockSizes(mat, viewer));
1197: PetscFunctionReturn(PETSC_SUCCESS);
1198: }
1200: PetscErrorCode MatView_MPIBAIJ(Mat mat, PetscViewer viewer)
1201: {
1202: PetscBool isascii, isdraw, issocket, isbinary;
1204: PetscFunctionBegin;
1205: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1206: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1207: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
1208: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1209: if (isascii || isdraw || issocket) PetscCall(MatView_MPIBAIJ_ASCIIorDraworSocket(mat, viewer));
1210: else if (isbinary) PetscCall(MatView_MPIBAIJ_Binary(mat, viewer));
1211: PetscFunctionReturn(PETSC_SUCCESS);
1212: }
1214: static PetscErrorCode MatMult_MPIBAIJ(Mat A, Vec xx, Vec yy)
1215: {
1216: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1217: PetscInt nt;
1219: PetscFunctionBegin;
1220: PetscCall(VecGetLocalSize(xx, &nt));
1221: PetscCheck(nt == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible partition of A and xx");
1222: PetscCall(VecGetLocalSize(yy, &nt));
1223: PetscCheck(nt == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible partition of A and yy");
1224: PetscCall(VecScatterBegin(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
1225: PetscUseTypeMethod(a->A, mult, xx, yy);
1226: PetscCall(VecScatterEnd(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
1227: PetscUseTypeMethod(a->B, multadd, a->lvec, yy, yy);
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: static PetscErrorCode MatMultAdd_MPIBAIJ(Mat A, Vec xx, Vec yy, Vec zz)
1232: {
1233: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1235: PetscFunctionBegin;
1236: PetscCall(VecScatterBegin(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
1237: PetscUseTypeMethod(a->A, multadd, xx, yy, zz);
1238: PetscCall(VecScatterEnd(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
1239: PetscUseTypeMethod(a->B, multadd, a->lvec, zz, zz);
1240: PetscFunctionReturn(PETSC_SUCCESS);
1241: }
1243: static PetscErrorCode MatMultTranspose_MPIBAIJ(Mat A, Vec xx, Vec yy)
1244: {
1245: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1247: PetscFunctionBegin;
1248: /* do nondiagonal part */
1249: PetscUseTypeMethod(a->B, multtranspose, xx, a->lvec);
1250: /* do local part */
1251: PetscUseTypeMethod(a->A, multtranspose, xx, yy);
1252: /* add partial results together */
1253: PetscCall(VecScatterBegin(a->Mvctx, a->lvec, yy, ADD_VALUES, SCATTER_REVERSE));
1254: PetscCall(VecScatterEnd(a->Mvctx, a->lvec, yy, ADD_VALUES, SCATTER_REVERSE));
1255: PetscFunctionReturn(PETSC_SUCCESS);
1256: }
1258: static PetscErrorCode MatMultTransposeAdd_MPIBAIJ(Mat A, Vec xx, Vec yy, Vec zz)
1259: {
1260: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1262: PetscFunctionBegin;
1263: /* do nondiagonal part */
1264: PetscUseTypeMethod(a->B, multtranspose, xx, a->lvec);
1265: /* do local part */
1266: PetscUseTypeMethod(a->A, multtransposeadd, xx, yy, zz);
1267: /* add partial results together */
1268: PetscCall(VecScatterBegin(a->Mvctx, a->lvec, zz, ADD_VALUES, SCATTER_REVERSE));
1269: PetscCall(VecScatterEnd(a->Mvctx, a->lvec, zz, ADD_VALUES, SCATTER_REVERSE));
1270: PetscFunctionReturn(PETSC_SUCCESS);
1271: }
1273: /*
1274: This only works correctly for square matrices where the subblock A->A is the
1275: diagonal block
1276: */
1277: static PetscErrorCode MatGetDiagonal_MPIBAIJ(Mat A, Vec v)
1278: {
1279: PetscFunctionBegin;
1280: PetscCheck(A->rmap->N == A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_SUP, "Supports only square matrix where A->A is diag block");
1281: PetscCall(MatGetDiagonal(((Mat_MPIBAIJ *)A->data)->A, v));
1282: PetscFunctionReturn(PETSC_SUCCESS);
1283: }
1285: static PetscErrorCode MatScale_MPIBAIJ(Mat A, PetscScalar aa)
1286: {
1287: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1289: PetscFunctionBegin;
1290: PetscCall(MatScale(a->A, aa));
1291: PetscCall(MatScale(a->B, aa));
1292: PetscFunctionReturn(PETSC_SUCCESS);
1293: }
1295: static PetscErrorCode MatGetRow_MPIBAIJ(Mat matin, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1296: {
1297: Mat_MPIBAIJ *mat = (Mat_MPIBAIJ *)matin->data;
1298: PetscScalar *vworkA, *vworkB, **pvA, **pvB, *v_p;
1299: PetscInt bs = matin->rmap->bs, bs2 = mat->bs2, i, *cworkA, *cworkB, **pcA, **pcB;
1300: PetscInt nztot, nzA, nzB, lrow, brstart = matin->rmap->rstart, brend = matin->rmap->rend;
1301: PetscInt *cmap, *idx_p, cstart = mat->cstartbs;
1303: PetscFunctionBegin;
1304: PetscCheck(row >= brstart && row < brend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local rows");
1305: PetscCheck(!mat->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Already active");
1306: mat->getrowactive = PETSC_TRUE;
1308: if (!mat->rowvalues && (idx || v)) {
1309: /*
1310: allocate enough space to hold information from the longest row.
1311: */
1312: Mat_SeqBAIJ *Aa = (Mat_SeqBAIJ *)mat->A->data, *Ba = (Mat_SeqBAIJ *)mat->B->data;
1313: PetscInt max = 1, mbs = mat->mbs, tmp;
1314: for (i = 0; i < mbs; i++) {
1315: tmp = Aa->i[i + 1] - Aa->i[i] + Ba->i[i + 1] - Ba->i[i];
1316: if (max < tmp) max = tmp;
1317: }
1318: PetscCall(PetscMalloc2(max * bs2, &mat->rowvalues, max * bs2, &mat->rowindices));
1319: }
1320: lrow = row - brstart;
1322: pvA = &vworkA;
1323: pcA = &cworkA;
1324: pvB = &vworkB;
1325: pcB = &cworkB;
1326: if (!v) {
1327: pvA = NULL;
1328: pvB = NULL;
1329: }
1330: if (!idx) {
1331: pcA = NULL;
1332: if (!v) pcB = NULL;
1333: }
1334: PetscUseTypeMethod(mat->A, getrow, lrow, &nzA, pcA, pvA);
1335: PetscUseTypeMethod(mat->B, getrow, lrow, &nzB, pcB, pvB);
1336: nztot = nzA + nzB;
1338: cmap = mat->garray;
1339: if (v || idx) {
1340: if (nztot) {
1341: /* Sort by increasing column numbers, assuming A and B already sorted */
1342: PetscInt imark = -1;
1343: if (v) {
1344: *v = v_p = mat->rowvalues;
1345: for (i = 0; i < nzB; i++) {
1346: if (cmap[cworkB[i] / bs] < cstart) v_p[i] = vworkB[i];
1347: else break;
1348: }
1349: imark = i;
1350: for (i = 0; i < nzA; i++) v_p[imark + i] = vworkA[i];
1351: for (i = imark; i < nzB; i++) v_p[nzA + i] = vworkB[i];
1352: }
1353: if (idx) {
1354: *idx = idx_p = mat->rowindices;
1355: if (imark > -1) {
1356: for (i = 0; i < imark; i++) idx_p[i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1357: } else {
1358: for (i = 0; i < nzB; i++) {
1359: if (cmap[cworkB[i] / bs] < cstart) idx_p[i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1360: else break;
1361: }
1362: imark = i;
1363: }
1364: for (i = 0; i < nzA; i++) idx_p[imark + i] = cstart * bs + cworkA[i];
1365: for (i = imark; i < nzB; i++) idx_p[nzA + i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1366: }
1367: } else {
1368: if (idx) *idx = NULL;
1369: if (v) *v = NULL;
1370: }
1371: }
1372: *nz = nztot;
1373: PetscUseTypeMethod(mat->A, restorerow, lrow, &nzA, pcA, pvA);
1374: PetscUseTypeMethod(mat->B, restorerow, lrow, &nzB, pcB, pvB);
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: static PetscErrorCode MatRestoreRow_MPIBAIJ(Mat mat, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1379: {
1380: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
1382: PetscFunctionBegin;
1383: PetscCheck(baij->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "MatGetRow not called");
1384: baij->getrowactive = PETSC_FALSE;
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: static PetscErrorCode MatZeroEntries_MPIBAIJ(Mat A)
1389: {
1390: Mat_MPIBAIJ *l = (Mat_MPIBAIJ *)A->data;
1392: PetscFunctionBegin;
1393: PetscCall(MatZeroEntries(l->A));
1394: PetscCall(MatZeroEntries(l->B));
1395: PetscFunctionReturn(PETSC_SUCCESS);
1396: }
1398: static PetscErrorCode MatGetInfo_MPIBAIJ(Mat matin, MatInfoType flag, MatInfo *info)
1399: {
1400: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)matin->data;
1401: Mat A = a->A, B = a->B;
1402: PetscLogDouble irecv[5];
1404: PetscFunctionBegin;
1405: info->block_size = (PetscReal)matin->rmap->bs;
1407: PetscCall(MatGetInfo(A, MAT_LOCAL, info));
1409: irecv[0] = info->nz_used;
1410: irecv[1] = info->nz_allocated;
1411: irecv[2] = info->nz_unneeded;
1412: irecv[3] = info->memory;
1413: irecv[4] = info->mallocs;
1415: PetscCall(MatGetInfo(B, MAT_LOCAL, info));
1417: irecv[0] += info->nz_used;
1418: irecv[1] += info->nz_allocated;
1419: irecv[2] += info->nz_unneeded;
1420: irecv[3] += info->memory;
1421: irecv[4] += info->mallocs;
1423: if (flag == MAT_LOCAL) {
1424: info->nz_used = irecv[0];
1425: info->nz_allocated = irecv[1];
1426: info->nz_unneeded = irecv[2];
1427: info->memory = irecv[3];
1428: info->mallocs = irecv[4];
1429: } else if (flag == MAT_GLOBAL_MAX) {
1430: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)matin)));
1432: info->nz_used = irecv[0];
1433: info->nz_allocated = irecv[1];
1434: info->nz_unneeded = irecv[2];
1435: info->memory = irecv[3];
1436: info->mallocs = irecv[4];
1437: } else if (flag == MAT_GLOBAL_SUM) {
1438: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)matin)));
1440: info->nz_used = irecv[0];
1441: info->nz_allocated = irecv[1];
1442: info->nz_unneeded = irecv[2];
1443: info->memory = irecv[3];
1444: info->mallocs = irecv[4];
1445: } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_ARG_WRONG, "Unknown MatInfoType argument %d", (int)flag);
1446: info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */
1447: info->fill_ratio_needed = 0;
1448: info->factor_mallocs = 0;
1449: PetscFunctionReturn(PETSC_SUCCESS);
1450: }
1452: static PetscErrorCode MatSetOption_MPIBAIJ(Mat A, MatOption op, PetscBool flg)
1453: {
1454: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1456: PetscFunctionBegin;
1457: switch (op) {
1458: case MAT_NEW_NONZERO_LOCATIONS:
1459: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1460: case MAT_UNUSED_NONZERO_LOCATION_ERR:
1461: case MAT_KEEP_NONZERO_PATTERN:
1462: case MAT_NEW_NONZERO_LOCATION_ERR:
1463: MatCheckPreallocated(A, 1);
1464: PetscCall(MatSetOption(a->A, op, flg));
1465: PetscCall(MatSetOption(a->B, op, flg));
1466: break;
1467: case MAT_ROW_ORIENTED:
1468: MatCheckPreallocated(A, 1);
1469: a->roworiented = flg;
1471: PetscCall(MatSetOption(a->A, op, flg));
1472: PetscCall(MatSetOption(a->B, op, flg));
1473: break;
1474: case MAT_IGNORE_OFF_PROC_ENTRIES:
1475: a->donotstash = flg;
1476: break;
1477: case MAT_USE_HASH_TABLE:
1478: a->ht_flag = flg;
1479: a->ht_fact = 1.39;
1480: break;
1481: case MAT_SPD:
1482: case MAT_SYMMETRIC:
1483: case MAT_STRUCTURALLY_SYMMETRIC:
1484: case MAT_HERMITIAN:
1485: case MAT_SYMMETRY_ETERNAL:
1486: case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
1487: case MAT_SPD_ETERNAL:
1488: /* if the diagonal matrix is square it inherits some of the properties above */
1489: if (a->A && A->rmap->n == A->cmap->n) PetscCall(MatSetOption(a->A, op, flg));
1490: break;
1491: default:
1492: break;
1493: }
1494: PetscFunctionReturn(PETSC_SUCCESS);
1495: }
1497: static PetscErrorCode MatTranspose_MPIBAIJ(Mat A, MatReuse reuse, Mat *matout)
1498: {
1499: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)A->data;
1500: Mat_SeqBAIJ *Aloc;
1501: Mat B;
1502: PetscInt M = A->rmap->N, N = A->cmap->N, *ai, *aj, i, *rvals, j, k, col;
1503: PetscInt bs = A->rmap->bs, mbs = baij->mbs;
1504: MatScalar *a;
1506: PetscFunctionBegin;
1507: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1508: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1509: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1510: PetscCall(MatSetSizes(B, A->cmap->n, A->rmap->n, N, M));
1511: PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
1512: /* Do not know preallocation information, but must set block size */
1513: PetscCall(MatMPIBAIJSetPreallocation(B, A->rmap->bs, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL));
1514: } else {
1515: B = *matout;
1516: }
1518: /* copy over the A part */
1519: Aloc = (Mat_SeqBAIJ *)baij->A->data;
1520: ai = Aloc->i;
1521: aj = Aloc->j;
1522: a = Aloc->a;
1523: PetscCall(PetscMalloc1(bs, &rvals));
1525: for (i = 0; i < mbs; i++) {
1526: rvals[0] = bs * (baij->rstartbs + i);
1527: for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
1528: for (j = ai[i]; j < ai[i + 1]; j++) {
1529: col = (baij->cstartbs + aj[j]) * bs;
1530: for (k = 0; k < bs; k++) {
1531: PetscCall(MatSetValues_MPIBAIJ(B, 1, &col, bs, rvals, a, INSERT_VALUES));
1533: col++;
1534: a += bs;
1535: }
1536: }
1537: }
1538: /* copy over the B part */
1539: Aloc = (Mat_SeqBAIJ *)baij->B->data;
1540: ai = Aloc->i;
1541: aj = Aloc->j;
1542: a = Aloc->a;
1543: for (i = 0; i < mbs; i++) {
1544: rvals[0] = bs * (baij->rstartbs + i);
1545: for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
1546: for (j = ai[i]; j < ai[i + 1]; j++) {
1547: col = baij->garray[aj[j]] * bs;
1548: for (k = 0; k < bs; k++) {
1549: PetscCall(MatSetValues_MPIBAIJ(B, 1, &col, bs, rvals, a, INSERT_VALUES));
1550: col++;
1551: a += bs;
1552: }
1553: }
1554: }
1555: PetscCall(PetscFree(rvals));
1556: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1557: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1559: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = B;
1560: else PetscCall(MatHeaderMerge(A, &B));
1561: PetscFunctionReturn(PETSC_SUCCESS);
1562: }
1564: static PetscErrorCode MatDiagonalScale_MPIBAIJ(Mat mat, Vec ll, Vec rr)
1565: {
1566: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
1567: Mat a = baij->A, b = baij->B;
1568: PetscInt s1, s2, s3;
1570: PetscFunctionBegin;
1571: PetscCall(MatGetLocalSize(mat, &s2, &s3));
1572: if (rr) {
1573: PetscCall(VecGetLocalSize(rr, &s1));
1574: PetscCheck(s1 == s3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "right vector non-conforming local size");
1575: /* Overlap communication with computation. */
1576: PetscCall(VecScatterBegin(baij->Mvctx, rr, baij->lvec, INSERT_VALUES, SCATTER_FORWARD));
1577: }
1578: if (ll) {
1579: PetscCall(VecGetLocalSize(ll, &s1));
1580: PetscCheck(s1 == s2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "left vector non-conforming local size");
1581: PetscUseTypeMethod(b, diagonalscale, ll, NULL);
1582: }
1583: /* scale the diagonal block */
1584: PetscUseTypeMethod(a, diagonalscale, ll, rr);
1586: if (rr) {
1587: /* Do a scatter end and then right scale the off-diagonal block */
1588: PetscCall(VecScatterEnd(baij->Mvctx, rr, baij->lvec, INSERT_VALUES, SCATTER_FORWARD));
1589: PetscUseTypeMethod(b, diagonalscale, NULL, baij->lvec);
1590: }
1591: PetscFunctionReturn(PETSC_SUCCESS);
1592: }
1594: static PetscErrorCode MatZeroRows_MPIBAIJ(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1595: {
1596: Mat_MPIBAIJ *l = (Mat_MPIBAIJ *)A->data;
1597: PetscInt *lrows;
1598: PetscInt r, len;
1599: PetscBool cong;
1601: PetscFunctionBegin;
1602: /* get locally owned rows */
1603: PetscCall(MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows));
1604: /* fix right-hand side if needed */
1605: if (x && b) {
1606: const PetscScalar *xx;
1607: PetscScalar *bb;
1609: PetscCall(VecGetArrayRead(x, &xx));
1610: PetscCall(VecGetArray(b, &bb));
1611: for (r = 0; r < len; ++r) bb[lrows[r]] = diag * xx[lrows[r]];
1612: PetscCall(VecRestoreArrayRead(x, &xx));
1613: PetscCall(VecRestoreArray(b, &bb));
1614: }
1616: /* actually zap the local rows */
1617: /*
1618: Zero the required rows. If the "diagonal block" of the matrix
1619: is square and the user wishes to set the diagonal we use separate
1620: code so that MatSetValues() is not called for each diagonal allocating
1621: new memory, thus calling lots of mallocs and slowing things down.
1623: */
1624: /* must zero l->B before l->A because the (diag) case below may put values into l->B*/
1625: PetscCall(MatZeroRows_SeqBAIJ(l->B, len, lrows, 0.0, NULL, NULL));
1626: PetscCall(MatHasCongruentLayouts(A, &cong));
1627: if ((diag != 0.0) && cong) {
1628: PetscCall(MatZeroRows_SeqBAIJ(l->A, len, lrows, diag, NULL, NULL));
1629: } else if (diag != 0.0) {
1630: PetscCall(MatZeroRows_SeqBAIJ(l->A, len, lrows, 0.0, NULL, NULL));
1631: PetscCheck(!((Mat_SeqBAIJ *)l->A->data)->nonew, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatZeroRows() on rectangular matrices cannot be used with the Mat options MAT_NEW_NONZERO_LOCATIONS, MAT_NEW_NONZERO_LOCATION_ERR, and MAT_NEW_NONZERO_ALLOCATION_ERR");
1632: for (r = 0; r < len; ++r) {
1633: const PetscInt row = lrows[r] + A->rmap->rstart;
1634: PetscCall(MatSetValues(A, 1, &row, 1, &row, &diag, INSERT_VALUES));
1635: }
1636: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1637: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1638: } else {
1639: PetscCall(MatZeroRows_SeqBAIJ(l->A, len, lrows, 0.0, NULL, NULL));
1640: }
1641: PetscCall(PetscFree(lrows));
1643: /* only change matrix nonzero state if pattern was allowed to be changed */
1644: if (!((Mat_SeqBAIJ *)l->A->data)->keepnonzeropattern || !((Mat_SeqBAIJ *)l->A->data)->nonew) {
1645: A->nonzerostate = l->A->nonzerostate + l->B->nonzerostate;
1646: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A)));
1647: }
1648: PetscFunctionReturn(PETSC_SUCCESS);
1649: }
1651: static PetscErrorCode MatZeroRowsColumns_MPIBAIJ(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1652: {
1653: Mat_MPIBAIJ *l = (Mat_MPIBAIJ *)A->data;
1654: PetscMPIInt n, p = 0;
1655: PetscInt i, j, k, r, len = 0, row, col, count;
1656: PetscInt *lrows, *owners = A->rmap->range;
1657: PetscSFNode *rrows;
1658: PetscSF sf;
1659: const PetscScalar *xx;
1660: PetscScalar *bb, *mask;
1661: Vec xmask, lmask;
1662: Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)l->B->data;
1663: PetscInt bs = A->rmap->bs, bs2 = baij->bs2;
1664: PetscScalar *aa;
1666: PetscFunctionBegin;
1667: PetscCall(PetscMPIIntCast(A->rmap->n, &n));
1668: /* create PetscSF where leaves are input rows and roots are owned rows */
1669: PetscCall(PetscMalloc1(n, &lrows));
1670: for (r = 0; r < n; ++r) lrows[r] = -1;
1671: PetscCall(PetscMalloc1(N, &rrows));
1672: for (r = 0; r < N; ++r) {
1673: const PetscInt idx = rows[r];
1674: PetscCheck(idx >= 0 && A->rmap->N > idx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " out of range [0,%" PetscInt_FMT ")", idx, A->rmap->N);
1675: if (idx < owners[p] || owners[p + 1] <= idx) { /* short-circuit the search if the last p owns this row too */
1676: PetscCall(PetscLayoutFindOwner(A->rmap, idx, &p));
1677: }
1678: rrows[r].rank = p;
1679: rrows[r].index = rows[r] - owners[p];
1680: }
1681: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
1682: PetscCall(PetscSFSetGraph(sf, n, N, NULL, PETSC_OWN_POINTER, rrows, PETSC_OWN_POINTER));
1683: /* collect flags for rows to be zeroed */
1684: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, (PetscInt *)rows, lrows, MPI_LOR));
1685: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, (PetscInt *)rows, lrows, MPI_LOR));
1686: PetscCall(PetscSFDestroy(&sf));
1687: /* compress and put in row numbers */
1688: for (r = 0; r < n; ++r)
1689: if (lrows[r] >= 0) lrows[len++] = r;
1690: /* zero diagonal part of matrix */
1691: PetscCall(MatZeroRowsColumns(l->A, len, lrows, diag, x, b));
1692: /* handle off-diagonal part of matrix */
1693: PetscCall(MatCreateVecs(A, &xmask, NULL));
1694: PetscCall(VecDuplicate(l->lvec, &lmask));
1695: PetscCall(VecGetArray(xmask, &bb));
1696: for (i = 0; i < len; i++) bb[lrows[i]] = 1;
1697: PetscCall(VecRestoreArray(xmask, &bb));
1698: PetscCall(VecScatterBegin(l->Mvctx, xmask, lmask, ADD_VALUES, SCATTER_FORWARD));
1699: PetscCall(VecScatterEnd(l->Mvctx, xmask, lmask, ADD_VALUES, SCATTER_FORWARD));
1700: PetscCall(VecDestroy(&xmask));
1701: if (x) {
1702: PetscCall(VecScatterBegin(l->Mvctx, x, l->lvec, INSERT_VALUES, SCATTER_FORWARD));
1703: PetscCall(VecScatterEnd(l->Mvctx, x, l->lvec, INSERT_VALUES, SCATTER_FORWARD));
1704: PetscCall(VecGetArrayRead(l->lvec, &xx));
1705: PetscCall(VecGetArray(b, &bb));
1706: }
1707: PetscCall(VecGetArray(lmask, &mask));
1708: /* remove zeroed rows of off-diagonal matrix */
1709: for (i = 0; i < len; ++i) {
1710: row = lrows[i];
1711: count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
1712: aa = baij->a + baij->i[row / bs] * bs2 + (row % bs);
1713: for (k = 0; k < count; ++k) {
1714: aa[0] = 0.0;
1715: aa += bs;
1716: }
1717: }
1718: /* loop over all elements of off process part of matrix zeroing removed columns */
1719: for (i = 0; i < l->B->rmap->N; ++i) {
1720: row = i / bs;
1721: for (j = baij->i[row]; j < baij->i[row + 1]; ++j) {
1722: for (k = 0; k < bs; ++k) {
1723: col = bs * baij->j[j] + k;
1724: if (PetscAbsScalar(mask[col])) {
1725: aa = baij->a + j * bs2 + (i % bs) + bs * k;
1726: if (x) bb[i] -= aa[0] * xx[col];
1727: aa[0] = 0.0;
1728: }
1729: }
1730: }
1731: }
1732: if (x) {
1733: PetscCall(VecRestoreArray(b, &bb));
1734: PetscCall(VecRestoreArrayRead(l->lvec, &xx));
1735: }
1736: PetscCall(VecRestoreArray(lmask, &mask));
1737: PetscCall(VecDestroy(&lmask));
1738: PetscCall(PetscFree(lrows));
1740: /* only change matrix nonzero state if pattern was allowed to be changed */
1741: if (!((Mat_SeqBAIJ *)l->A->data)->nonew) {
1742: A->nonzerostate = l->A->nonzerostate + l->B->nonzerostate;
1743: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A)));
1744: }
1745: PetscFunctionReturn(PETSC_SUCCESS);
1746: }
1748: static PetscErrorCode MatSetUnfactored_MPIBAIJ(Mat A)
1749: {
1750: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1752: PetscFunctionBegin;
1753: PetscCall(MatSetUnfactored(a->A));
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }
1757: static PetscErrorCode MatDuplicate_MPIBAIJ(Mat, MatDuplicateOption, Mat *);
1759: static PetscErrorCode MatEqual_MPIBAIJ(Mat A, Mat B, PetscBool *flag)
1760: {
1761: Mat_MPIBAIJ *matB = (Mat_MPIBAIJ *)B->data, *matA = (Mat_MPIBAIJ *)A->data;
1762: Mat a, b, c, d;
1764: PetscFunctionBegin;
1765: a = matA->A;
1766: b = matA->B;
1767: c = matB->A;
1768: d = matB->B;
1770: PetscCall(MatEqual(a, c, flag));
1771: if (*flag) PetscCall(MatEqual(b, d, flag));
1772: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, flag, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
1773: PetscFunctionReturn(PETSC_SUCCESS);
1774: }
1776: static PetscErrorCode MatCopy_MPIBAIJ(Mat A, Mat B, MatStructure str)
1777: {
1778: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1779: Mat_MPIBAIJ *b = (Mat_MPIBAIJ *)B->data;
1781: PetscFunctionBegin;
1782: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1783: if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
1784: PetscCall(MatCopy_Basic(A, B, str));
1785: } else {
1786: PetscCall(MatCopy(a->A, b->A, str));
1787: PetscCall(MatCopy(a->B, b->B, str));
1788: }
1789: PetscCall(PetscObjectStateIncrease((PetscObject)B));
1790: PetscFunctionReturn(PETSC_SUCCESS);
1791: }
1793: PetscErrorCode MatAXPYGetPreallocation_MPIBAIJ(Mat Y, const PetscInt *yltog, Mat X, const PetscInt *xltog, PetscInt *nnz)
1794: {
1795: PetscInt bs = Y->rmap->bs, m = Y->rmap->N / bs;
1796: Mat_SeqBAIJ *x = (Mat_SeqBAIJ *)X->data;
1797: Mat_SeqBAIJ *y = (Mat_SeqBAIJ *)Y->data;
1799: PetscFunctionBegin;
1800: PetscCall(MatAXPYGetPreallocation_MPIX_private(m, x->i, x->j, xltog, y->i, y->j, yltog, nnz));
1801: PetscFunctionReturn(PETSC_SUCCESS);
1802: }
1804: static PetscErrorCode MatAXPY_MPIBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str)
1805: {
1806: Mat_MPIBAIJ *xx = (Mat_MPIBAIJ *)X->data, *yy = (Mat_MPIBAIJ *)Y->data;
1807: PetscBLASInt bnz, one = 1;
1808: Mat_SeqBAIJ *x, *y;
1809: PetscInt bs2 = Y->rmap->bs * Y->rmap->bs;
1811: PetscFunctionBegin;
1812: if (str == SAME_NONZERO_PATTERN) {
1813: PetscScalar alpha = a;
1814: x = (Mat_SeqBAIJ *)xx->A->data;
1815: y = (Mat_SeqBAIJ *)yy->A->data;
1816: PetscCall(PetscBLASIntCast(x->nz * bs2, &bnz));
1817: PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one));
1818: x = (Mat_SeqBAIJ *)xx->B->data;
1819: y = (Mat_SeqBAIJ *)yy->B->data;
1820: PetscCall(PetscBLASIntCast(x->nz * bs2, &bnz));
1821: PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one));
1822: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
1823: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1824: PetscCall(MatAXPY_Basic(Y, a, X, str));
1825: } else {
1826: Mat B;
1827: PetscInt *nnz_d, *nnz_o, bs = Y->rmap->bs;
1828: PetscCall(PetscMalloc1(yy->A->rmap->N, &nnz_d));
1829: PetscCall(PetscMalloc1(yy->B->rmap->N, &nnz_o));
1830: PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B));
1831: PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name));
1832: PetscCall(MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N));
1833: PetscCall(MatSetBlockSizesFromMats(B, Y, Y));
1834: PetscCall(MatSetType(B, MATMPIBAIJ));
1835: PetscCall(MatAXPYGetPreallocation_SeqBAIJ(yy->A, xx->A, nnz_d));
1836: PetscCall(MatAXPYGetPreallocation_MPIBAIJ(yy->B, yy->garray, xx->B, xx->garray, nnz_o));
1837: PetscCall(MatMPIBAIJSetPreallocation(B, bs, 0, nnz_d, 0, nnz_o));
1838: /* MatAXPY_BasicWithPreallocation() for BAIJ matrix is much slower than AIJ, even for bs=1 ! */
1839: PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
1840: PetscCall(MatHeaderMerge(Y, &B));
1841: PetscCall(PetscFree(nnz_d));
1842: PetscCall(PetscFree(nnz_o));
1843: }
1844: PetscFunctionReturn(PETSC_SUCCESS);
1845: }
1847: static PetscErrorCode MatConjugate_MPIBAIJ(Mat mat)
1848: {
1849: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)mat->data;
1851: PetscFunctionBegin;
1852: PetscCall(MatConjugate_SeqBAIJ(a->A));
1853: PetscCall(MatConjugate_SeqBAIJ(a->B));
1854: PetscFunctionReturn(PETSC_SUCCESS);
1855: }
1857: static PetscErrorCode MatRealPart_MPIBAIJ(Mat A)
1858: {
1859: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1861: PetscFunctionBegin;
1862: PetscCall(MatRealPart(a->A));
1863: PetscCall(MatRealPart(a->B));
1864: PetscFunctionReturn(PETSC_SUCCESS);
1865: }
1867: static PetscErrorCode MatImaginaryPart_MPIBAIJ(Mat A)
1868: {
1869: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1871: PetscFunctionBegin;
1872: PetscCall(MatImaginaryPart(a->A));
1873: PetscCall(MatImaginaryPart(a->B));
1874: PetscFunctionReturn(PETSC_SUCCESS);
1875: }
1877: static PetscErrorCode MatCreateSubMatrix_MPIBAIJ(Mat mat, IS isrow, IS iscol, MatReuse call, Mat *newmat)
1878: {
1879: IS iscol_local;
1880: PetscInt csize;
1882: PetscFunctionBegin;
1883: PetscCall(ISGetLocalSize(iscol, &csize));
1884: if (call == MAT_REUSE_MATRIX) {
1885: PetscCall(PetscObjectQuery((PetscObject)*newmat, "ISAllGather", (PetscObject *)&iscol_local));
1886: PetscCheck(iscol_local, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Submatrix passed in was not used before, cannot reuse");
1887: } else {
1888: PetscCall(ISAllGather(iscol, &iscol_local));
1889: }
1890: PetscCall(MatCreateSubMatrix_MPIBAIJ_Private(mat, isrow, iscol_local, csize, call, newmat, PETSC_FALSE));
1891: if (call == MAT_INITIAL_MATRIX) {
1892: PetscCall(PetscObjectCompose((PetscObject)*newmat, "ISAllGather", (PetscObject)iscol_local));
1893: PetscCall(ISDestroy(&iscol_local));
1894: }
1895: PetscFunctionReturn(PETSC_SUCCESS);
1896: }
1898: /*
1899: Not great since it makes two copies of the submatrix, first an SeqBAIJ
1900: in local and then by concatenating the local matrices the end result.
1901: Writing it directly would be much like MatCreateSubMatrices_MPIBAIJ().
1902: This routine is used for BAIJ and SBAIJ matrices (unfortunate dependency).
1903: */
1904: PetscErrorCode MatCreateSubMatrix_MPIBAIJ_Private(Mat mat, IS isrow, IS iscol, PetscInt csize, MatReuse call, Mat *newmat, PetscBool sym)
1905: {
1906: PetscMPIInt rank, size;
1907: PetscInt i, m, n, rstart, row, rend, nz, *cwork, j, bs;
1908: PetscInt *ii, *jj, nlocal, *dlens, *olens, dlen, olen, jend, mglobal;
1909: Mat M, Mreuse;
1910: MatScalar *vwork, *aa;
1911: MPI_Comm comm;
1912: IS isrow_new, iscol_new;
1913: Mat_SeqBAIJ *aij;
1915: PetscFunctionBegin;
1916: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
1917: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1918: PetscCallMPI(MPI_Comm_size(comm, &size));
1919: /* The compression and expansion should be avoided. Doesn't point
1920: out errors, might change the indices, hence buggey */
1921: PetscCall(ISCompressIndicesGeneral(mat->rmap->N, mat->rmap->n, mat->rmap->bs, 1, &isrow, &isrow_new));
1922: if (isrow == iscol) {
1923: iscol_new = isrow_new;
1924: PetscCall(PetscObjectReference((PetscObject)iscol_new));
1925: } else PetscCall(ISCompressIndicesGeneral(mat->cmap->N, mat->cmap->n, mat->cmap->bs, 1, &iscol, &iscol_new));
1927: if (call == MAT_REUSE_MATRIX) {
1928: PetscCall(PetscObjectQuery((PetscObject)*newmat, "SubMatrix", (PetscObject *)&Mreuse));
1929: PetscCheck(Mreuse, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Submatrix passed in was not used before, cannot reuse");
1930: PetscCall(MatCreateSubMatrices_MPIBAIJ_local(mat, 1, &isrow_new, &iscol_new, MAT_REUSE_MATRIX, &Mreuse, sym));
1931: } else {
1932: PetscCall(MatCreateSubMatrices_MPIBAIJ_local(mat, 1, &isrow_new, &iscol_new, MAT_INITIAL_MATRIX, &Mreuse, sym));
1933: }
1934: PetscCall(ISDestroy(&isrow_new));
1935: PetscCall(ISDestroy(&iscol_new));
1936: /*
1937: m - number of local rows
1938: n - number of columns (same on all processors)
1939: rstart - first row in new global matrix generated
1940: */
1941: PetscCall(MatGetBlockSize(mat, &bs));
1942: PetscCall(MatGetSize(Mreuse, &m, &n));
1943: m = m / bs;
1944: n = n / bs;
1946: if (call == MAT_INITIAL_MATRIX) {
1947: aij = (Mat_SeqBAIJ *)Mreuse->data;
1948: ii = aij->i;
1949: jj = aij->j;
1951: /*
1952: Determine the number of non-zeros in the diagonal and off-diagonal
1953: portions of the matrix in order to do correct preallocation
1954: */
1956: /* first get start and end of "diagonal" columns */
1957: if (csize == PETSC_DECIDE) {
1958: PetscCall(ISGetSize(isrow, &mglobal));
1959: if (mglobal == n * bs) { /* square matrix */
1960: nlocal = m;
1961: } else {
1962: nlocal = n / size + ((n % size) > rank);
1963: }
1964: } else {
1965: nlocal = csize / bs;
1966: }
1967: PetscCallMPI(MPI_Scan(&nlocal, &rend, 1, MPIU_INT, MPI_SUM, comm));
1968: rstart = rend - nlocal;
1969: PetscCheck(rank != size - 1 || rend == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local column sizes %" PetscInt_FMT " do not add up to total number of columns %" PetscInt_FMT, rend, n);
1971: /* next, compute all the lengths */
1972: PetscCall(PetscMalloc2(m + 1, &dlens, m + 1, &olens));
1973: for (i = 0; i < m; i++) {
1974: jend = ii[i + 1] - ii[i];
1975: olen = 0;
1976: dlen = 0;
1977: for (j = 0; j < jend; j++) {
1978: if (*jj < rstart || *jj >= rend) olen++;
1979: else dlen++;
1980: jj++;
1981: }
1982: olens[i] = olen;
1983: dlens[i] = dlen;
1984: }
1985: PetscCall(MatCreate(comm, &M));
1986: PetscCall(MatSetSizes(M, bs * m, bs * nlocal, PETSC_DECIDE, bs * n));
1987: PetscCall(MatSetType(M, sym ? ((PetscObject)mat)->type_name : MATMPIBAIJ));
1988: PetscCall(MatMPIBAIJSetPreallocation(M, bs, 0, dlens, 0, olens));
1989: PetscCall(MatMPISBAIJSetPreallocation(M, bs, 0, dlens, 0, olens));
1990: PetscCall(PetscFree2(dlens, olens));
1991: } else {
1992: PetscInt ml, nl;
1994: M = *newmat;
1995: PetscCall(MatGetLocalSize(M, &ml, &nl));
1996: PetscCheck(ml == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Previous matrix must be same size/layout as request");
1997: PetscCall(MatZeroEntries(M));
1998: /*
1999: The next two lines are needed so we may call MatSetValues_MPIAIJ() below directly,
2000: rather than the slower MatSetValues().
2001: */
2002: M->was_assembled = PETSC_TRUE;
2003: M->assembled = PETSC_FALSE;
2004: }
2005: PetscCall(MatSetOption(M, MAT_ROW_ORIENTED, PETSC_FALSE));
2006: PetscCall(MatGetOwnershipRange(M, &rstart, &rend));
2007: aij = (Mat_SeqBAIJ *)Mreuse->data;
2008: ii = aij->i;
2009: jj = aij->j;
2010: aa = aij->a;
2011: for (i = 0; i < m; i++) {
2012: row = rstart / bs + i;
2013: nz = ii[i + 1] - ii[i];
2014: cwork = jj;
2015: jj = PetscSafePointerPlusOffset(jj, nz);
2016: vwork = aa;
2017: aa = PetscSafePointerPlusOffset(aa, nz * bs * bs);
2018: PetscCall(MatSetValuesBlocked_MPIBAIJ(M, 1, &row, nz, cwork, vwork, INSERT_VALUES));
2019: }
2021: PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
2022: PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
2023: *newmat = M;
2025: /* save submatrix used in processor for next request */
2026: if (call == MAT_INITIAL_MATRIX) {
2027: PetscCall(PetscObjectCompose((PetscObject)M, "SubMatrix", (PetscObject)Mreuse));
2028: PetscCall(PetscObjectDereference((PetscObject)Mreuse));
2029: }
2030: PetscFunctionReturn(PETSC_SUCCESS);
2031: }
2033: static PetscErrorCode MatPermute_MPIBAIJ(Mat A, IS rowp, IS colp, Mat *B)
2034: {
2035: MPI_Comm comm, pcomm;
2036: PetscInt clocal_size, nrows;
2037: const PetscInt *rows;
2038: PetscMPIInt size;
2039: IS crowp, lcolp;
2041: PetscFunctionBegin;
2042: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
2043: /* make a collective version of 'rowp' */
2044: PetscCall(PetscObjectGetComm((PetscObject)rowp, &pcomm));
2045: if (pcomm == comm) {
2046: crowp = rowp;
2047: } else {
2048: PetscCall(ISGetSize(rowp, &nrows));
2049: PetscCall(ISGetIndices(rowp, &rows));
2050: PetscCall(ISCreateGeneral(comm, nrows, rows, PETSC_COPY_VALUES, &crowp));
2051: PetscCall(ISRestoreIndices(rowp, &rows));
2052: }
2053: PetscCall(ISSetPermutation(crowp));
2054: /* make a local version of 'colp' */
2055: PetscCall(PetscObjectGetComm((PetscObject)colp, &pcomm));
2056: PetscCallMPI(MPI_Comm_size(pcomm, &size));
2057: if (size == 1) {
2058: lcolp = colp;
2059: } else {
2060: PetscCall(ISAllGather(colp, &lcolp));
2061: }
2062: PetscCall(ISSetPermutation(lcolp));
2063: /* now we just get the submatrix */
2064: PetscCall(MatGetLocalSize(A, NULL, &clocal_size));
2065: PetscCall(MatCreateSubMatrix_MPIBAIJ_Private(A, crowp, lcolp, clocal_size, MAT_INITIAL_MATRIX, B, PETSC_FALSE));
2066: /* clean up */
2067: if (pcomm != comm) PetscCall(ISDestroy(&crowp));
2068: if (size > 1) PetscCall(ISDestroy(&lcolp));
2069: PetscFunctionReturn(PETSC_SUCCESS);
2070: }
2072: static PetscErrorCode MatGetGhosts_MPIBAIJ(Mat mat, PetscInt *nghosts, const PetscInt *ghosts[])
2073: {
2074: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
2075: Mat_SeqBAIJ *B = (Mat_SeqBAIJ *)baij->B->data;
2077: PetscFunctionBegin;
2078: if (nghosts) *nghosts = B->nbs;
2079: if (ghosts) *ghosts = baij->garray;
2080: PetscFunctionReturn(PETSC_SUCCESS);
2081: }
2083: static PetscErrorCode MatGetSeqNonzeroStructure_MPIBAIJ(Mat A, Mat *newmat)
2084: {
2085: Mat B;
2086: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
2087: Mat_SeqBAIJ *ad = (Mat_SeqBAIJ *)a->A->data, *bd = (Mat_SeqBAIJ *)a->B->data;
2088: Mat_SeqAIJ *b;
2089: PetscMPIInt size, rank, *recvcounts = NULL, *displs = NULL;
2090: PetscInt sendcount, i, *rstarts = A->rmap->range, n, cnt, j, bs = A->rmap->bs;
2091: PetscInt m, *garray = a->garray, *lens, *jsendbuf, *a_jsendbuf, *b_jsendbuf;
2093: PetscFunctionBegin;
2094: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
2095: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
2097: /* Tell every processor the number of nonzeros per row */
2098: PetscCall(PetscMalloc1(A->rmap->N / bs, &lens));
2099: for (i = A->rmap->rstart / bs; i < A->rmap->rend / bs; i++) lens[i] = ad->i[i - A->rmap->rstart / bs + 1] - ad->i[i - A->rmap->rstart / bs] + bd->i[i - A->rmap->rstart / bs + 1] - bd->i[i - A->rmap->rstart / bs];
2100: PetscCall(PetscMalloc1(2 * size, &recvcounts));
2101: displs = recvcounts + size;
2102: for (i = 0; i < size; i++) {
2103: PetscCall(PetscMPIIntCast(A->rmap->range[i + 1] / bs - A->rmap->range[i] / bs, &recvcounts[i]));
2104: PetscCall(PetscMPIIntCast(A->rmap->range[i] / bs, &displs[i]));
2105: }
2106: PetscCallMPI(MPI_Allgatherv(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, lens, recvcounts, displs, MPIU_INT, PetscObjectComm((PetscObject)A)));
2107: /* Create the sequential matrix of the same type as the local block diagonal */
2108: PetscCall(MatCreate(PETSC_COMM_SELF, &B));
2109: PetscCall(MatSetSizes(B, A->rmap->N / bs, A->cmap->N / bs, PETSC_DETERMINE, PETSC_DETERMINE));
2110: PetscCall(MatSetType(B, MATSEQAIJ));
2111: PetscCall(MatSeqAIJSetPreallocation(B, 0, lens));
2112: b = (Mat_SeqAIJ *)B->data;
2114: /* Copy my part of matrix column indices over */
2115: sendcount = ad->nz + bd->nz;
2116: jsendbuf = b->j + b->i[rstarts[rank] / bs];
2117: a_jsendbuf = ad->j;
2118: b_jsendbuf = bd->j;
2119: n = A->rmap->rend / bs - A->rmap->rstart / bs;
2120: cnt = 0;
2121: for (i = 0; i < n; i++) {
2122: /* put in lower diagonal portion */
2123: m = bd->i[i + 1] - bd->i[i];
2124: while (m > 0) {
2125: /* is it above diagonal (in bd (compressed) numbering) */
2126: if (garray[*b_jsendbuf] > A->rmap->rstart / bs + i) break;
2127: jsendbuf[cnt++] = garray[*b_jsendbuf++];
2128: m--;
2129: }
2131: /* put in diagonal portion */
2132: for (j = ad->i[i]; j < ad->i[i + 1]; j++) jsendbuf[cnt++] = A->rmap->rstart / bs + *a_jsendbuf++;
2134: /* put in upper diagonal portion */
2135: while (m-- > 0) jsendbuf[cnt++] = garray[*b_jsendbuf++];
2136: }
2137: PetscCheck(cnt == sendcount, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted PETSc matrix: nz given %" PetscInt_FMT " actual nz %" PetscInt_FMT, sendcount, cnt);
2139: /* Gather all column indices to all processors */
2140: for (i = 0; i < size; i++) {
2141: recvcounts[i] = 0;
2142: for (j = A->rmap->range[i] / bs; j < A->rmap->range[i + 1] / bs; j++) recvcounts[i] += lens[j];
2143: }
2144: displs[0] = 0;
2145: for (i = 1; i < size; i++) displs[i] = displs[i - 1] + recvcounts[i - 1];
2146: PetscCallMPI(MPI_Allgatherv(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, b->j, recvcounts, displs, MPIU_INT, PetscObjectComm((PetscObject)A)));
2147: /* Assemble the matrix into usable form (note numerical values not yet set) */
2148: /* set the b->ilen (length of each row) values */
2149: PetscCall(PetscArraycpy(b->ilen, lens, A->rmap->N / bs));
2150: /* set the b->i indices */
2151: b->i[0] = 0;
2152: for (i = 1; i <= A->rmap->N / bs; i++) b->i[i] = b->i[i - 1] + lens[i - 1];
2153: PetscCall(PetscFree(lens));
2154: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2155: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2156: PetscCall(PetscFree(recvcounts));
2158: PetscCall(MatPropagateSymmetryOptions(A, B));
2159: *newmat = B;
2160: PetscFunctionReturn(PETSC_SUCCESS);
2161: }
2163: static PetscErrorCode MatSOR_MPIBAIJ(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
2164: {
2165: Mat_MPIBAIJ *mat = (Mat_MPIBAIJ *)matin->data;
2166: Vec bb1 = NULL;
2168: PetscFunctionBegin;
2169: if (flag == SOR_APPLY_UPPER) {
2170: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
2171: PetscFunctionReturn(PETSC_SUCCESS);
2172: }
2174: if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS) PetscCall(VecDuplicate(bb, &bb1));
2176: if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
2177: if (flag & SOR_ZERO_INITIAL_GUESS) {
2178: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
2179: its--;
2180: }
2182: while (its--) {
2183: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2184: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2186: /* update rhs: bb1 = bb - B*x */
2187: PetscCall(VecScale(mat->lvec, -1.0));
2188: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
2190: /* local sweep */
2191: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, 1, xx);
2192: }
2193: } else if (flag & SOR_LOCAL_FORWARD_SWEEP) {
2194: if (flag & SOR_ZERO_INITIAL_GUESS) {
2195: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
2196: its--;
2197: }
2198: while (its--) {
2199: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2200: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2202: /* update rhs: bb1 = bb - B*x */
2203: PetscCall(VecScale(mat->lvec, -1.0));
2204: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
2206: /* local sweep */
2207: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_FORWARD_SWEEP, fshift, lits, 1, xx);
2208: }
2209: } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
2210: if (flag & SOR_ZERO_INITIAL_GUESS) {
2211: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
2212: its--;
2213: }
2214: while (its--) {
2215: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2216: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2218: /* update rhs: bb1 = bb - B*x */
2219: PetscCall(VecScale(mat->lvec, -1.0));
2220: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
2222: /* local sweep */
2223: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_BACKWARD_SWEEP, fshift, lits, 1, xx);
2224: }
2225: } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_SUP, "Parallel version of SOR requested not supported");
2227: PetscCall(VecDestroy(&bb1));
2228: PetscFunctionReturn(PETSC_SUCCESS);
2229: }
2231: static PetscErrorCode MatGetColumnReductions_MPIBAIJ(Mat A, PetscInt type, PetscReal *reductions)
2232: {
2233: Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)A->data;
2234: PetscInt m, N, i, *garray = aij->garray;
2235: PetscInt ib, jb, bs = A->rmap->bs;
2236: Mat_SeqBAIJ *a_aij = (Mat_SeqBAIJ *)aij->A->data;
2237: MatScalar *a_val = a_aij->a;
2238: Mat_SeqBAIJ *b_aij = (Mat_SeqBAIJ *)aij->B->data;
2239: MatScalar *b_val = b_aij->a;
2241: PetscFunctionBegin;
2242: PetscCall(MatGetSize(A, &m, &N));
2243: PetscCall(PetscArrayzero(reductions, N));
2244: if (type == NORM_2) {
2245: for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2246: for (jb = 0; jb < bs; jb++) {
2247: for (ib = 0; ib < bs; ib++) {
2248: reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val * *a_val);
2249: a_val++;
2250: }
2251: }
2252: }
2253: for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2254: for (jb = 0; jb < bs; jb++) {
2255: for (ib = 0; ib < bs; ib++) {
2256: reductions[garray[b_aij->j[i]] * bs + jb] += PetscAbsScalar(*b_val * *b_val);
2257: b_val++;
2258: }
2259: }
2260: }
2261: } else if (type == NORM_1) {
2262: for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2263: for (jb = 0; jb < bs; jb++) {
2264: for (ib = 0; ib < bs; ib++) {
2265: reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val);
2266: a_val++;
2267: }
2268: }
2269: }
2270: for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2271: for (jb = 0; jb < bs; jb++) {
2272: for (ib = 0; ib < bs; ib++) {
2273: reductions[garray[b_aij->j[i]] * bs + jb] += PetscAbsScalar(*b_val);
2274: b_val++;
2275: }
2276: }
2277: }
2278: } else if (type == NORM_INFINITY) {
2279: for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2280: for (jb = 0; jb < bs; jb++) {
2281: for (ib = 0; ib < bs; ib++) {
2282: PetscInt col = A->cmap->rstart + a_aij->j[i] * bs + jb;
2283: reductions[col] = PetscMax(PetscAbsScalar(*a_val), reductions[col]);
2284: a_val++;
2285: }
2286: }
2287: }
2288: for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2289: for (jb = 0; jb < bs; jb++) {
2290: for (ib = 0; ib < bs; ib++) {
2291: PetscInt col = garray[b_aij->j[i]] * bs + jb;
2292: reductions[col] = PetscMax(PetscAbsScalar(*b_val), reductions[col]);
2293: b_val++;
2294: }
2295: }
2296: }
2297: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
2298: for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2299: for (jb = 0; jb < bs; jb++) {
2300: for (ib = 0; ib < bs; ib++) {
2301: reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscRealPart(*a_val);
2302: a_val++;
2303: }
2304: }
2305: }
2306: for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2307: for (jb = 0; jb < bs; jb++) {
2308: for (ib = 0; ib < bs; ib++) {
2309: reductions[garray[b_aij->j[i]] * bs + jb] += PetscRealPart(*b_val);
2310: b_val++;
2311: }
2312: }
2313: }
2314: } else {
2315: PetscCheck(type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
2316: for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2317: for (jb = 0; jb < bs; jb++) {
2318: for (ib = 0; ib < bs; ib++) {
2319: reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscImaginaryPart(*a_val);
2320: a_val++;
2321: }
2322: }
2323: }
2324: for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2325: for (jb = 0; jb < bs; jb++) {
2326: for (ib = 0; ib < bs; ib++) {
2327: reductions[garray[b_aij->j[i]] * bs + jb] += PetscImaginaryPart(*b_val);
2328: b_val++;
2329: }
2330: }
2331: }
2332: }
2333: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, reductions, N, MPIU_REAL, type == NORM_INFINITY ? MPIU_MAX : MPIU_SUM, PetscObjectComm((PetscObject)A)));
2334: if (type == NORM_2) {
2335: for (i = 0; i < N; i++) reductions[i] = PetscSqrtReal(reductions[i]);
2336: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2337: for (i = 0; i < N; i++) reductions[i] /= m;
2338: }
2339: PetscFunctionReturn(PETSC_SUCCESS);
2340: }
2342: static PetscErrorCode MatInvertBlockDiagonal_MPIBAIJ(Mat A, const PetscScalar **values)
2343: {
2344: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
2346: PetscFunctionBegin;
2347: PetscCall(MatInvertBlockDiagonal(a->A, values));
2348: A->factorerrortype = a->A->factorerrortype;
2349: A->factorerror_zeropivot_value = a->A->factorerror_zeropivot_value;
2350: A->factorerror_zeropivot_row = a->A->factorerror_zeropivot_row;
2351: PetscFunctionReturn(PETSC_SUCCESS);
2352: }
2354: static PetscErrorCode MatShift_MPIBAIJ(Mat Y, PetscScalar a)
2355: {
2356: Mat_MPIBAIJ *maij = (Mat_MPIBAIJ *)Y->data;
2357: Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)maij->A->data;
2359: PetscFunctionBegin;
2360: if (!Y->preallocated) {
2361: PetscCall(MatMPIBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL, 0, NULL));
2362: } else if (!aij->nz) {
2363: PetscInt nonew = aij->nonew;
2364: PetscCall(MatSeqBAIJSetPreallocation(maij->A, Y->rmap->bs, 1, NULL));
2365: aij->nonew = nonew;
2366: }
2367: PetscCall(MatShift_Basic(Y, a));
2368: PetscFunctionReturn(PETSC_SUCCESS);
2369: }
2371: static PetscErrorCode MatGetDiagonalBlock_MPIBAIJ(Mat A, Mat *a)
2372: {
2373: PetscFunctionBegin;
2374: *a = ((Mat_MPIBAIJ *)A->data)->A;
2375: PetscFunctionReturn(PETSC_SUCCESS);
2376: }
2378: static PetscErrorCode MatEliminateZeros_MPIBAIJ(Mat A, PetscBool keep)
2379: {
2380: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
2382: PetscFunctionBegin;
2383: PetscCall(MatEliminateZeros_SeqBAIJ(a->A, keep)); // possibly keep zero diagonal coefficients
2384: PetscCall(MatEliminateZeros_SeqBAIJ(a->B, PETSC_FALSE)); // never keep zero diagonal coefficients
2385: PetscFunctionReturn(PETSC_SUCCESS);
2386: }
2388: static struct _MatOps MatOps_Values = {MatSetValues_MPIBAIJ,
2389: MatGetRow_MPIBAIJ,
2390: MatRestoreRow_MPIBAIJ,
2391: MatMult_MPIBAIJ,
2392: /* 4*/ MatMultAdd_MPIBAIJ,
2393: MatMultTranspose_MPIBAIJ,
2394: MatMultTransposeAdd_MPIBAIJ,
2395: NULL,
2396: NULL,
2397: NULL,
2398: /*10*/ NULL,
2399: NULL,
2400: NULL,
2401: MatSOR_MPIBAIJ,
2402: MatTranspose_MPIBAIJ,
2403: /*15*/ MatGetInfo_MPIBAIJ,
2404: MatEqual_MPIBAIJ,
2405: MatGetDiagonal_MPIBAIJ,
2406: MatDiagonalScale_MPIBAIJ,
2407: MatNorm_MPIBAIJ,
2408: /*20*/ MatAssemblyBegin_MPIBAIJ,
2409: MatAssemblyEnd_MPIBAIJ,
2410: MatSetOption_MPIBAIJ,
2411: MatZeroEntries_MPIBAIJ,
2412: /*24*/ MatZeroRows_MPIBAIJ,
2413: NULL,
2414: NULL,
2415: NULL,
2416: NULL,
2417: /*29*/ MatSetUp_MPI_Hash,
2418: NULL,
2419: NULL,
2420: MatGetDiagonalBlock_MPIBAIJ,
2421: NULL,
2422: /*34*/ MatDuplicate_MPIBAIJ,
2423: NULL,
2424: NULL,
2425: NULL,
2426: NULL,
2427: /*39*/ MatAXPY_MPIBAIJ,
2428: MatCreateSubMatrices_MPIBAIJ,
2429: MatIncreaseOverlap_MPIBAIJ,
2430: MatGetValues_MPIBAIJ,
2431: MatCopy_MPIBAIJ,
2432: /*44*/ NULL,
2433: MatScale_MPIBAIJ,
2434: MatShift_MPIBAIJ,
2435: NULL,
2436: MatZeroRowsColumns_MPIBAIJ,
2437: /*49*/ NULL,
2438: NULL,
2439: NULL,
2440: NULL,
2441: NULL,
2442: /*54*/ MatFDColoringCreate_MPIXAIJ,
2443: NULL,
2444: MatSetUnfactored_MPIBAIJ,
2445: MatPermute_MPIBAIJ,
2446: MatSetValuesBlocked_MPIBAIJ,
2447: /*59*/ MatCreateSubMatrix_MPIBAIJ,
2448: MatDestroy_MPIBAIJ,
2449: MatView_MPIBAIJ,
2450: NULL,
2451: NULL,
2452: /*64*/ NULL,
2453: NULL,
2454: NULL,
2455: NULL,
2456: MatGetRowMaxAbs_MPIBAIJ,
2457: /*69*/ NULL,
2458: NULL,
2459: NULL,
2460: MatFDColoringApply_BAIJ,
2461: NULL,
2462: /*74*/ NULL,
2463: NULL,
2464: NULL,
2465: NULL,
2466: MatLoad_MPIBAIJ,
2467: /*79*/ NULL,
2468: NULL,
2469: NULL,
2470: NULL,
2471: NULL,
2472: /*84*/ NULL,
2473: NULL,
2474: NULL,
2475: NULL,
2476: NULL,
2477: /*89*/ NULL,
2478: NULL,
2479: NULL,
2480: NULL,
2481: MatConjugate_MPIBAIJ,
2482: /*94*/ NULL,
2483: NULL,
2484: MatRealPart_MPIBAIJ,
2485: MatImaginaryPart_MPIBAIJ,
2486: NULL,
2487: /*99*/ NULL,
2488: NULL,
2489: NULL,
2490: NULL,
2491: NULL,
2492: /*104*/ MatGetSeqNonzeroStructure_MPIBAIJ,
2493: NULL,
2494: MatGetGhosts_MPIBAIJ,
2495: NULL,
2496: NULL,
2497: /*109*/ NULL,
2498: NULL,
2499: NULL,
2500: NULL,
2501: MatGetMultiProcBlock_MPIBAIJ,
2502: /*114*/ NULL,
2503: MatGetColumnReductions_MPIBAIJ,
2504: MatInvertBlockDiagonal_MPIBAIJ,
2505: NULL,
2506: NULL,
2507: /*119*/ NULL,
2508: NULL,
2509: NULL,
2510: NULL,
2511: NULL,
2512: /*124*/ NULL,
2513: MatSetBlockSizes_Default,
2514: NULL,
2515: MatFDColoringSetUp_MPIXAIJ,
2516: NULL,
2517: /*129*/ MatCreateMPIMatConcatenateSeqMat_MPIBAIJ,
2518: NULL,
2519: NULL,
2520: NULL,
2521: NULL,
2522: /*134*/ NULL,
2523: MatEliminateZeros_MPIBAIJ,
2524: MatGetRowSumAbs_MPIBAIJ,
2525: NULL,
2526: NULL,
2527: /*139*/ NULL,
2528: MatCopyHashToXAIJ_MPI_Hash,
2529: NULL,
2530: NULL,
2531: MatADot_Default,
2532: /*144*/ MatANorm_Default,
2533: NULL,
2534: NULL,
2535: NULL};
2537: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPISBAIJ(Mat, MatType, MatReuse, Mat *);
2538: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat, MatType, MatReuse, Mat *);
2540: static PetscErrorCode MatMPIBAIJSetPreallocationCSR_MPIBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[])
2541: {
2542: PetscInt m, rstart, cstart, cend;
2543: PetscInt i, j, dlen, olen, nz, nz_max = 0, *d_nnz = NULL, *o_nnz = NULL;
2544: const PetscInt *JJ = NULL;
2545: PetscScalar *values = NULL;
2546: PetscBool roworiented = ((Mat_MPIBAIJ *)B->data)->roworiented;
2547: PetscBool nooffprocentries;
2549: PetscFunctionBegin;
2550: PetscCall(PetscLayoutSetBlockSize(B->rmap, bs));
2551: PetscCall(PetscLayoutSetBlockSize(B->cmap, bs));
2552: PetscCall(PetscLayoutSetUp(B->rmap));
2553: PetscCall(PetscLayoutSetUp(B->cmap));
2554: PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
2555: m = B->rmap->n / bs;
2556: rstart = B->rmap->rstart / bs;
2557: cstart = B->cmap->rstart / bs;
2558: cend = B->cmap->rend / bs;
2560: PetscCheck(!ii[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]);
2561: PetscCall(PetscMalloc2(m, &d_nnz, m, &o_nnz));
2562: for (i = 0; i < m; i++) {
2563: nz = ii[i + 1] - ii[i];
2564: PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz);
2565: nz_max = PetscMax(nz_max, nz);
2566: dlen = 0;
2567: olen = 0;
2568: JJ = jj + ii[i];
2569: for (j = 0; j < nz; j++) {
2570: if (*JJ < cstart || *JJ >= cend) olen++;
2571: else dlen++;
2572: JJ++;
2573: }
2574: d_nnz[i] = dlen;
2575: o_nnz[i] = olen;
2576: }
2577: PetscCall(MatMPIBAIJSetPreallocation(B, bs, 0, d_nnz, 0, o_nnz));
2578: PetscCall(PetscFree2(d_nnz, o_nnz));
2580: values = (PetscScalar *)V;
2581: if (!values) PetscCall(PetscCalloc1(bs * bs * nz_max, &values));
2582: for (i = 0; i < m; i++) {
2583: PetscInt row = i + rstart;
2584: PetscInt ncols = ii[i + 1] - ii[i];
2585: const PetscInt *icols = jj + ii[i];
2586: if (bs == 1 || !roworiented) { /* block ordering matches the non-nested layout of MatSetValues so we can insert entire rows */
2587: const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0);
2588: PetscCall(MatSetValuesBlocked_MPIBAIJ(B, 1, &row, ncols, icols, svals, INSERT_VALUES));
2589: } else { /* block ordering does not match so we can only insert one block at a time. */
2590: PetscInt j;
2591: for (j = 0; j < ncols; j++) {
2592: const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0);
2593: PetscCall(MatSetValuesBlocked_MPIBAIJ(B, 1, &row, 1, &icols[j], svals, INSERT_VALUES));
2594: }
2595: }
2596: }
2598: if (!V) PetscCall(PetscFree(values));
2599: nooffprocentries = B->nooffprocentries;
2600: B->nooffprocentries = PETSC_TRUE;
2601: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2602: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2603: B->nooffprocentries = nooffprocentries;
2605: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
2606: PetscFunctionReturn(PETSC_SUCCESS);
2607: }
2609: /*@C
2610: MatMPIBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATBAIJ` format using the given nonzero structure and (optional) numerical values
2612: Collective
2614: Input Parameters:
2615: + B - the matrix
2616: . bs - the block size
2617: . i - the indices into `j` for the start of each local row (starts with zero)
2618: . j - the column indices for each local row (starts with zero) these must be sorted for each row
2619: - v - optional values in the matrix, use `NULL` if not provided
2621: Level: advanced
2623: Notes:
2624: The `i`, `j`, and `v` arrays ARE copied by this routine into the internal format used by PETSc;
2625: thus you CANNOT change the matrix entries by changing the values of `v` after you have
2626: called this routine.
2628: The order of the entries in values is specified by the `MatOption` `MAT_ROW_ORIENTED`. For example, C programs
2629: may want to use the default `MAT_ROW_ORIENTED` with value `PETSC_TRUE` and use an array v[nnz][bs][bs] where the second index is
2630: over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
2631: `MAT_ROW_ORIENTED` with value `PETSC_FALSE` and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2632: block column and the second index is over columns within a block.
2634: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries and usually the numerical values as well
2636: .seealso: `Mat`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatMPIBAIJSetPreallocation()`, `MatCreateAIJ()`, `MATMPIAIJ`, `MatCreateMPIBAIJWithArrays()`, `MATMPIBAIJ`
2637: @*/
2638: PetscErrorCode MatMPIBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[])
2639: {
2640: PetscFunctionBegin;
2644: PetscTryMethod(B, "MatMPIBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v));
2645: PetscFunctionReturn(PETSC_SUCCESS);
2646: }
2648: PetscErrorCode MatMPIBAIJSetPreallocation_MPIBAIJ(Mat B, PetscInt bs, PetscInt d_nz, const PetscInt *d_nnz, PetscInt o_nz, const PetscInt *o_nnz)
2649: {
2650: Mat_MPIBAIJ *b = (Mat_MPIBAIJ *)B->data;
2651: PetscInt i;
2652: PetscMPIInt size;
2654: PetscFunctionBegin;
2655: if (B->hash_active) {
2656: B->ops[0] = b->cops;
2657: B->hash_active = PETSC_FALSE;
2658: }
2659: if (!B->preallocated) PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), bs, &B->bstash));
2660: PetscCall(MatSetBlockSize(B, bs));
2661: PetscCall(PetscLayoutSetUp(B->rmap));
2662: PetscCall(PetscLayoutSetUp(B->cmap));
2663: PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
2665: if (d_nnz) {
2666: for (i = 0; i < B->rmap->n / bs; i++) PetscCheck(d_nnz[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "d_nnz cannot be less than -1: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, d_nnz[i]);
2667: }
2668: if (o_nnz) {
2669: for (i = 0; i < B->rmap->n / bs; i++) PetscCheck(o_nnz[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "o_nnz cannot be less than -1: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, o_nnz[i]);
2670: }
2672: b->bs2 = bs * bs;
2673: b->mbs = B->rmap->n / bs;
2674: b->nbs = B->cmap->n / bs;
2675: b->Mbs = B->rmap->N / bs;
2676: b->Nbs = B->cmap->N / bs;
2678: for (i = 0; i <= b->size; i++) b->rangebs[i] = B->rmap->range[i] / bs;
2679: b->rstartbs = B->rmap->rstart / bs;
2680: b->rendbs = B->rmap->rend / bs;
2681: b->cstartbs = B->cmap->rstart / bs;
2682: b->cendbs = B->cmap->rend / bs;
2684: #if PetscDefined(USE_CTABLE)
2685: PetscCall(PetscHMapIDestroy(&b->colmap));
2686: #else
2687: PetscCall(PetscFree(b->colmap));
2688: #endif
2689: PetscCall(PetscFree(b->garray));
2690: PetscCall(VecDestroy(&b->lvec));
2691: PetscCall(VecScatterDestroy(&b->Mvctx));
2693: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
2695: MatSeqXAIJGetOptions_Private(b->B);
2696: PetscCall(MatDestroy(&b->B));
2697: PetscCall(MatCreate(PETSC_COMM_SELF, &b->B));
2698: PetscCall(MatSetSizes(b->B, B->rmap->n, size > 1 ? B->cmap->N : 0, B->rmap->n, size > 1 ? B->cmap->N : 0));
2699: PetscCall(MatSetType(b->B, MATSEQBAIJ));
2700: MatSeqXAIJRestoreOptions_Private(b->B);
2702: MatSeqXAIJGetOptions_Private(b->A);
2703: PetscCall(MatDestroy(&b->A));
2704: PetscCall(MatCreate(PETSC_COMM_SELF, &b->A));
2705: PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n));
2706: PetscCall(MatSetType(b->A, MATSEQBAIJ));
2707: MatSeqXAIJRestoreOptions_Private(b->A);
2709: PetscCall(MatSeqBAIJSetPreallocation(b->A, bs, d_nz, d_nnz));
2710: PetscCall(MatSeqBAIJSetPreallocation(b->B, bs, o_nz, o_nnz));
2711: B->preallocated = PETSC_TRUE;
2712: B->was_assembled = PETSC_FALSE;
2713: B->assembled = PETSC_FALSE;
2714: PetscFunctionReturn(PETSC_SUCCESS);
2715: }
2717: extern PetscErrorCode MatDiagonalScaleLocal_MPIBAIJ(Mat, Vec);
2718: extern PetscErrorCode MatSetHashTableFactor_MPIBAIJ(Mat, PetscReal);
2720: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIAdj(Mat B, MatType newtype, MatReuse reuse, Mat *adj)
2721: {
2722: Mat_MPIBAIJ *b = (Mat_MPIBAIJ *)B->data;
2723: Mat_SeqBAIJ *d = (Mat_SeqBAIJ *)b->A->data, *o = (Mat_SeqBAIJ *)b->B->data;
2724: PetscInt M = B->rmap->n / B->rmap->bs, i, *ii, *jj, cnt, j, k, rstart = B->rmap->rstart / B->rmap->bs;
2725: const PetscInt *id = d->i, *jd = d->j, *io = o->i, *jo = o->j, *garray = b->garray;
2727: PetscFunctionBegin;
2728: PetscCall(PetscMalloc1(M + 1, &ii));
2729: ii[0] = 0;
2730: for (i = 0; i < M; i++) {
2731: PetscCheck((id[i + 1] - id[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Indices wrong %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, i, id[i], id[i + 1]);
2732: PetscCheck((io[i + 1] - io[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Indices wrong %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, i, io[i], io[i + 1]);
2733: ii[i + 1] = ii[i] + id[i + 1] - id[i] + io[i + 1] - io[i];
2734: /* remove one from count of matrix has diagonal */
2735: for (j = id[i]; j < id[i + 1]; j++) {
2736: if (jd[j] == i) {
2737: ii[i + 1]--;
2738: break;
2739: }
2740: }
2741: }
2742: PetscCall(PetscMalloc1(ii[M], &jj));
2743: cnt = 0;
2744: for (i = 0; i < M; i++) {
2745: for (j = io[i]; j < io[i + 1]; j++) {
2746: if (garray[jo[j]] > rstart) break;
2747: jj[cnt++] = garray[jo[j]];
2748: }
2749: for (k = id[i]; k < id[i + 1]; k++) {
2750: if (jd[k] != i) jj[cnt++] = rstart + jd[k];
2751: }
2752: for (; j < io[i + 1]; j++) jj[cnt++] = garray[jo[j]];
2753: }
2754: PetscCall(MatCreateMPIAdj(PetscObjectComm((PetscObject)B), M, B->cmap->N / B->rmap->bs, ii, jj, NULL, adj));
2755: PetscFunctionReturn(PETSC_SUCCESS);
2756: }
2758: #include <../src/mat/impls/aij/mpi/mpiaij.h>
2760: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *);
2762: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2763: {
2764: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
2765: Mat_MPIAIJ *b;
2766: Mat B;
2768: PetscFunctionBegin;
2769: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled");
2771: if (reuse == MAT_REUSE_MATRIX) {
2772: B = *newmat;
2773: } else {
2774: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
2775: PetscCall(MatSetType(B, MATMPIAIJ));
2776: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
2777: PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs));
2778: PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));
2779: PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));
2780: }
2781: b = (Mat_MPIAIJ *)B->data;
2783: if (reuse == MAT_REUSE_MATRIX) {
2784: PetscCall(MatConvert_SeqBAIJ_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A));
2785: PetscCall(MatConvert_SeqBAIJ_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B));
2786: } else {
2787: PetscInt *garray = a->garray;
2788: Mat_SeqAIJ *bB;
2789: PetscInt bs, nnz;
2790: PetscCall(MatDestroy(&b->A));
2791: PetscCall(MatDestroy(&b->B));
2792: /* just clear out the data structure */
2793: PetscCall(MatDisAssemble_MPIAIJ(B, PETSC_FALSE));
2794: PetscCall(MatConvert_SeqBAIJ_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A));
2795: PetscCall(MatConvert_SeqBAIJ_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B));
2797: /* Global numbering for b->B columns */
2798: bB = (Mat_SeqAIJ *)b->B->data;
2799: bs = A->rmap->bs;
2800: nnz = bB->i[A->rmap->n];
2801: for (PetscInt k = 0; k < nnz; k++) {
2802: PetscInt bj = bB->j[k] / bs;
2803: PetscInt br = bB->j[k] % bs;
2804: bB->j[k] = garray[bj] * bs + br;
2805: }
2806: }
2807: PetscCall(MatSetOption(B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2808: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2809: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2810: PetscCall(MatSetOption(B, MAT_NO_OFF_PROC_ENTRIES, PETSC_FALSE));
2812: if (reuse == MAT_INPLACE_MATRIX) {
2813: PetscCall(MatHeaderReplace(A, &B));
2814: } else {
2815: *newmat = B;
2816: }
2817: PetscFunctionReturn(PETSC_SUCCESS);
2818: }
2820: /*MC
2821: MATMPIBAIJ - MATMPIBAIJ = "mpibaij" - A matrix type to be used for distributed block sparse matrices.
2823: Options Database Keys:
2824: + -mat_type mpibaij - sets the matrix type to `MATMPIBAIJ` during a call to `MatSetFromOptions()`
2825: . -mat_block_size bs - set the blocksize used to store the matrix
2826: . -mat_baij_mult_version version - indicate the version of the matrix-vector product to use (0 often indicates using BLAS)
2827: - -mat_use_hash_table fact - set hash table factor
2829: Level: beginner
2831: Note:
2832: `MatSetOption(A, MAT_STRUCTURE_ONLY, PETSC_TRUE)` may be called for this matrix type. In this no
2833: space is allocated for the nonzero entries and any entries passed with `MatSetValues()` are ignored
2835: .seealso: `Mat`, `MATBAIJ`, `MATSEQBAIJ`, `MatCreateBAIJ`
2836: M*/
2838: typedef struct {
2839: MPIAIJ_MPIDense scatter;
2840: Mat workC;
2841: } MPIBAIJ_MPIDense;
2843: static PetscErrorCode MatMPIBAIJ_MPIDenseDestroy(PetscCtxRt ctx)
2844: {
2845: MPIBAIJ_MPIDense *data = *(MPIBAIJ_MPIDense **)ctx;
2847: PetscFunctionBegin;
2848: PetscCall(MatDestroy(&data->workC));
2849: PetscCall(MatMPIDenseScatterDestroy_Private(&data->scatter));
2850: PetscCall(PetscFree(data));
2851: PetscFunctionReturn(PETSC_SUCCESS);
2852: }
2854: static PetscErrorCode MatMPIDenseScatter_MPIBAIJ(Mat A, Mat B, Mat workB, Mat C)
2855: {
2856: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)A->data;
2857: MPIBAIJ_MPIDense *data = (MPIBAIJ_MPIDense *)C->product->data;
2858: PetscInt bs;
2860: PetscFunctionBegin;
2861: PetscCall(MatGetBlockSize(A, &bs));
2862: PetscCall(MatMPIDenseScatter_Private(baij->Mvctx, baij->B->cmap->n, bs, workB, &data->scatter, B, C));
2863: PetscFunctionReturn(PETSC_SUCCESS);
2864: }
2866: static PetscErrorCode MatMatMultNumeric_MPIBAIJ_MPIDense(Mat A, Mat B, Mat C)
2867: {
2868: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)A->data;
2869: Mat_MPIDense *bdense = (Mat_MPIDense *)B->data;
2870: Mat_MPIDense *cdense = (Mat_MPIDense *)C->data;
2871: Mat workB;
2872: MPIBAIJ_MPIDense *data;
2874: PetscFunctionBegin;
2875: MatCheckProduct(C, 3);
2876: PetscCheck(C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
2877: data = (MPIBAIJ_MPIDense *)C->product->data;
2878: if (!cdense->A->product) {
2879: PetscCall(MatProductCreateWithMat(baij->A, bdense->A, NULL, cdense->A));
2880: PetscCall(MatProductSetType(cdense->A, MATPRODUCT_AB));
2881: PetscCall(MatProductSetFromOptions(cdense->A));
2882: PetscCall(MatProductSymbolic(cdense->A));
2883: } else PetscCall(MatProductReplaceMats(baij->A, bdense->A, NULL, cdense->A));
2884: PetscCall(MatProductNumeric(cdense->A));
2886: if (data->scatter.workB->cmap->n == B->cmap->N) {
2887: workB = data->scatter.workB;
2888: PetscCall(MatMPIDenseScatter_MPIBAIJ(A, B, workB, C));
2889: if (data->workC) {
2890: PetscCall(MatProductReplaceMats(baij->B, workB, NULL, data->workC));
2891: PetscCall(MatProductNumeric(data->workC));
2892: PetscCall(MatAXPY(cdense->A, 1.0, data->workC, SAME_NONZERO_PATTERN));
2893: }
2894: } else {
2895: Mat Bb, Cb, workC;
2896: Mat_MPIDense *cbdense;
2897: PetscInt BN = B->cmap->N, n = data->scatter.workB->cmap->n, cols;
2899: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Column batch size must be positive");
2900: for (PetscInt i = 0; i < BN; i += n) {
2901: cols = PetscMin(n, BN - i);
2902: workB = data->scatter.workB;
2903: workC = data->workC;
2904: if (cols != n) {
2905: PetscCall(MatDenseGetSubMatrix(data->scatter.workB, PETSC_DECIDE, PETSC_DECIDE, 0, cols, &workB));
2906: if (workC) PetscCall(MatDenseGetSubMatrix(data->workC, PETSC_DECIDE, PETSC_DECIDE, 0, cols, &workC));
2907: }
2908: PetscCall(MatDenseGetSubMatrix(B, PETSC_DECIDE, PETSC_DECIDE, i, i + cols, &Bb));
2909: PetscCall(MatDenseGetSubMatrix(C, PETSC_DECIDE, PETSC_DECIDE, i, i + cols, &Cb));
2910: PetscCall(MatMPIDenseScatter_MPIBAIJ(A, Bb, workB, C));
2911: if (workC) {
2912: cbdense = (Mat_MPIDense *)Cb->data;
2913: PetscCall(MatProductReplaceMats(baij->B, workB, NULL, workC));
2914: PetscCall(MatProductNumeric(workC));
2915: PetscCall(MatAXPY(cbdense->A, 1.0, workC, SAME_NONZERO_PATTERN));
2916: }
2917: if (cols != n) {
2918: if (workC) PetscCall(MatDenseRestoreSubMatrix(data->workC, &workC));
2919: PetscCall(MatDenseRestoreSubMatrix(data->scatter.workB, &workB));
2920: }
2921: PetscCall(MatDenseRestoreSubMatrix(B, &Bb));
2922: PetscCall(MatDenseRestoreSubMatrix(C, &Cb));
2923: }
2924: }
2925: PetscFunctionReturn(PETSC_SUCCESS);
2926: }
2928: static PetscErrorCode MatMatMultSymbolic_MPIBAIJ_MPIDense(Mat A, Mat B, PetscReal fill, Mat C)
2929: {
2930: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)A->data;
2931: MPIBAIJ_MPIDense *data;
2932: VecScatter ctx = baij->Mvctx;
2933: PetscInt nz = baij->B->cmap->n, m, M, n, N, bs;
2934: PetscInt Am = A->rmap->n, BN = B->cmap->N, Bbn, numBb;
2935: Mat workB1, workC1;
2936: PetscBool cisdense;
2938: PetscFunctionBegin;
2939: MatCheckProduct(C, 4);
2940: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
2941: PetscCall(PetscObjectBaseTypeCompare((PetscObject)C, MATMPIDENSE, &cisdense));
2942: if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
2943: PetscCall(MatGetLocalSize(C, &m, &n));
2944: PetscCall(MatGetSize(C, &M, &N));
2945: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) PetscCall(MatSetSizes(C, Am, B->cmap->n, A->rmap->N, BN));
2946: PetscCall(MatSetBlockSizesFromMats(C, A, B));
2947: PetscCall(MatSetUp(C));
2948: PetscCall(MatGetBlockSize(A, &bs));
2949: PetscCall(PetscNew(&data));
2950: PetscCall(MatMPIDenseScatterSetUp_Private(ctx, nz, bs, Am, B, C, &data->scatter, &Bbn, &numBb));
2952: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2953: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
2954: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
2955: PetscCall(MatProductClear(baij->A));
2956: PetscCall(MatProductClear(((Mat_MPIDense *)B->data)->A));
2957: PetscCall(MatProductClear(((Mat_MPIDense *)C->data)->A));
2958: PetscCall(MatProductCreateWithMat(baij->A, ((Mat_MPIDense *)B->data)->A, NULL, ((Mat_MPIDense *)C->data)->A));
2959: PetscCall(MatProductSetType(((Mat_MPIDense *)C->data)->A, MATPRODUCT_AB));
2960: PetscCall(MatProductSetFromOptions(((Mat_MPIDense *)C->data)->A));
2961: PetscCall(MatProductSymbolic(((Mat_MPIDense *)C->data)->A));
2963: if (nz) {
2964: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Am, Bbn ? Bbn : BN, NULL, &data->workC));
2965: PetscCall(MatProductCreateWithMat(baij->B, data->scatter.workB, NULL, data->workC));
2966: PetscCall(MatProductSetType(data->workC, MATPRODUCT_AB));
2967: PetscCall(MatProductSetFromOptions(data->workC));
2968: PetscCall(MatProductSymbolic(data->workC));
2969: if (numBb && BN % Bbn) {
2970: PetscCall(MatDenseGetSubMatrix(data->scatter.workB, PETSC_DECIDE, PETSC_DECIDE, 0, BN % Bbn, &workB1));
2971: PetscCall(MatDenseGetSubMatrix(data->workC, PETSC_DECIDE, PETSC_DECIDE, 0, BN % Bbn, &workC1));
2972: PetscCall(MatProductCreateWithMat(baij->B, workB1, NULL, workC1));
2973: PetscCall(MatProductSetType(workC1, MATPRODUCT_AB));
2974: PetscCall(MatProductSetFromOptions(workC1));
2975: PetscCall(MatProductSymbolic(workC1));
2976: PetscCall(MatDenseRestoreSubMatrix(data->workC, &workC1));
2977: PetscCall(MatDenseRestoreSubMatrix(data->scatter.workB, &workB1));
2978: }
2979: }
2981: C->product->data = data;
2982: C->product->destroy = MatMPIBAIJ_MPIDenseDestroy;
2983: C->ops->matmultnumeric = MatMatMultNumeric_MPIBAIJ_MPIDense;
2984: PetscFunctionReturn(PETSC_SUCCESS);
2985: }
2987: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_MPIBAIJ_MPIDense(Mat C)
2988: {
2989: Mat_Product *product = C->product;
2991: PetscFunctionBegin;
2992: MatCheckProduct(C, 1);
2993: if (product->type == MATPRODUCT_AB) {
2994: C->ops->matmultsymbolic = MatMatMultSymbolic_MPIBAIJ_MPIDense;
2995: C->ops->productsymbolic = MatProductSymbolic_AB;
2996: }
2997: PetscFunctionReturn(PETSC_SUCCESS);
2998: }
3000: static PetscErrorCode MatGetMultPetscSF_MPIBAIJ(Mat A, PetscSF *sf)
3001: {
3002: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
3004: PetscFunctionBegin;
3005: *sf = a->Mvctx;
3006: PetscFunctionReturn(PETSC_SUCCESS);
3007: }
3009: PETSC_EXTERN PetscErrorCode MatCreate_MPIBAIJ(Mat B)
3010: {
3011: Mat_MPIBAIJ *b;
3012: PetscBool flg = PETSC_FALSE;
3014: PetscFunctionBegin;
3015: PetscCall(PetscNew(&b));
3016: B->data = (void *)b;
3017: B->ops[0] = MatOps_Values;
3018: B->assembled = PETSC_FALSE;
3020: B->insertmode = NOT_SET_VALUES;
3021: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank));
3022: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &b->size));
3024: /* build local table of row and column ownerships */
3025: PetscCall(PetscMalloc1(b->size + 1, &b->rangebs));
3027: /* build cache for off array entries formed */
3028: PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash));
3030: b->donotstash = PETSC_FALSE;
3031: b->colmap = NULL;
3032: b->garray = NULL;
3033: b->roworiented = PETSC_TRUE;
3035: /* stuff used in block assembly */
3036: b->barray = NULL;
3038: /* stuff used for matrix vector multiply */
3039: b->lvec = NULL;
3040: b->Mvctx = NULL;
3042: /* stuff for MatGetRow() */
3043: b->rowindices = NULL;
3044: b->rowvalues = NULL;
3045: b->getrowactive = PETSC_FALSE;
3047: /* hash table stuff */
3048: b->ht = NULL;
3049: b->hd = NULL;
3050: b->ht_size = 0;
3051: b->ht_flag = PETSC_FALSE;
3052: b->ht_fact = 0;
3053: b->ht_total_ct = 0;
3054: b->ht_insert_ct = 0;
3056: /* stuff for MatCreateSubMatrices_MPIBAIJ_local() */
3057: b->ijonly = PETSC_FALSE;
3059: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_mpiadj_C", MatConvert_MPIBAIJ_MPIAdj));
3060: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_mpiaij_C", MatConvert_MPIBAIJ_MPIAIJ));
3061: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_mpisbaij_C", MatConvert_MPIBAIJ_MPISBAIJ));
3062: #if PetscDefined(HAVE_HYPRE)
3063: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_hypre_C", MatConvert_AIJ_HYPRE));
3064: #endif
3065: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPIBAIJ));
3066: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPIBAIJ));
3067: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIBAIJSetPreallocation_C", MatMPIBAIJSetPreallocation_MPIBAIJ));
3068: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIBAIJSetPreallocationCSR_C", MatMPIBAIJSetPreallocationCSR_MPIBAIJ));
3069: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDiagonalScaleLocal_C", MatDiagonalScaleLocal_MPIBAIJ));
3070: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetHashTableFactor_C", MatSetHashTableFactor_MPIBAIJ));
3071: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_is_C", MatConvert_XAIJ_IS));
3072: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatGetMultPetscSF_C", MatGetMultPetscSF_MPIBAIJ));
3073: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpibaij_mpidense_C", MatProductSetFromOptions_MPIBAIJ_MPIDense));
3074: #if PetscDefined(HAVE_LIBXSMM)
3075: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_mpibaijlibxsmm_C", MatConvert_MPIBAIJ_MPIBAIJLIBXSMM));
3076: #endif
3077: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPIBAIJ));
3079: PetscOptionsBegin(PetscObjectComm((PetscObject)B), NULL, "Options for loading MPIBAIJ matrix 1", "Mat");
3080: PetscCall(PetscOptionsName("-mat_use_hash_table", "Use hash table to save time in constructing matrix", "MatSetOption", &flg));
3081: if (flg) {
3082: PetscReal fact = 1.39;
3083: PetscCall(MatSetOption(B, MAT_USE_HASH_TABLE, PETSC_TRUE));
3084: PetscCall(PetscOptionsReal("-mat_use_hash_table", "Use hash table factor", "MatMPIBAIJSetHashTableFactor", fact, &fact, NULL));
3085: if (fact <= 1.0) fact = 1.39;
3086: PetscCall(MatMPIBAIJSetHashTableFactor(B, fact));
3087: PetscCall(PetscInfo(B, "Hash table Factor used %5.2g\n", (double)fact));
3088: }
3089: PetscOptionsEnd();
3090: PetscFunctionReturn(PETSC_SUCCESS);
3091: }
3093: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
3094: /*MC
3095: MATBAIJ - MATBAIJ = "baij" - A matrix type to be used for block sparse matrices.
3097: This matrix type is identical to `MATSEQBAIJ` when constructed with a single process communicator,
3098: and `MATMPIBAIJ` otherwise.
3100: Options Database Keys:
3101: . -mat_type baij - sets the matrix type to `MATBAIJ` during a call to `MatSetFromOptions()`
3103: Level: beginner
3105: .seealso: `Mat`, `MatCreateBAIJ()`, `MATSEQBAIJ`, `MATMPIBAIJ`, `MatMPIBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocationCSR()`
3106: M*/
3108: /*@
3109: MatMPIBAIJSetPreallocation - Allocates memory for a sparse parallel matrix in `MATMPIBAIJ` format
3110: (block compressed row).
3112: Collective
3114: Input Parameters:
3115: + B - the matrix
3116: . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
3117: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
3118: . d_nz - number of block nonzeros per block row in diagonal portion of local
3119: submatrix (same for all local rows)
3120: . d_nnz - array containing the number of block nonzeros in the various block rows
3121: of the in diagonal portion of the local (possibly different for each block
3122: row) or `NULL`. If you plan to factor the matrix you must leave room for the diagonal entry and
3123: set it even if it is zero.
3124: . o_nz - number of block nonzeros per block row in the off-diagonal portion of local
3125: submatrix (same for all local rows).
3126: - o_nnz - array containing the number of nonzeros in the various block rows of the
3127: off-diagonal portion of the local submatrix (possibly different for
3128: each block row) or `NULL`.
3130: If the *_nnz parameter is given then the *_nz parameter is ignored
3132: Options Database Keys:
3133: + -mat_block_size - size of the blocks to use
3134: - -mat_use_hash_table fact - set hash table factor
3136: Level: intermediate
3138: Notes:
3139: For good matrix assembly performance
3140: the user should preallocate the matrix storage by setting the parameters
3141: `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`). By setting these parameters accurately,
3142: performance can be increased by more than a factor of 50.
3144: If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one processor
3145: than it must be used on all processors that share the object for that argument.
3147: Storage Information:
3148: For a square global matrix we define each processor's diagonal portion
3149: to be its local rows and the corresponding columns (a square submatrix);
3150: each processor's off-diagonal portion encompasses the remainder of the
3151: local matrix (a rectangular submatrix).
3153: The user can specify preallocated storage for the diagonal part of
3154: the local submatrix with either `d_nz` or `d_nnz` (not both). Set
3155: `d_nz` = `PETSC_DEFAULT` and `d_nnz` = `NULL` for PETSc to control dynamic
3156: memory allocation. Likewise, specify preallocated storage for the
3157: off-diagonal part of the local submatrix with `o_nz` or `o_nnz` (not both).
3159: Consider a processor that owns rows 3, 4 and 5 of a parallel matrix. In
3160: the figure below we depict these three local rows and all columns (0-11).
3162: .vb
3163: 0 1 2 3 4 5 6 7 8 9 10 11
3164: --------------------------
3165: row 3 |o o o d d d o o o o o o
3166: row 4 |o o o d d d o o o o o o
3167: row 5 |o o o d d d o o o o o o
3168: --------------------------
3169: .ve
3171: Thus, any entries in the d locations are stored in the d (diagonal)
3172: submatrix, and any entries in the o locations are stored in the
3173: o (off-diagonal) submatrix. Note that the d and the o submatrices are
3174: stored simply in the `MATSEQBAIJ` format for compressed row storage.
3176: Now `d_nz` should indicate the number of block nonzeros per row in the d matrix,
3177: and `o_nz` should indicate the number of block nonzeros per row in the o matrix.
3178: In general, for PDE problems in which most nonzeros are near the diagonal,
3179: one expects `d_nz` >> `o_nz`.
3181: You can call `MatGetInfo()` to get information on how effective the preallocation was;
3182: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3183: You can also run with the option `-info` and look for messages with the string
3184: malloc in them to see if additional memory allocation was needed.
3186: .seealso: `Mat`, `MATMPIBAIJ`, `MatCreate()`, `MatCreateSeqBAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`, `MatMPIBAIJSetPreallocationCSR()`, `PetscSplitOwnership()`
3187: @*/
3188: PetscErrorCode MatMPIBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
3189: {
3190: PetscFunctionBegin;
3194: PetscTryMethod(B, "MatMPIBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, bs, d_nz, d_nnz, o_nz, o_nnz));
3195: PetscFunctionReturn(PETSC_SUCCESS);
3196: }
3198: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
3199: /*@
3200: MatCreateBAIJ - Creates a sparse parallel matrix in `MATBAIJ` format
3201: (block compressed row).
3203: Collective
3205: Input Parameters:
3206: + comm - MPI communicator
3207: . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
3208: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
3209: . m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given)
3210: This value should be the same as the local size used in creating the
3211: y vector for the matrix-vector product y = Ax.
3212: . n - number of local columns (or `PETSC_DECIDE` to have calculated if N is given)
3213: This value should be the same as the local size used in creating the
3214: x vector for the matrix-vector product y = Ax.
3215: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
3216: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
3217: . d_nz - number of nonzero blocks per block row in diagonal portion of local
3218: submatrix (same for all local rows)
3219: . d_nnz - array containing the number of nonzero blocks in the various block rows
3220: of the in diagonal portion of the local (possibly different for each block
3221: row) or NULL. If you plan to factor the matrix you must leave room for the diagonal entry
3222: and set it even if it is zero.
3223: . o_nz - number of nonzero blocks per block row in the off-diagonal portion of local
3224: submatrix (same for all local rows).
3225: - o_nnz - array containing the number of nonzero blocks in the various block rows of the
3226: off-diagonal portion of the local submatrix (possibly different for
3227: each block row) or NULL.
3229: Output Parameter:
3230: . A - the matrix
3232: Options Database Keys:
3233: + -mat_block_size - size of the blocks to use
3234: - -mat_use_hash_table fact - set hash table factor
3236: Level: intermediate
3238: Notes:
3239: It is recommended that one use `MatCreateFromOptions()` or the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
3240: MatXXXXSetPreallocation() paradigm instead of this routine directly.
3241: [MatXXXXSetPreallocation() is, for example, `MatSeqBAIJSetPreallocation()`]
3243: For good matrix assembly performance
3244: the user should preallocate the matrix storage by setting the parameters
3245: `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`). By setting these parameters accurately,
3246: performance can be increased by more than a factor of 50.
3248: If the *_nnz parameter is given then the *_nz parameter is ignored
3250: A nonzero block is any block that as 1 or more nonzeros in it
3252: The user MUST specify either the local or global matrix dimensions
3253: (possibly both).
3255: If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one processor
3256: than it must be used on all processors that share the object for that argument.
3258: If `m` and `n` are not `PETSC_DECIDE`, then the values determine the `PetscLayout` of the matrix and the ranges returned by
3259: `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `MatGetOwnershipRangeColumn()`, and `MatGetOwnershipRangesColumn()`.
3261: Storage Information:
3262: For a square global matrix we define each processor's diagonal portion
3263: to be its local rows and the corresponding columns (a square submatrix);
3264: each processor's off-diagonal portion encompasses the remainder of the
3265: local matrix (a rectangular submatrix).
3267: The user can specify preallocated storage for the diagonal part of
3268: the local submatrix with either d_nz or d_nnz (not both). Set
3269: `d_nz` = `PETSC_DEFAULT` and `d_nnz` = `NULL` for PETSc to control dynamic
3270: memory allocation. Likewise, specify preallocated storage for the
3271: off-diagonal part of the local submatrix with `o_nz` or `o_nnz` (not both).
3273: Consider a processor that owns rows 3, 4 and 5 of a parallel matrix. In
3274: the figure below we depict these three local rows and all columns (0-11).
3276: .vb
3277: 0 1 2 3 4 5 6 7 8 9 10 11
3278: --------------------------
3279: row 3 |o o o d d d o o o o o o
3280: row 4 |o o o d d d o o o o o o
3281: row 5 |o o o d d d o o o o o o
3282: --------------------------
3283: .ve
3285: Thus, any entries in the d locations are stored in the d (diagonal)
3286: submatrix, and any entries in the o locations are stored in the
3287: o (off-diagonal) submatrix. Note that the d and the o submatrices are
3288: stored simply in the `MATSEQBAIJ` format for compressed row storage.
3290: Now `d_nz` should indicate the number of block nonzeros per row in the d matrix,
3291: and `o_nz` should indicate the number of block nonzeros per row in the o matrix.
3292: In general, for PDE problems in which most nonzeros are near the diagonal,
3293: one expects `d_nz` >> `o_nz`.
3295: .seealso: `Mat`, `MatCreate()`, `MatCreateSeqBAIJ()`, `MatSetValues()`, `MatMPIBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocationCSR()`,
3296: `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `MatGetOwnershipRangeColumn()`, `MatGetOwnershipRangesColumn()`, `PetscLayout`
3297: @*/
3298: PetscErrorCode MatCreateBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[], Mat *A)
3299: {
3300: PetscMPIInt size;
3302: PetscFunctionBegin;
3303: PetscCall(MatCreate(comm, A));
3304: PetscCall(MatSetSizes(*A, m, n, M, N));
3305: PetscCallMPI(MPI_Comm_size(comm, &size));
3306: if (size > 1) {
3307: PetscCall(MatSetType(*A, MATMPIBAIJ));
3308: PetscCall(MatMPIBAIJSetPreallocation(*A, bs, d_nz, d_nnz, o_nz, o_nnz));
3309: } else {
3310: PetscCall(MatSetType(*A, MATSEQBAIJ));
3311: PetscCall(MatSeqBAIJSetPreallocation(*A, bs, d_nz, d_nnz));
3312: }
3313: PetscFunctionReturn(PETSC_SUCCESS);
3314: }
3316: static PetscErrorCode MatDuplicate_MPIBAIJ(Mat matin, MatDuplicateOption cpvalues, Mat *newmat)
3317: {
3318: Mat mat;
3319: Mat_MPIBAIJ *a, *oldmat = (Mat_MPIBAIJ *)matin->data;
3320: PetscInt len = 0;
3322: PetscFunctionBegin;
3323: *newmat = NULL;
3324: PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat));
3325: PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N));
3326: PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name));
3328: PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap));
3329: PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap));
3330: if (matin->hash_active) PetscCall(MatSetUp(mat));
3331: else {
3332: mat->factortype = matin->factortype;
3333: mat->preallocated = PETSC_TRUE;
3334: mat->assembled = PETSC_TRUE;
3335: mat->insertmode = NOT_SET_VALUES;
3337: a = (Mat_MPIBAIJ *)mat->data;
3338: mat->rmap->bs = matin->rmap->bs;
3339: a->bs2 = oldmat->bs2;
3340: a->mbs = oldmat->mbs;
3341: a->nbs = oldmat->nbs;
3342: a->Mbs = oldmat->Mbs;
3343: a->Nbs = oldmat->Nbs;
3345: a->size = oldmat->size;
3346: a->rank = oldmat->rank;
3347: a->donotstash = oldmat->donotstash;
3348: a->roworiented = oldmat->roworiented;
3349: a->rowindices = NULL;
3350: a->rowvalues = NULL;
3351: a->getrowactive = PETSC_FALSE;
3352: a->barray = NULL;
3353: a->rstartbs = oldmat->rstartbs;
3354: a->rendbs = oldmat->rendbs;
3355: a->cstartbs = oldmat->cstartbs;
3356: a->cendbs = oldmat->cendbs;
3358: /* hash table stuff */
3359: a->ht = NULL;
3360: a->hd = NULL;
3361: a->ht_size = 0;
3362: a->ht_flag = oldmat->ht_flag;
3363: a->ht_fact = oldmat->ht_fact;
3364: a->ht_total_ct = 0;
3365: a->ht_insert_ct = 0;
3367: PetscCall(PetscArraycpy(a->rangebs, oldmat->rangebs, a->size + 1));
3368: if (oldmat->colmap) {
3369: #if PetscDefined(USE_CTABLE)
3370: PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap));
3371: #else
3372: PetscCall(PetscMalloc1(a->Nbs, &a->colmap));
3373: PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, a->Nbs));
3374: #endif
3375: } else a->colmap = NULL;
3377: if (oldmat->garray && (len = ((Mat_SeqBAIJ *)oldmat->B->data)->nbs)) {
3378: PetscCall(PetscMalloc1(len, &a->garray));
3379: PetscCall(PetscArraycpy(a->garray, oldmat->garray, len));
3380: } else a->garray = NULL;
3382: PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)matin), matin->rmap->bs, &mat->bstash));
3383: PetscCall(VecDuplicate(oldmat->lvec, &a->lvec));
3384: PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx));
3386: PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A));
3387: PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B));
3388: }
3389: PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist));
3390: *newmat = mat;
3391: PetscFunctionReturn(PETSC_SUCCESS);
3392: }
3394: /* Used for both MPIBAIJ and MPISBAIJ matrices */
3395: PetscErrorCode MatLoad_MPIBAIJ_Binary(Mat mat, PetscViewer viewer)
3396: {
3397: PetscInt header[4], M, N, nz, bs, m, n, mbs, nbs, rows, cols, sum, i, j, k;
3398: PetscInt *rowidxs, *colidxs, rs, cs, ce;
3399: PetscScalar *matvals;
3401: PetscFunctionBegin;
3402: PetscCall(PetscViewerSetUp(viewer));
3404: /* read in matrix header */
3405: PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
3406: PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
3407: M = header[1];
3408: N = header[2];
3409: nz = header[3];
3410: PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
3411: PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
3412: PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix stored in special format on disk, cannot load as MPIBAIJ");
3414: /* set block sizes from the viewer's .info file */
3415: PetscCall(MatLoad_Binary_BlockSizes(mat, viewer));
3416: /* set local sizes if not set already */
3417: if (mat->rmap->n < 0 && M == N) mat->rmap->n = mat->cmap->n;
3418: if (mat->cmap->n < 0 && M == N) mat->cmap->n = mat->rmap->n;
3419: /* set global sizes if not set already */
3420: if (mat->rmap->N < 0) mat->rmap->N = M;
3421: if (mat->cmap->N < 0) mat->cmap->N = N;
3422: PetscCall(PetscLayoutSetUp(mat->rmap));
3423: PetscCall(PetscLayoutSetUp(mat->cmap));
3425: /* check if the matrix sizes are correct */
3426: PetscCall(MatGetSize(mat, &rows, &cols));
3427: PetscCheck(M == rows && N == cols, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols);
3428: PetscCall(MatGetBlockSize(mat, &bs));
3429: PetscCall(MatGetLocalSize(mat, &m, &n));
3430: PetscCall(PetscLayoutGetRange(mat->rmap, &rs, NULL));
3431: PetscCall(PetscLayoutGetRange(mat->cmap, &cs, &ce));
3432: mbs = m / bs;
3433: nbs = n / bs;
3435: /* read in row lengths and build row indices */
3436: PetscCall(PetscMalloc1(m + 1, &rowidxs));
3437: PetscCall(PetscViewerBinaryReadAll(viewer, rowidxs + 1, m, PETSC_DECIDE, M, PETSC_INT));
3438: rowidxs[0] = 0;
3439: for (i = 0; i < m; i++) rowidxs[i + 1] += rowidxs[i];
3440: PetscCallMPI(MPIU_Allreduce(&rowidxs[m], &sum, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)viewer)));
3441: PetscCheck(sum == nz, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT, nz, sum);
3443: /* read in column indices and matrix values */
3444: PetscCall(PetscMalloc2(rowidxs[m], &colidxs, rowidxs[m], &matvals));
3445: PetscCall(PetscViewerBinaryReadAll(viewer, colidxs, rowidxs[m], PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
3446: PetscCall(PetscViewerBinaryReadAll(viewer, matvals, rowidxs[m], PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
3448: { /* preallocate matrix storage */
3449: PetscBT bt; /* helper bit set to count diagonal nonzeros */
3450: PetscHSetI ht; /* helper hash set to count off-diagonal nonzeros */
3451: PetscBool sbaij, done;
3452: PetscInt *d_nnz, *o_nnz;
3454: PetscCall(PetscBTCreate(nbs, &bt));
3455: PetscCall(PetscHSetICreate(&ht));
3456: PetscCall(PetscCalloc2(mbs, &d_nnz, mbs, &o_nnz));
3457: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPISBAIJ, &sbaij));
3458: for (i = 0; i < mbs; i++) {
3459: PetscCall(PetscBTMemzero(nbs, bt));
3460: PetscCall(PetscHSetIClear(ht));
3461: for (k = 0; k < bs; k++) {
3462: PetscInt row = bs * i + k;
3463: for (j = rowidxs[row]; j < rowidxs[row + 1]; j++) {
3464: PetscInt col = colidxs[j];
3465: if (!sbaij || col >= row) {
3466: if (col >= cs && col < ce) {
3467: if (!PetscBTLookupSet(bt, (col - cs) / bs)) d_nnz[i]++;
3468: } else {
3469: PetscCall(PetscHSetIQueryAdd(ht, col / bs, &done));
3470: if (done) o_nnz[i]++;
3471: }
3472: }
3473: }
3474: }
3475: }
3476: PetscCall(PetscBTDestroy(&bt));
3477: PetscCall(PetscHSetIDestroy(&ht));
3478: PetscCall(MatMPIBAIJSetPreallocation(mat, bs, 0, d_nnz, 0, o_nnz));
3479: PetscCall(MatMPISBAIJSetPreallocation(mat, bs, 0, d_nnz, 0, o_nnz));
3480: PetscCall(PetscFree2(d_nnz, o_nnz));
3481: }
3483: /* store matrix values */
3484: for (i = 0; i < m; i++) {
3485: PetscInt row = rs + i, s = rowidxs[i], e = rowidxs[i + 1];
3486: PetscUseTypeMethod(mat, setvalues, 1, &row, e - s, colidxs + s, matvals + s, INSERT_VALUES);
3487: }
3489: PetscCall(PetscFree(rowidxs));
3490: PetscCall(PetscFree2(colidxs, matvals));
3491: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3492: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3493: PetscFunctionReturn(PETSC_SUCCESS);
3494: }
3496: PetscErrorCode MatLoad_MPIBAIJ(Mat mat, PetscViewer viewer)
3497: {
3498: PetscBool isbinary;
3500: PetscFunctionBegin;
3501: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
3502: PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not yet supported for reading %s matrices", ((PetscObject)viewer)->type_name, ((PetscObject)mat)->type_name);
3503: PetscCall(MatLoad_MPIBAIJ_Binary(mat, viewer));
3504: PetscFunctionReturn(PETSC_SUCCESS);
3505: }
3507: /*@
3508: MatMPIBAIJSetHashTableFactor - Sets the factor required to compute the size of the matrices hash table
3510: Input Parameters:
3511: + mat - the matrix
3512: - fact - factor
3514: Options Database Key:
3515: . -mat_use_hash_table fact - provide the factor
3517: Level: advanced
3519: .seealso: `Mat`, `MATMPIBAIJ`, `MatSetOption()`
3520: @*/
3521: PetscErrorCode MatMPIBAIJSetHashTableFactor(Mat mat, PetscReal fact)
3522: {
3523: PetscFunctionBegin;
3524: PetscTryMethod(mat, "MatSetHashTableFactor_C", (Mat, PetscReal), (mat, fact));
3525: PetscFunctionReturn(PETSC_SUCCESS);
3526: }
3528: PetscErrorCode MatSetHashTableFactor_MPIBAIJ(Mat mat, PetscReal fact)
3529: {
3530: Mat_MPIBAIJ *baij;
3532: PetscFunctionBegin;
3533: baij = (Mat_MPIBAIJ *)mat->data;
3534: baij->ht_fact = fact;
3535: PetscFunctionReturn(PETSC_SUCCESS);
3536: }
3538: /*@C
3539: MatMPIBAIJGetSeqBAIJ - Get the on-process (diagonal block) and off-process (off-diagonal block) `MATSEQBAIJ`
3540: matrices that make up an `MATMPIBAIJ` matrix, together with the local-to-global column map for the off-diagonal block.
3542: Not Collective
3544: Input Parameter:
3545: . A - the `MATMPIBAIJ` matrix
3547: Output Parameters:
3548: + Ad - the diagonal block `MATSEQBAIJ`, or `NULL` if not needed
3549: . Ao - the off-diagonal block `MATSEQBAIJ`, or `NULL` if not needed
3550: - colmap - the local-to-global column index map for `Ao`, or `NULL` if not needed
3552: Level: advanced
3554: .seealso: `Mat`, `MATMPIBAIJ`, `MATSEQBAIJ`, `MatMPIAIJGetSeqAIJ()`
3555: @*/
3556: PetscErrorCode MatMPIBAIJGetSeqBAIJ(Mat A, Mat *Ad, Mat *Ao, const PetscInt *colmap[])
3557: {
3558: Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
3559: PetscBool flg;
3561: PetscFunctionBegin;
3562: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIBAIJ, &flg));
3563: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "This function requires a MATMPIBAIJ matrix as input");
3564: if (Ad) *Ad = a->A;
3565: if (Ao) *Ao = a->B;
3566: if (colmap) *colmap = a->garray;
3567: PetscFunctionReturn(PETSC_SUCCESS);
3568: }
3570: /*
3571: Special version for direct calls from Fortran (to eliminate two function call overheads
3572: */
3573: #if PetscDefined(HAVE_FORTRAN_CAPS)
3574: #define matmpibaijsetvaluesblocked_ MATMPIBAIJSETVALUESBLOCKED
3575: #elif !PetscDefined(HAVE_FORTRAN_UNDERSCORE)
3576: #define matmpibaijsetvaluesblocked_ matmpibaijsetvaluesblocked
3577: #endif
3579: // PetscClangLinter pragma disable: -fdoc-synopsis-matching-symbol-name
3580: /*@C
3581: MatMPIBAIJSetValuesBlocked - Direct Fortran call to replace call to `MatSetValuesBlocked()`
3583: Collective
3585: Input Parameters:
3586: + matin - the matrix
3587: . min - number of input rows
3588: . im - input rows
3589: . nin - number of input columns
3590: . in - input columns
3591: . v - numerical values input
3592: - addvin - `INSERT_VALUES` or `ADD_VALUES`
3594: Level: advanced
3596: Developer Notes:
3597: This has a complete copy of `MatSetValuesBlocked_MPIBAIJ()` which is terrible code un-reuse.
3599: .seealso: `Mat`, `MatSetValuesBlocked()`
3600: @*/
3601: PETSC_EXTERN PetscErrorCode matmpibaijsetvaluesblocked_(Mat *matin, PetscInt *min, const PetscInt im[], PetscInt *nin, const PetscInt in[], const MatScalar v[], InsertMode *addvin)
3602: {
3603: /* convert input arguments to C version */
3604: Mat mat = *matin;
3605: PetscInt m = *min, n = *nin;
3606: InsertMode addv = *addvin;
3608: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
3609: const MatScalar *value;
3610: MatScalar *barray = baij->barray;
3611: PetscBool roworiented = baij->roworiented;
3612: PetscInt i, j, ii, jj, row, col, rstart = baij->rstartbs;
3613: PetscInt rend = baij->rendbs, cstart = baij->cstartbs, stepval;
3614: PetscInt cend = baij->cendbs, bs = mat->rmap->bs, bs2 = baij->bs2;
3616: PetscFunctionBegin;
3617: /* tasks normally handled by MatSetValuesBlocked() */
3618: if (mat->insertmode == NOT_SET_VALUES) mat->insertmode = addv;
3619: else PetscCheck(mat->insertmode == addv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot mix add values and insert values");
3620: PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3621: if (mat->assembled) {
3622: mat->was_assembled = PETSC_TRUE;
3623: mat->assembled = PETSC_FALSE;
3624: }
3625: PetscCall(PetscLogEventBegin(MAT_SetValues, mat, 0, 0, 0));
3627: if (!barray) {
3628: PetscCall(PetscMalloc1(bs2, &barray));
3629: baij->barray = barray;
3630: }
3632: if (roworiented) stepval = (n - 1) * bs;
3633: else stepval = (m - 1) * bs;
3635: for (i = 0; i < m; i++) {
3636: if (im[i] < 0) continue;
3637: PetscCheck(im[i] < baij->Mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large, row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], baij->Mbs - 1);
3638: if (im[i] >= rstart && im[i] < rend) {
3639: row = im[i] - rstart;
3640: for (j = 0; j < n; j++) {
3641: /* If NumCol = 1 then a copy is not required */
3642: if ((roworiented) && (n == 1)) {
3643: barray = (MatScalar *)v + i * bs2;
3644: } else if ((!roworiented) && (m == 1)) {
3645: barray = (MatScalar *)v + j * bs2;
3646: } else { /* Here a copy is required */
3647: if (roworiented) {
3648: value = v + i * (stepval + bs) * bs + j * bs;
3649: } else {
3650: value = v + j * (stepval + bs) * bs + i * bs;
3651: }
3652: for (ii = 0; ii < bs; ii++, value += stepval) {
3653: for (jj = 0; jj < bs; jj++) *barray++ = *value++;
3654: }
3655: barray -= bs2;
3656: }
3658: if (in[j] >= cstart && in[j] < cend) {
3659: col = in[j] - cstart;
3660: PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->A, row, col, barray, addv, im[i], in[j]));
3661: } else if (in[j] < 0) {
3662: continue;
3663: } else {
3664: PetscCheck(in[j] < baij->Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large, col %" PetscInt_FMT " max %" PetscInt_FMT, in[j], baij->Nbs - 1);
3665: if (mat->was_assembled) {
3666: if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
3668: #if PetscDefined(USE_CTABLE)
3669: if (PetscDefined(USE_DEBUG)) {
3670: PetscInt data;
3671: PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] + 1, 0, &data));
3672: PetscCheck((data - 1) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect colmap");
3673: }
3674: #else
3675: if (PetscDefined(USE_DEBUG)) PetscCheck((baij->colmap[in[j]] - 1) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect colmap");
3676: #endif
3677: #if PetscDefined(USE_CTABLE)
3678: PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] + 1, 0, &col));
3679: col = (col - 1) / bs;
3680: #else
3681: col = (baij->colmap[in[j]] - 1) / bs;
3682: #endif
3683: if (col < 0 && !((Mat_SeqBAIJ *)baij->A->data)->nonew) {
3684: PetscCall(MatDisAssemble_MPIBAIJ(mat));
3685: col = in[j];
3686: }
3687: } else col = in[j];
3688: PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->B, row, col, barray, addv, im[i], in[j]));
3689: }
3690: }
3691: } else {
3692: if (!baij->donotstash) {
3693: if (roworiented) {
3694: PetscCall(MatStashValuesRowBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
3695: } else {
3696: PetscCall(MatStashValuesColBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
3697: }
3698: }
3699: }
3700: }
3702: /* task normally handled by MatSetValuesBlocked() */
3703: PetscCall(PetscLogEventEnd(MAT_SetValues, mat, 0, 0, 0));
3704: PetscFunctionReturn(PETSC_SUCCESS);
3705: }
3707: /*@
3708: MatCreateMPIBAIJWithArrays - creates a `MATMPIBAIJ` matrix using arrays that contain in standard block CSR format for the local rows.
3710: Collective
3712: Input Parameters:
3713: + comm - MPI communicator
3714: . bs - the block size, only a block size of 1 is supported
3715: . m - number of local rows (Cannot be `PETSC_DECIDE`)
3716: . n - This value should be the same as the local size used in creating the
3717: x vector for the matrix-vector product $ y = Ax $. (or `PETSC_DECIDE` to have
3718: calculated if `N` is given) For square matrices `n` is almost always `m`.
3719: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
3720: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
3721: . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that rowth block row of the matrix
3722: . j - column indices
3723: - a - matrix values
3725: Output Parameter:
3726: . mat - the matrix
3728: Level: intermediate
3730: Notes:
3731: The `i`, `j`, and `a` arrays ARE copied by this routine into the internal format used by PETSc;
3732: thus you CANNOT change the matrix entries by changing the values of a[] after you have
3733: called this routine. Use `MatCreateMPIAIJWithSplitArrays()` to avoid needing to copy the arrays.
3735: The order of the entries in values is the same as the block compressed sparse row storage format; that is, it is
3736: the same as a three dimensional array in Fortran values(bs,bs,nnz) that contains the first column of the first
3737: block, followed by the second column of the first block etc etc. That is, the blocks are contiguous in memory
3738: with column-major ordering within blocks.
3740: The `i` and `j` indices are 0 based, and `i` indices are indices corresponding to the local `j` array.
3742: .seealso: `Mat`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatMPIAIJSetPreallocation()`, `MatMPIAIJSetPreallocationCSR()`,
3743: `MATMPIAIJ`, `MatCreateAIJ()`, `MatCreateMPIAIJWithSplitArrays()`
3744: @*/
3745: PetscErrorCode MatCreateMPIBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, const PetscInt i[], const PetscInt j[], const PetscScalar a[], Mat *mat)
3746: {
3747: PetscFunctionBegin;
3748: PetscCheck(!i[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");
3749: PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "local number of rows (m) cannot be PETSC_DECIDE, or negative");
3750: PetscCall(MatCreate(comm, mat));
3751: PetscCall(MatSetSizes(*mat, m, n, M, N));
3752: PetscCall(MatSetType(*mat, MATMPIBAIJ));
3753: PetscCall(MatSetBlockSize(*mat, bs));
3754: PetscCall(MatSetUp(*mat));
3755: PetscCall(MatSetOption(*mat, MAT_ROW_ORIENTED, PETSC_FALSE));
3756: PetscCall(MatMPIBAIJSetPreallocationCSR(*mat, bs, i, j, a));
3757: PetscCall(MatSetOption(*mat, MAT_ROW_ORIENTED, PETSC_TRUE));
3758: PetscFunctionReturn(PETSC_SUCCESS);
3759: }
3761: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_MPIBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3762: {
3763: PetscInt m, N, i, rstart, nnz, Ii, bs, cbs;
3764: PetscInt *indx;
3765: PetscScalar *values;
3767: PetscFunctionBegin;
3768: PetscCall(MatGetSize(inmat, &m, &N));
3769: if (scall == MAT_INITIAL_MATRIX) { /* symbolic phase */
3770: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)inmat->data;
3771: PetscInt *dnz, *onz, mbs, Nbs, nbs;
3772: PetscInt *bindx, rmax = a->rmax, j;
3773: PetscMPIInt rank, size;
3775: PetscCall(MatGetBlockSizes(inmat, &bs, &cbs));
3776: mbs = m / bs;
3777: Nbs = N / cbs;
3778: if (n == PETSC_DECIDE) PetscCall(PetscSplitOwnershipBlock(comm, cbs, &n, &N));
3779: nbs = n / cbs;
3781: PetscCall(PetscMalloc1(rmax, &bindx));
3782: MatPreallocateBegin(comm, mbs, nbs, dnz, onz); /* inline function, output __end and __rstart are used below */
3784: PetscCallMPI(MPI_Comm_rank(comm, &rank));
3785: PetscCallMPI(MPI_Comm_size(comm, &size));
3786: if (rank == size - 1) {
3787: /* Check sum(nbs) = Nbs */
3788: PetscCheck(__end == Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Sum of local block columns %" PetscInt_FMT " != global block columns %" PetscInt_FMT, __end, Nbs);
3789: }
3791: rstart = __rstart; /* block rstart of *outmat; see inline function MatPreallocateBegin */
3792: for (i = 0; i < mbs; i++) {
3793: PetscCall(MatGetRow_SeqBAIJ(inmat, i * bs, &nnz, &indx, NULL)); /* non-blocked nnz and indx */
3794: nnz = nnz / bs;
3795: for (j = 0; j < nnz; j++) bindx[j] = indx[j * bs] / bs;
3796: PetscCall(MatPreallocateSet(i + rstart, nnz, bindx, dnz, onz));
3797: PetscCall(MatRestoreRow_SeqBAIJ(inmat, i * bs, &nnz, &indx, NULL));
3798: }
3799: PetscCall(PetscFree(bindx));
3801: PetscCall(MatCreate(comm, outmat));
3802: PetscCall(MatSetSizes(*outmat, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
3803: PetscCall(MatSetBlockSizes(*outmat, bs, cbs));
3804: PetscCall(MatSetType(*outmat, MATBAIJ));
3805: PetscCall(MatSeqBAIJSetPreallocation(*outmat, bs, 0, dnz));
3806: PetscCall(MatMPIBAIJSetPreallocation(*outmat, bs, 0, dnz, 0, onz));
3807: MatPreallocateEnd(dnz, onz);
3808: PetscCall(MatSetOption(*outmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
3809: }
3811: /* numeric phase */
3812: PetscCall(MatGetBlockSizes(inmat, &bs, &cbs));
3813: PetscCall(MatGetOwnershipRange(*outmat, &rstart, NULL));
3815: for (i = 0; i < m; i++) {
3816: PetscCall(MatGetRow_SeqBAIJ(inmat, i, &nnz, &indx, &values));
3817: Ii = i + rstart;
3818: PetscCall(MatSetValues(*outmat, 1, &Ii, nnz, indx, values, INSERT_VALUES));
3819: PetscCall(MatRestoreRow_SeqBAIJ(inmat, i, &nnz, &indx, &values));
3820: }
3821: PetscCall(MatAssemblyBegin(*outmat, MAT_FINAL_ASSEMBLY));
3822: PetscCall(MatAssemblyEnd(*outmat, MAT_FINAL_ASSEMBLY));
3823: PetscFunctionReturn(PETSC_SUCCESS);
3824: }