Actual source code: mpisbaij.c
1: #include <../src/mat/impls/baij/mpi/mpibaij.h>
2: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
3: #include <../src/mat/impls/sbaij/seq/sbaij.h>
4: #include <petscblaslapack.h>
5: #include <petscsf.h>
7: static PetscErrorCode MatDestroy_MPISBAIJ(Mat mat)
8: {
9: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
11: PetscFunctionBegin;
12: PetscCall(PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ",Cols=%" PetscInt_FMT, mat->rmap->N, mat->cmap->N));
13: PetscCall(MatStashDestroy_Private(&mat->stash));
14: PetscCall(MatStashDestroy_Private(&mat->bstash));
15: PetscCall(MatDestroy(&baij->A));
16: PetscCall(MatDestroy(&baij->B));
17: #if defined(PETSC_USE_CTABLE)
18: PetscCall(PetscHMapIDestroy(&baij->colmap));
19: #else
20: PetscCall(PetscFree(baij->colmap));
21: #endif
22: PetscCall(PetscFree(baij->garray));
23: PetscCall(VecDestroy(&baij->lvec));
24: PetscCall(VecScatterDestroy(&baij->Mvctx));
25: PetscCall(VecDestroy(&baij->slvec0));
26: PetscCall(VecDestroy(&baij->slvec0b));
27: PetscCall(VecDestroy(&baij->slvec1));
28: PetscCall(VecDestroy(&baij->slvec1a));
29: PetscCall(VecDestroy(&baij->slvec1b));
30: PetscCall(VecScatterDestroy(&baij->sMvctx));
31: PetscCall(PetscFree2(baij->rowvalues, baij->rowindices));
32: PetscCall(PetscFree(baij->barray));
33: PetscCall(PetscFree(baij->hd));
34: PetscCall(VecDestroy(&baij->diag));
35: PetscCall(VecDestroy(&baij->bb1));
36: PetscCall(VecDestroy(&baij->xx1));
37: #if defined(PETSC_USE_REAL_MAT_SINGLE)
38: PetscCall(PetscFree(baij->setvaluescopy));
39: #endif
40: PetscCall(PetscFree(baij->in_loc));
41: PetscCall(PetscFree(baij->v_loc));
42: PetscCall(PetscFree(baij->rangebs));
43: PetscCall(PetscFree(mat->data));
45: PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
46: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatStoreValues_C", NULL));
47: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatRetrieveValues_C", NULL));
48: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPISBAIJSetPreallocation_C", NULL));
49: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPISBAIJSetPreallocationCSR_C", NULL));
50: #if defined(PETSC_HAVE_ELEMENTAL)
51: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisbaij_elemental_C", NULL));
52: #endif
53: #if defined(PETSC_HAVE_SCALAPACK) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
54: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisbaij_scalapack_C", NULL));
55: #endif
56: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisbaij_mpiaij_C", NULL));
57: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisbaij_mpibaij_C", NULL));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: /* defines MatSetValues_MPI_Hash(), MatAssemblyBegin_MPI_Hash(), MatAssemblyEnd_MPI_Hash(), MatSetUp_MPI_Hash() */
62: #define TYPE SBAIJ
63: #define TYPE_SBAIJ
64: #include "../src/mat/impls/aij/mpi/mpihashmat.h"
65: #undef TYPE
66: #undef TYPE_SBAIJ
68: #if defined(PETSC_HAVE_ELEMENTAL)
69: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Elemental(Mat, MatType, MatReuse, Mat *);
70: #endif
71: #if defined(PETSC_HAVE_SCALAPACK) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
72: PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *);
73: #endif
75: /* This could be moved to matimpl.h */
76: static PetscErrorCode MatPreallocateWithMats_Private(Mat B, PetscInt nm, Mat X[], PetscBool symm[], PetscBool fill)
77: {
78: Mat preallocator;
79: PetscInt r, rstart, rend;
80: PetscInt bs, i, m, n, M, N;
81: PetscBool cong = PETSC_TRUE;
83: PetscFunctionBegin;
86: for (i = 0; i < nm; i++) {
88: PetscCall(PetscLayoutCompare(B->rmap, X[i]->rmap, &cong));
89: PetscCheck(cong, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Not for different layouts");
90: }
92: PetscCall(MatGetBlockSize(B, &bs));
93: PetscCall(MatGetSize(B, &M, &N));
94: PetscCall(MatGetLocalSize(B, &m, &n));
95: PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &preallocator));
96: PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
97: PetscCall(MatSetBlockSize(preallocator, bs));
98: PetscCall(MatSetSizes(preallocator, m, n, M, N));
99: PetscCall(MatSetUp(preallocator));
100: PetscCall(MatGetOwnershipRange(preallocator, &rstart, &rend));
101: for (r = rstart; r < rend; ++r) {
102: PetscInt ncols;
103: const PetscInt *row;
104: const PetscScalar *vals;
106: for (i = 0; i < nm; i++) {
107: PetscCall(MatGetRow(X[i], r, &ncols, &row, &vals));
108: PetscCall(MatSetValues(preallocator, 1, &r, ncols, row, vals, INSERT_VALUES));
109: if (symm && symm[i]) PetscCall(MatSetValues(preallocator, ncols, row, 1, &r, vals, INSERT_VALUES));
110: PetscCall(MatRestoreRow(X[i], r, &ncols, &row, &vals));
111: }
112: }
113: PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
114: PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
115: PetscCall(MatPreallocatorPreallocate(preallocator, fill, B));
116: PetscCall(MatDestroy(&preallocator));
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
121: {
122: Mat B;
123: PetscInt r;
125: PetscFunctionBegin;
126: if (reuse != MAT_REUSE_MATRIX) {
127: PetscBool symm = PETSC_TRUE, isdense;
128: PetscInt bs;
130: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
131: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
132: PetscCall(MatSetType(B, newtype));
133: PetscCall(MatGetBlockSize(A, &bs));
134: PetscCall(MatSetBlockSize(B, bs));
135: PetscCall(PetscLayoutSetUp(B->rmap));
136: PetscCall(PetscLayoutSetUp(B->cmap));
137: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, ""));
138: if (!isdense) {
139: PetscCall(MatGetRowUpperTriangular(A));
140: PetscCall(MatPreallocateWithMats_Private(B, 1, &A, &symm, PETSC_TRUE));
141: PetscCall(MatRestoreRowUpperTriangular(A));
142: } else {
143: PetscCall(MatSetUp(B));
144: }
145: } else {
146: B = *newmat;
147: PetscCall(MatZeroEntries(B));
148: }
150: PetscCall(MatGetRowUpperTriangular(A));
151: for (r = A->rmap->rstart; r < A->rmap->rend; r++) {
152: PetscInt ncols;
153: const PetscInt *row;
154: const PetscScalar *vals;
156: PetscCall(MatGetRow(A, r, &ncols, &row, &vals));
157: PetscCall(MatSetValues(B, 1, &r, ncols, row, vals, INSERT_VALUES));
158: #if defined(PETSC_USE_COMPLEX)
159: if (A->hermitian == PETSC_BOOL3_TRUE) {
160: PetscInt i;
161: for (i = 0; i < ncols; i++) PetscCall(MatSetValue(B, row[i], r, PetscConj(vals[i]), INSERT_VALUES));
162: } else {
163: PetscCall(MatSetValues(B, ncols, row, 1, &r, vals, INSERT_VALUES));
164: }
165: #else
166: PetscCall(MatSetValues(B, ncols, row, 1, &r, vals, INSERT_VALUES));
167: #endif
168: PetscCall(MatRestoreRow(A, r, &ncols, &row, &vals));
169: }
170: PetscCall(MatRestoreRowUpperTriangular(A));
171: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
172: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
174: if (reuse == MAT_INPLACE_MATRIX) {
175: PetscCall(MatHeaderReplace(A, &B));
176: } else {
177: *newmat = B;
178: }
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: static PetscErrorCode MatStoreValues_MPISBAIJ(Mat mat)
183: {
184: Mat_MPISBAIJ *aij = (Mat_MPISBAIJ *)mat->data;
186: PetscFunctionBegin;
187: PetscCall(MatStoreValues(aij->A));
188: PetscCall(MatStoreValues(aij->B));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode MatRetrieveValues_MPISBAIJ(Mat mat)
193: {
194: Mat_MPISBAIJ *aij = (Mat_MPISBAIJ *)mat->data;
196: PetscFunctionBegin;
197: PetscCall(MatRetrieveValues(aij->A));
198: PetscCall(MatRetrieveValues(aij->B));
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: #define MatSetValues_SeqSBAIJ_A_Private(row, col, value, addv, orow, ocol) \
203: do { \
204: brow = row / bs; \
205: rp = aj + ai[brow]; \
206: ap = aa + bs2 * ai[brow]; \
207: rmax = aimax[brow]; \
208: nrow = ailen[brow]; \
209: bcol = col / bs; \
210: ridx = row % bs; \
211: cidx = col % bs; \
212: low = 0; \
213: high = nrow; \
214: while (high - low > 3) { \
215: t = (low + high) / 2; \
216: if (rp[t] > bcol) high = t; \
217: else low = t; \
218: } \
219: for (_i = low; _i < high; _i++) { \
220: if (rp[_i] > bcol) break; \
221: if (rp[_i] == bcol) { \
222: bap = ap + bs2 * _i + bs * cidx + ridx; \
223: if (addv == ADD_VALUES) *bap += value; \
224: else *bap = value; \
225: goto a_noinsert; \
226: } \
227: } \
228: if (a->nonew == 1) goto a_noinsert; \
229: 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); \
230: MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, aimax, a->nonew, MatScalar); \
231: N = nrow++ - 1; \
232: /* shift up all the later entries in this row */ \
233: PetscCall(PetscArraymove(rp + _i + 1, rp + _i, N - _i + 1)); \
234: PetscCall(PetscArraymove(ap + bs2 * (_i + 1), ap + bs2 * _i, bs2 * (N - _i + 1))); \
235: PetscCall(PetscArrayzero(ap + bs2 * _i, bs2)); \
236: rp[_i] = bcol; \
237: ap[bs2 * _i + bs * cidx + ridx] = value; \
238: a_noinsert:; \
239: ailen[brow] = nrow; \
240: } while (0)
242: #define MatSetValues_SeqSBAIJ_B_Private(row, col, value, addv, orow, ocol) \
243: do { \
244: brow = row / bs; \
245: rp = bj + bi[brow]; \
246: ap = ba + bs2 * bi[brow]; \
247: rmax = bimax[brow]; \
248: nrow = bilen[brow]; \
249: bcol = col / bs; \
250: ridx = row % bs; \
251: cidx = col % bs; \
252: low = 0; \
253: high = nrow; \
254: while (high - low > 3) { \
255: t = (low + high) / 2; \
256: if (rp[t] > bcol) high = t; \
257: else low = t; \
258: } \
259: for (_i = low; _i < high; _i++) { \
260: if (rp[_i] > bcol) break; \
261: if (rp[_i] == bcol) { \
262: bap = ap + bs2 * _i + bs * cidx + ridx; \
263: if (addv == ADD_VALUES) *bap += value; \
264: else *bap = value; \
265: goto b_noinsert; \
266: } \
267: } \
268: if (b->nonew == 1) goto b_noinsert; \
269: 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); \
270: MatSeqXAIJReallocateAIJ(B, b->mbs, bs2, nrow, brow, bcol, rmax, ba, bi, bj, rp, ap, bimax, b->nonew, MatScalar); \
271: N = nrow++ - 1; \
272: /* shift up all the later entries in this row */ \
273: PetscCall(PetscArraymove(rp + _i + 1, rp + _i, N - _i + 1)); \
274: PetscCall(PetscArraymove(ap + bs2 * (_i + 1), ap + bs2 * _i, bs2 * (N - _i + 1))); \
275: PetscCall(PetscArrayzero(ap + bs2 * _i, bs2)); \
276: rp[_i] = bcol; \
277: ap[bs2 * _i + bs * cidx + ridx] = value; \
278: b_noinsert:; \
279: bilen[brow] = nrow; \
280: } while (0)
282: /* Only add/insert a(i,j) with i<=j (blocks).
283: Any a(i,j) with i>j input by user is ignored or generates an error
284: */
285: static PetscErrorCode MatSetValues_MPISBAIJ(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
286: {
287: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
288: MatScalar value;
289: PetscBool roworiented = baij->roworiented;
290: PetscInt i, j, row, col;
291: PetscInt rstart_orig = mat->rmap->rstart;
292: PetscInt rend_orig = mat->rmap->rend, cstart_orig = mat->cmap->rstart;
293: PetscInt cend_orig = mat->cmap->rend, bs = mat->rmap->bs;
295: /* Some Variables required in the macro */
296: Mat A = baij->A;
297: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
298: PetscInt *aimax = a->imax, *ai = a->i, *ailen = a->ilen, *aj = a->j;
299: MatScalar *aa = a->a;
301: Mat B = baij->B;
302: Mat_SeqBAIJ *b = (Mat_SeqBAIJ *)B->data;
303: PetscInt *bimax = b->imax, *bi = b->i, *bilen = b->ilen, *bj = b->j;
304: MatScalar *ba = b->a;
306: PetscInt *rp, ii, nrow, _i, rmax, N, brow, bcol;
307: PetscInt low, high, t, ridx, cidx, bs2 = a->bs2;
308: MatScalar *ap, *bap;
310: /* for stash */
311: PetscInt n_loc, *in_loc = NULL;
312: MatScalar *v_loc = NULL;
314: PetscFunctionBegin;
315: if (!baij->donotstash) {
316: if (n > baij->n_loc) {
317: PetscCall(PetscFree(baij->in_loc));
318: PetscCall(PetscFree(baij->v_loc));
319: PetscCall(PetscMalloc1(n, &baij->in_loc));
320: PetscCall(PetscMalloc1(n, &baij->v_loc));
322: baij->n_loc = n;
323: }
324: in_loc = baij->in_loc;
325: v_loc = baij->v_loc;
326: }
328: for (i = 0; i < m; i++) {
329: if (im[i] < 0) continue;
330: 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);
331: if (im[i] >= rstart_orig && im[i] < rend_orig) { /* this processor entry */
332: row = im[i] - rstart_orig; /* local row index */
333: for (j = 0; j < n; j++) {
334: if (im[i] / bs > in[j] / bs) {
335: PetscCheck(a->ignore_ltriangular, PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
336: continue; /* ignore lower triangular blocks */
337: }
338: if (in[j] >= cstart_orig && in[j] < cend_orig) { /* diag entry (A) */
339: col = in[j] - cstart_orig; /* local col index */
340: brow = row / bs;
341: bcol = col / bs;
342: if (brow > bcol) continue; /* ignore lower triangular blocks of A */
343: if (roworiented) value = v[i * n + j];
344: else value = v[i + j * m];
345: MatSetValues_SeqSBAIJ_A_Private(row, col, value, addv, im[i], in[j]);
346: /* PetscCall(MatSetValues_SeqBAIJ(baij->A,1,&row,1,&col,&value,addv)); */
347: } else if (in[j] < 0) {
348: continue;
349: } else {
350: 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);
351: /* off-diag entry (B) */
352: if (mat->was_assembled) {
353: if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
354: #if defined(PETSC_USE_CTABLE)
355: PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] / bs + 1, 0, &col));
356: col = col - 1;
357: #else
358: col = baij->colmap[in[j] / bs] - 1;
359: #endif
360: if (col < 0 && !((Mat_SeqSBAIJ *)baij->A->data)->nonew) {
361: PetscCall(MatDisAssemble_MPISBAIJ(mat));
362: col = in[j];
363: /* Reinitialize the variables required by MatSetValues_SeqBAIJ_B_Private() */
364: B = baij->B;
365: b = (Mat_SeqBAIJ *)B->data;
366: bimax = b->imax;
367: bi = b->i;
368: bilen = b->ilen;
369: bj = b->j;
370: ba = b->a;
371: } else col += in[j] % bs;
372: } else col = in[j];
373: if (roworiented) value = v[i * n + j];
374: else value = v[i + j * m];
375: MatSetValues_SeqSBAIJ_B_Private(row, col, value, addv, im[i], in[j]);
376: /* PetscCall(MatSetValues_SeqBAIJ(baij->B,1,&row,1,&col,&value,addv)); */
377: }
378: }
379: } else { /* off processor entry */
380: 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]);
381: if (!baij->donotstash) {
382: mat->assembled = PETSC_FALSE;
383: n_loc = 0;
384: for (j = 0; j < n; j++) {
385: if (im[i] / bs > in[j] / bs) continue; /* ignore lower triangular blocks */
386: in_loc[n_loc] = in[j];
387: if (roworiented) {
388: v_loc[n_loc] = v[i * n + j];
389: } else {
390: v_loc[n_loc] = v[j * m + i];
391: }
392: n_loc++;
393: }
394: PetscCall(MatStashValuesRow_Private(&mat->stash, im[i], n_loc, in_loc, v_loc, PETSC_FALSE));
395: }
396: }
397: }
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: static inline PetscErrorCode MatSetValuesBlocked_SeqSBAIJ_Inlined(Mat A, PetscInt row, PetscInt col, const PetscScalar v[], InsertMode is, PetscInt orow, PetscInt ocol)
402: {
403: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
404: PetscInt *rp, low, high, t, ii, jj, nrow, i, rmax, N;
405: PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen;
406: PetscInt *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs;
407: PetscBool roworiented = a->roworiented;
408: const PetscScalar *value = v;
409: MatScalar *ap, *aa = a->a, *bap;
411: PetscFunctionBegin;
412: if (col < row) {
413: PetscCheck(a->ignore_ltriangular, PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
414: PetscFunctionReturn(PETSC_SUCCESS); /* ignore lower triangular block */
415: }
416: rp = aj + ai[row];
417: ap = aa + bs2 * ai[row];
418: rmax = imax[row];
419: nrow = ailen[row];
420: value = v;
421: low = 0;
422: high = nrow;
424: while (high - low > 7) {
425: t = (low + high) / 2;
426: if (rp[t] > col) high = t;
427: else low = t;
428: }
429: for (i = low; i < high; i++) {
430: if (rp[i] > col) break;
431: if (rp[i] == col) {
432: bap = ap + bs2 * i;
433: if (roworiented) {
434: if (is == ADD_VALUES) {
435: for (ii = 0; ii < bs; ii++) {
436: for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
437: }
438: } else {
439: for (ii = 0; ii < bs; ii++) {
440: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
441: }
442: }
443: } else {
444: if (is == ADD_VALUES) {
445: for (ii = 0; ii < bs; ii++) {
446: for (jj = 0; jj < bs; jj++) *bap++ += *value++;
447: }
448: } else {
449: for (ii = 0; ii < bs; ii++) {
450: for (jj = 0; jj < bs; jj++) *bap++ = *value++;
451: }
452: }
453: }
454: goto noinsert2;
455: }
456: }
457: if (nonew == 1) goto noinsert2;
458: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new block index nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", orow, ocol);
459: MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
460: N = nrow++ - 1;
461: high++;
462: /* shift up all the later entries in this row */
463: PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
464: PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
465: rp[i] = col;
466: bap = ap + bs2 * i;
467: if (roworiented) {
468: for (ii = 0; ii < bs; ii++) {
469: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
470: }
471: } else {
472: for (ii = 0; ii < bs; ii++) {
473: for (jj = 0; jj < bs; jj++) *bap++ = *value++;
474: }
475: }
476: noinsert2:;
477: ailen[row] = nrow;
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*
482: This routine is exactly duplicated in mpibaij.c
483: */
484: static inline PetscErrorCode MatSetValuesBlocked_SeqBAIJ_Inlined(Mat A, PetscInt row, PetscInt col, const PetscScalar v[], InsertMode is, PetscInt orow, PetscInt ocol)
485: {
486: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
487: PetscInt *rp, low, high, t, ii, jj, nrow, i, rmax, N;
488: PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen;
489: PetscInt *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs;
490: PetscBool roworiented = a->roworiented;
491: const PetscScalar *value = v;
492: MatScalar *ap, *aa = a->a, *bap;
494: PetscFunctionBegin;
495: rp = aj + ai[row];
496: ap = aa + bs2 * ai[row];
497: rmax = imax[row];
498: nrow = ailen[row];
499: low = 0;
500: high = nrow;
501: value = v;
502: while (high - low > 7) {
503: t = (low + high) / 2;
504: if (rp[t] > col) high = t;
505: else low = t;
506: }
507: for (i = low; i < high; i++) {
508: if (rp[i] > col) break;
509: if (rp[i] == col) {
510: bap = ap + bs2 * i;
511: if (roworiented) {
512: if (is == ADD_VALUES) {
513: for (ii = 0; ii < bs; ii++) {
514: for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
515: }
516: } else {
517: for (ii = 0; ii < bs; ii++) {
518: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
519: }
520: }
521: } else {
522: if (is == ADD_VALUES) {
523: for (ii = 0; ii < bs; ii++, value += bs) {
524: for (jj = 0; jj < bs; jj++) bap[jj] += value[jj];
525: bap += bs;
526: }
527: } else {
528: for (ii = 0; ii < bs; ii++, value += bs) {
529: for (jj = 0; jj < bs; jj++) bap[jj] = value[jj];
530: bap += bs;
531: }
532: }
533: }
534: goto noinsert2;
535: }
536: }
537: if (nonew == 1) goto noinsert2;
538: 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);
539: MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
540: N = nrow++ - 1;
541: high++;
542: /* shift up all the later entries in this row */
543: PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
544: PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
545: rp[i] = col;
546: bap = ap + bs2 * i;
547: if (roworiented) {
548: for (ii = 0; ii < bs; ii++) {
549: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
550: }
551: } else {
552: for (ii = 0; ii < bs; ii++) {
553: for (jj = 0; jj < bs; jj++) *bap++ = *value++;
554: }
555: }
556: noinsert2:;
557: ailen[row] = nrow;
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: /*
562: This routine could be optimized by removing the need for the block copy below and passing stride information
563: to the above inline routines; similarly in MatSetValuesBlocked_MPIBAIJ()
564: */
565: static PetscErrorCode MatSetValuesBlocked_MPISBAIJ(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const MatScalar v[], InsertMode addv)
566: {
567: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
568: const MatScalar *value;
569: MatScalar *barray = baij->barray;
570: PetscBool roworiented = baij->roworiented, ignore_ltriangular = ((Mat_SeqSBAIJ *)baij->A->data)->ignore_ltriangular;
571: PetscInt i, j, ii, jj, row, col, rstart = baij->rstartbs;
572: PetscInt rend = baij->rendbs, cstart = baij->cstartbs, stepval;
573: PetscInt cend = baij->cendbs, bs = mat->rmap->bs, bs2 = baij->bs2;
575: PetscFunctionBegin;
576: if (!barray) {
577: PetscCall(PetscMalloc1(bs2, &barray));
578: baij->barray = barray;
579: }
581: if (roworiented) {
582: stepval = (n - 1) * bs;
583: } else {
584: stepval = (m - 1) * bs;
585: }
586: for (i = 0; i < m; i++) {
587: if (im[i] < 0) continue;
588: 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);
589: if (im[i] >= rstart && im[i] < rend) {
590: row = im[i] - rstart;
591: for (j = 0; j < n; j++) {
592: if (im[i] > in[j]) {
593: PetscCheck(ignore_ltriangular, PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
594: continue; /* ignore lower triangular blocks */
595: }
596: /* If NumCol = 1 then a copy is not required */
597: if (roworiented && n == 1) {
598: barray = (MatScalar *)v + i * bs2;
599: } else if ((!roworiented) && (m == 1)) {
600: barray = (MatScalar *)v + j * bs2;
601: } else { /* Here a copy is required */
602: if (roworiented) {
603: value = v + i * (stepval + bs) * bs + j * bs;
604: } else {
605: value = v + j * (stepval + bs) * bs + i * bs;
606: }
607: for (ii = 0; ii < bs; ii++, value += stepval) {
608: for (jj = 0; jj < bs; jj++) *barray++ = *value++;
609: }
610: barray -= bs2;
611: }
613: if (in[j] >= cstart && in[j] < cend) {
614: col = in[j] - cstart;
615: PetscCall(MatSetValuesBlocked_SeqSBAIJ_Inlined(baij->A, row, col, barray, addv, im[i], in[j]));
616: } else if (in[j] < 0) {
617: continue;
618: } else {
619: 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);
620: if (mat->was_assembled) {
621: if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
623: #if defined(PETSC_USE_CTABLE)
624: PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] + 1, 0, &col));
625: col = col < 1 ? -1 : (col - 1) / bs;
626: #else
627: col = baij->colmap[in[j]] < 1 ? -1 : (baij->colmap[in[j]] - 1) / bs;
628: #endif
629: if (col < 0 && !((Mat_SeqBAIJ *)baij->A->data)->nonew) {
630: PetscCall(MatDisAssemble_MPISBAIJ(mat));
631: col = in[j];
632: }
633: } else col = in[j];
634: PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->B, row, col, barray, addv, im[i], in[j]));
635: }
636: }
637: } else {
638: 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]);
639: if (!baij->donotstash) {
640: if (roworiented) {
641: PetscCall(MatStashValuesRowBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
642: } else {
643: PetscCall(MatStashValuesColBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
644: }
645: }
646: }
647: }
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
651: static PetscErrorCode MatGetValues_MPISBAIJ(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
652: {
653: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
654: PetscInt bs = mat->rmap->bs, i, j, bsrstart = mat->rmap->rstart, bsrend = mat->rmap->rend;
655: PetscInt bscstart = mat->cmap->rstart, bscend = mat->cmap->rend, row, col, data;
657: PetscFunctionBegin;
658: for (i = 0; i < m; i++) {
659: if (idxm[i] < 0) continue; /* negative row */
660: 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);
661: PetscCheck(idxm[i] >= bsrstart && idxm[i] < bsrend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local values currently supported");
662: row = idxm[i] - bsrstart;
663: for (j = 0; j < n; j++) {
664: if (idxn[j] < 0) continue; /* negative column */
665: 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);
666: if (idxn[j] >= bscstart && idxn[j] < bscend) {
667: col = idxn[j] - bscstart;
668: PetscCall(MatGetValues_SeqSBAIJ(baij->A, 1, &row, 1, &col, v + i * n + j));
669: } else {
670: if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
671: #if defined(PETSC_USE_CTABLE)
672: PetscCall(PetscHMapIGetWithDefault(baij->colmap, idxn[j] / bs + 1, 0, &data));
673: data--;
674: #else
675: data = baij->colmap[idxn[j] / bs] - 1;
676: #endif
677: if (data < 0 || baij->garray[data / bs] != idxn[j] / bs) *(v + i * n + j) = 0.0;
678: else {
679: col = data + idxn[j] % bs;
680: PetscCall(MatGetValues_SeqBAIJ(baij->B, 1, &row, 1, &col, v + i * n + j));
681: }
682: }
683: }
684: }
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: static PetscErrorCode MatNorm_MPISBAIJ(Mat mat, NormType type, PetscReal *norm)
689: {
690: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
691: PetscReal sum[2], *lnorm2;
693: PetscFunctionBegin;
694: if (baij->size == 1) {
695: PetscCall(MatNorm(baij->A, type, norm));
696: } else {
697: if (type == NORM_FROBENIUS) {
698: PetscCall(PetscMalloc1(2, &lnorm2));
699: PetscCall(MatNorm(baij->A, type, lnorm2));
700: *lnorm2 = (*lnorm2) * (*lnorm2);
701: lnorm2++; /* square power of norm(A) */
702: PetscCall(MatNorm(baij->B, type, lnorm2));
703: *lnorm2 = (*lnorm2) * (*lnorm2);
704: lnorm2--; /* square power of norm(B) */
705: PetscCallMPI(MPIU_Allreduce(lnorm2, sum, 2, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)mat)));
706: *norm = PetscSqrtReal(sum[0] + 2 * sum[1]);
707: PetscCall(PetscFree(lnorm2));
708: } else if (type == NORM_INFINITY || type == NORM_1) { /* max row/column sum */
709: Mat_SeqSBAIJ *amat = (Mat_SeqSBAIJ *)baij->A->data;
710: Mat_SeqBAIJ *bmat = (Mat_SeqBAIJ *)baij->B->data;
711: PetscReal *rsum, vabs;
712: PetscInt *jj, *garray = baij->garray, rstart = baij->rstartbs, nz;
713: PetscInt brow, bcol, col, bs = baij->A->rmap->bs, row, grow, gcol, mbs = amat->mbs;
714: MatScalar *v;
716: PetscCall(PetscCalloc1(mat->cmap->N, &rsum));
717: /* Amat */
718: v = amat->a;
719: jj = amat->j;
720: for (brow = 0; brow < mbs; brow++) {
721: grow = bs * (rstart + brow);
722: nz = amat->i[brow + 1] - amat->i[brow];
723: for (bcol = 0; bcol < nz; bcol++) {
724: gcol = bs * (rstart + *jj);
725: jj++;
726: for (col = 0; col < bs; col++) {
727: for (row = 0; row < bs; row++) {
728: vabs = PetscAbsScalar(*v);
729: v++;
730: rsum[gcol + col] += vabs;
731: /* non-diagonal block */
732: if (bcol > 0 && vabs > 0.0) rsum[grow + row] += vabs;
733: }
734: }
735: }
736: PetscCall(PetscLogFlops(nz * bs * bs));
737: }
738: /* Bmat */
739: v = bmat->a;
740: jj = bmat->j;
741: for (brow = 0; brow < mbs; brow++) {
742: grow = bs * (rstart + brow);
743: nz = bmat->i[brow + 1] - bmat->i[brow];
744: for (bcol = 0; bcol < nz; bcol++) {
745: gcol = bs * garray[*jj];
746: jj++;
747: for (col = 0; col < bs; col++) {
748: for (row = 0; row < bs; row++) {
749: vabs = PetscAbsScalar(*v);
750: v++;
751: rsum[gcol + col] += vabs;
752: rsum[grow + row] += vabs;
753: }
754: }
755: }
756: PetscCall(PetscLogFlops(nz * bs * bs));
757: }
758: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, rsum, mat->cmap->N, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)mat)));
759: *norm = 0.0;
760: for (col = 0; col < mat->cmap->N; col++) {
761: if (rsum[col] > *norm) *norm = rsum[col];
762: }
763: PetscCall(PetscFree(rsum));
764: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for this norm yet");
765: }
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: static PetscErrorCode MatAssemblyBegin_MPISBAIJ(Mat mat, MatAssemblyType mode)
770: {
771: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
772: PetscInt nstash, reallocs;
774: PetscFunctionBegin;
775: if (baij->donotstash || mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);
777: PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
778: PetscCall(MatStashScatterBegin_Private(mat, &mat->bstash, baij->rangebs));
779: PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs));
780: PetscCall(PetscInfo(mat, "Stash has %" PetscInt_FMT " entries,uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
781: PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs));
782: PetscCall(PetscInfo(mat, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }
786: static PetscErrorCode MatAssemblyEnd_MPISBAIJ(Mat mat, MatAssemblyType mode)
787: {
788: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
789: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)baij->A->data;
790: PetscInt i, j, rstart, ncols, flg, bs2 = baij->bs2;
791: PetscInt *row, *col;
792: PetscBool all_assembled;
793: PetscMPIInt n;
794: PetscBool r1, r2, r3;
795: MatScalar *val;
797: /* do not use 'b=(Mat_SeqBAIJ*)baij->B->data' as B can be reset in disassembly */
798: PetscFunctionBegin;
799: if (!baij->donotstash && !mat->nooffprocentries) {
800: while (1) {
801: PetscCall(MatStashScatterGetMesg_Private(&mat->stash, &n, &row, &col, &val, &flg));
802: if (!flg) break;
804: for (i = 0; i < n;) {
805: /* Now identify the consecutive vals belonging to the same row */
806: for (j = i, rstart = row[j]; j < n; j++) {
807: if (row[j] != rstart) break;
808: }
809: if (j < n) ncols = j - i;
810: else ncols = n - i;
811: /* Now assemble all these values with a single function call */
812: PetscCall(MatSetValues_MPISBAIJ(mat, 1, row + i, ncols, col + i, val + i, mat->insertmode));
813: i = j;
814: }
815: }
816: PetscCall(MatStashScatterEnd_Private(&mat->stash));
817: /* Now process the block-stash. Since the values are stashed column-oriented,
818: set the row-oriented flag to column-oriented, and after MatSetValues()
819: restore the original flags */
820: r1 = baij->roworiented;
821: r2 = a->roworiented;
822: r3 = ((Mat_SeqBAIJ *)baij->B->data)->roworiented;
824: baij->roworiented = PETSC_FALSE;
825: a->roworiented = PETSC_FALSE;
827: ((Mat_SeqBAIJ *)baij->B->data)->roworiented = PETSC_FALSE; /* b->roworiented */
828: while (1) {
829: PetscCall(MatStashScatterGetMesg_Private(&mat->bstash, &n, &row, &col, &val, &flg));
830: if (!flg) break;
832: for (i = 0; i < n;) {
833: /* Now identify the consecutive vals belonging to the same row */
834: for (j = i, rstart = row[j]; j < n; j++) {
835: if (row[j] != rstart) break;
836: }
837: if (j < n) ncols = j - i;
838: else ncols = n - i;
839: PetscCall(MatSetValuesBlocked_MPISBAIJ(mat, 1, row + i, ncols, col + i, val + i * bs2, mat->insertmode));
840: i = j;
841: }
842: }
843: PetscCall(MatStashScatterEnd_Private(&mat->bstash));
845: baij->roworiented = r1;
846: a->roworiented = r2;
848: ((Mat_SeqBAIJ *)baij->B->data)->roworiented = r3; /* b->roworiented */
849: }
851: PetscCall(MatAssemblyBegin(baij->A, mode));
852: PetscCall(MatAssemblyEnd(baij->A, mode));
854: /* determine if any process has disassembled, if so we must
855: also disassemble ourselves, in order that we may reassemble. */
856: /*
857: if nonzero structure of submatrix B cannot change then we know that
858: no process disassembled thus we can skip this stuff
859: */
860: if (!((Mat_SeqBAIJ *)baij->B->data)->nonew) {
861: PetscCallMPI(MPIU_Allreduce(&mat->was_assembled, &all_assembled, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
862: if (mat->was_assembled && !all_assembled) PetscCall(MatDisAssemble_MPISBAIJ(mat));
863: }
865: if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) PetscCall(MatSetUpMultiply_MPISBAIJ(mat)); /* setup Mvctx and sMvctx */
866: PetscCall(MatAssemblyBegin(baij->B, mode));
867: PetscCall(MatAssemblyEnd(baij->B, mode));
869: PetscCall(PetscFree2(baij->rowvalues, baij->rowindices));
871: baij->rowvalues = NULL;
873: /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
874: if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqBAIJ *)baij->A->data)->nonew) {
875: PetscObjectState state = baij->A->nonzerostate + baij->B->nonzerostate;
876: PetscCallMPI(MPIU_Allreduce(&state, &mat->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)mat)));
877: }
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: extern PetscErrorCode MatSetValues_MPIBAIJ(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode);
882: #include <petscdraw.h>
883: static PetscErrorCode MatView_MPISBAIJ_ASCIIorDraworSocket(Mat mat, PetscViewer viewer)
884: {
885: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
886: PetscInt bs = mat->rmap->bs;
887: PetscMPIInt rank = baij->rank;
888: PetscBool isascii, isdraw;
889: PetscViewer sviewer;
890: PetscViewerFormat format;
892: PetscFunctionBegin;
893: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
894: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
895: if (isascii) {
896: PetscCall(PetscViewerGetFormat(viewer, &format));
897: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
898: MatInfo info;
899: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank));
900: PetscCall(MatGetInfo(mat, MAT_LOCAL, &info));
901: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
902: 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,
903: mat->rmap->bs, info.memory));
904: PetscCall(MatGetInfo(baij->A, MAT_LOCAL, &info));
905: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] on-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
906: PetscCall(MatGetInfo(baij->B, MAT_LOCAL, &info));
907: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] off-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
908: PetscCall(PetscViewerFlush(viewer));
909: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
910: PetscCall(PetscViewerASCIIPrintf(viewer, "Information on VecScatter used in matrix-vector product: \n"));
911: PetscCall(VecScatterView(baij->Mvctx, viewer));
912: PetscFunctionReturn(PETSC_SUCCESS);
913: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_FACTOR_INFO) PetscFunctionReturn(PETSC_SUCCESS);
914: }
916: if (isdraw) {
917: PetscDraw draw;
918: PetscBool isnull;
919: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
920: PetscCall(PetscDrawIsNull(draw, &isnull));
921: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
922: }
924: {
925: /* assemble the entire matrix onto first processor. */
926: Mat A;
927: Mat_SeqSBAIJ *Aloc;
928: Mat_SeqBAIJ *Bloc;
929: PetscInt M = mat->rmap->N, N = mat->cmap->N, *ai, *aj, col, i, j, k, *rvals, mbs = baij->mbs;
930: MatScalar *a;
931: const char *matname;
933: /* Should this be the same type as mat? */
934: PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &A));
935: if (rank == 0) {
936: PetscCall(MatSetSizes(A, M, N, M, N));
937: } else {
938: PetscCall(MatSetSizes(A, 0, 0, M, N));
939: }
940: PetscCall(MatSetType(A, MATMPISBAIJ));
941: PetscCall(MatMPISBAIJSetPreallocation(A, mat->rmap->bs, 0, NULL, 0, NULL));
942: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
944: /* copy over the A part */
945: Aloc = (Mat_SeqSBAIJ *)baij->A->data;
946: ai = Aloc->i;
947: aj = Aloc->j;
948: a = Aloc->a;
949: PetscCall(PetscMalloc1(bs, &rvals));
951: for (i = 0; i < mbs; i++) {
952: rvals[0] = bs * (baij->rstartbs + i);
953: for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
954: for (j = ai[i]; j < ai[i + 1]; j++) {
955: col = (baij->cstartbs + aj[j]) * bs;
956: for (k = 0; k < bs; k++) {
957: PetscCall(MatSetValues_MPISBAIJ(A, bs, rvals, 1, &col, a, INSERT_VALUES));
958: col++;
959: a += bs;
960: }
961: }
962: }
963: /* copy over the B part */
964: Bloc = (Mat_SeqBAIJ *)baij->B->data;
965: ai = Bloc->i;
966: aj = Bloc->j;
967: a = Bloc->a;
968: for (i = 0; i < mbs; i++) {
969: rvals[0] = bs * (baij->rstartbs + i);
970: for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
971: for (j = ai[i]; j < ai[i + 1]; j++) {
972: col = baij->garray[aj[j]] * bs;
973: for (k = 0; k < bs; k++) {
974: PetscCall(MatSetValues_MPIBAIJ(A, bs, rvals, 1, &col, a, INSERT_VALUES));
975: col++;
976: a += bs;
977: }
978: }
979: }
980: PetscCall(PetscFree(rvals));
981: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
982: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
983: /*
984: Everyone has to call to draw the matrix since the graphics waits are
985: synchronized across all processors that share the PetscDraw object
986: */
987: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
988: if (((PetscObject)mat)->name) PetscCall(PetscObjectGetName((PetscObject)mat, &matname));
989: if (rank == 0) {
990: if (((PetscObject)mat)->name) PetscCall(PetscObjectSetName((PetscObject)((Mat_MPISBAIJ *)A->data)->A, matname));
991: PetscCall(MatView_SeqSBAIJ(((Mat_MPISBAIJ *)A->data)->A, sviewer));
992: }
993: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
994: PetscCall(MatDestroy(&A));
995: }
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: /* Used for both MPIBAIJ and MPISBAIJ matrices */
1000: #define MatView_MPISBAIJ_Binary MatView_MPIBAIJ_Binary
1002: static PetscErrorCode MatView_MPISBAIJ(Mat mat, PetscViewer viewer)
1003: {
1004: PetscBool isascii, isdraw, issocket, isbinary;
1006: PetscFunctionBegin;
1007: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1008: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1009: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
1010: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1011: if (isascii || isdraw || issocket) PetscCall(MatView_MPISBAIJ_ASCIIorDraworSocket(mat, viewer));
1012: else if (isbinary) PetscCall(MatView_MPISBAIJ_Binary(mat, viewer));
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: #if defined(PETSC_USE_COMPLEX)
1017: static PetscErrorCode MatMult_MPISBAIJ_Hermitian(Mat A, Vec xx, Vec yy)
1018: {
1019: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1020: PetscInt mbs = a->mbs, bs = A->rmap->bs;
1021: PetscScalar *from;
1022: const PetscScalar *x;
1024: PetscFunctionBegin;
1025: /* diagonal part */
1026: PetscCall((*a->A->ops->mult)(a->A, xx, a->slvec1a));
1027: /* since a->slvec1b shares memory (dangerously) with a->slec1 changes to a->slec1 will affect it */
1028: PetscCall(PetscObjectStateIncrease((PetscObject)a->slvec1b));
1029: PetscCall(VecZeroEntries(a->slvec1b));
1031: /* subdiagonal part */
1032: PetscCheck(a->B->ops->multhermitiantranspose, PetscObjectComm((PetscObject)a->B), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)a->B)->type_name);
1033: PetscCall((*a->B->ops->multhermitiantranspose)(a->B, xx, a->slvec0b));
1035: /* copy x into the vec slvec0 */
1036: PetscCall(VecGetArray(a->slvec0, &from));
1037: PetscCall(VecGetArrayRead(xx, &x));
1039: PetscCall(PetscArraycpy(from, x, bs * mbs));
1040: PetscCall(VecRestoreArray(a->slvec0, &from));
1041: PetscCall(VecRestoreArrayRead(xx, &x));
1043: PetscCall(VecScatterBegin(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
1044: PetscCall(VecScatterEnd(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
1045: /* supperdiagonal part */
1046: PetscCall((*a->B->ops->multadd)(a->B, a->slvec1b, a->slvec1a, yy));
1047: PetscFunctionReturn(PETSC_SUCCESS);
1048: }
1049: #endif
1051: static PetscErrorCode MatMult_MPISBAIJ(Mat A, Vec xx, Vec yy)
1052: {
1053: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1054: PetscInt mbs = a->mbs, bs = A->rmap->bs;
1055: PetscScalar *from;
1056: const PetscScalar *x;
1058: PetscFunctionBegin;
1059: /* diagonal part */
1060: PetscCall((*a->A->ops->mult)(a->A, xx, a->slvec1a));
1061: /* since a->slvec1b shares memory (dangerously) with a->slec1 changes to a->slec1 will affect it */
1062: PetscCall(PetscObjectStateIncrease((PetscObject)a->slvec1b));
1063: PetscCall(VecZeroEntries(a->slvec1b));
1065: /* subdiagonal part */
1066: PetscCall((*a->B->ops->multtranspose)(a->B, xx, a->slvec0b));
1068: /* copy x into the vec slvec0 */
1069: PetscCall(VecGetArray(a->slvec0, &from));
1070: PetscCall(VecGetArrayRead(xx, &x));
1072: PetscCall(PetscArraycpy(from, x, bs * mbs));
1073: PetscCall(VecRestoreArray(a->slvec0, &from));
1074: PetscCall(VecRestoreArrayRead(xx, &x));
1076: PetscCall(VecScatterBegin(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
1077: PetscCall(VecScatterEnd(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
1078: /* supperdiagonal part */
1079: PetscCall((*a->B->ops->multadd)(a->B, a->slvec1b, a->slvec1a, yy));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: #if PetscDefined(USE_COMPLEX)
1084: static PetscErrorCode MatMultAdd_MPISBAIJ_Hermitian(Mat A, Vec xx, Vec yy, Vec zz)
1085: {
1086: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1087: PetscInt mbs = a->mbs, bs = A->rmap->bs;
1088: PetscScalar *from;
1089: const PetscScalar *x;
1091: PetscFunctionBegin;
1092: /* diagonal part */
1093: PetscCall((*a->A->ops->multadd)(a->A, xx, yy, a->slvec1a));
1094: PetscCall(PetscObjectStateIncrease((PetscObject)a->slvec1b));
1095: PetscCall(VecZeroEntries(a->slvec1b));
1097: /* subdiagonal part */
1098: PetscCheck(a->B->ops->multhermitiantranspose, PetscObjectComm((PetscObject)a->B), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)a->B)->type_name);
1099: PetscCall((*a->B->ops->multhermitiantranspose)(a->B, xx, a->slvec0b));
1101: /* copy x into the vec slvec0 */
1102: PetscCall(VecGetArray(a->slvec0, &from));
1103: PetscCall(VecGetArrayRead(xx, &x));
1104: PetscCall(PetscArraycpy(from, x, bs * mbs));
1105: PetscCall(VecRestoreArray(a->slvec0, &from));
1107: PetscCall(VecScatterBegin(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
1108: PetscCall(VecRestoreArrayRead(xx, &x));
1109: PetscCall(VecScatterEnd(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
1111: /* supperdiagonal part */
1112: PetscCall((*a->B->ops->multadd)(a->B, a->slvec1b, a->slvec1a, zz));
1113: PetscFunctionReturn(PETSC_SUCCESS);
1114: }
1115: #endif
1117: static PetscErrorCode MatMultAdd_MPISBAIJ(Mat A, Vec xx, Vec yy, Vec zz)
1118: {
1119: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1120: PetscInt mbs = a->mbs, bs = A->rmap->bs;
1121: PetscScalar *from;
1122: const PetscScalar *x;
1124: PetscFunctionBegin;
1125: /* diagonal part */
1126: PetscCall((*a->A->ops->multadd)(a->A, xx, yy, a->slvec1a));
1127: PetscCall(PetscObjectStateIncrease((PetscObject)a->slvec1b));
1128: PetscCall(VecZeroEntries(a->slvec1b));
1130: /* subdiagonal part */
1131: PetscCall((*a->B->ops->multtranspose)(a->B, xx, a->slvec0b));
1133: /* copy x into the vec slvec0 */
1134: PetscCall(VecGetArray(a->slvec0, &from));
1135: PetscCall(VecGetArrayRead(xx, &x));
1136: PetscCall(PetscArraycpy(from, x, bs * mbs));
1137: PetscCall(VecRestoreArray(a->slvec0, &from));
1139: PetscCall(VecScatterBegin(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
1140: PetscCall(VecRestoreArrayRead(xx, &x));
1141: PetscCall(VecScatterEnd(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
1143: /* supperdiagonal part */
1144: PetscCall((*a->B->ops->multadd)(a->B, a->slvec1b, a->slvec1a, zz));
1145: PetscFunctionReturn(PETSC_SUCCESS);
1146: }
1148: /*
1149: This only works correctly for square matrices where the subblock A->A is the
1150: diagonal block
1151: */
1152: static PetscErrorCode MatGetDiagonal_MPISBAIJ(Mat A, Vec v)
1153: {
1154: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1156: PetscFunctionBegin;
1157: /* PetscCheck(a->rmap->N == a->cmap->N,PETSC_COMM_SELF,PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block"); */
1158: PetscCall(MatGetDiagonal(a->A, v));
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: static PetscErrorCode MatScale_MPISBAIJ(Mat A, PetscScalar aa)
1163: {
1164: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1166: PetscFunctionBegin;
1167: PetscCall(MatScale(a->A, aa));
1168: PetscCall(MatScale(a->B, aa));
1169: PetscFunctionReturn(PETSC_SUCCESS);
1170: }
1172: static PetscErrorCode MatGetRow_MPISBAIJ(Mat matin, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1173: {
1174: Mat_MPISBAIJ *mat = (Mat_MPISBAIJ *)matin->data;
1175: PetscScalar *vworkA, *vworkB, **pvA, **pvB, *v_p;
1176: PetscInt bs = matin->rmap->bs, bs2 = mat->bs2, i, *cworkA, *cworkB, **pcA, **pcB;
1177: PetscInt nztot, nzA, nzB, lrow, brstart = matin->rmap->rstart, brend = matin->rmap->rend;
1178: PetscInt *cmap, *idx_p, cstart = mat->rstartbs;
1180: PetscFunctionBegin;
1181: PetscCheck(!mat->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Already active");
1182: mat->getrowactive = PETSC_TRUE;
1184: if (!mat->rowvalues && (idx || v)) {
1185: /*
1186: allocate enough space to hold information from the longest row.
1187: */
1188: Mat_SeqSBAIJ *Aa = (Mat_SeqSBAIJ *)mat->A->data;
1189: Mat_SeqBAIJ *Ba = (Mat_SeqBAIJ *)mat->B->data;
1190: PetscInt max = 1, mbs = mat->mbs, tmp;
1191: for (i = 0; i < mbs; i++) {
1192: tmp = Aa->i[i + 1] - Aa->i[i] + Ba->i[i + 1] - Ba->i[i]; /* row length */
1193: if (max < tmp) max = tmp;
1194: }
1195: PetscCall(PetscMalloc2(max * bs2, &mat->rowvalues, max * bs2, &mat->rowindices));
1196: }
1198: PetscCheck(row >= brstart && row < brend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local rows");
1199: lrow = row - brstart; /* local row index */
1201: pvA = &vworkA;
1202: pcA = &cworkA;
1203: pvB = &vworkB;
1204: pcB = &cworkB;
1205: if (!v) {
1206: pvA = NULL;
1207: pvB = NULL;
1208: }
1209: if (!idx) {
1210: pcA = NULL;
1211: if (!v) pcB = NULL;
1212: }
1213: PetscCall((*mat->A->ops->getrow)(mat->A, lrow, &nzA, pcA, pvA));
1214: PetscCall((*mat->B->ops->getrow)(mat->B, lrow, &nzB, pcB, pvB));
1215: nztot = nzA + nzB;
1217: cmap = mat->garray;
1218: if (v || idx) {
1219: if (nztot) {
1220: /* Sort by increasing column numbers, assuming A and B already sorted */
1221: PetscInt imark = -1;
1222: if (v) {
1223: *v = v_p = mat->rowvalues;
1224: for (i = 0; i < nzB; i++) {
1225: if (cmap[cworkB[i] / bs] < cstart) v_p[i] = vworkB[i];
1226: else break;
1227: }
1228: imark = i;
1229: for (i = 0; i < nzA; i++) v_p[imark + i] = vworkA[i];
1230: for (i = imark; i < nzB; i++) v_p[nzA + i] = vworkB[i];
1231: }
1232: if (idx) {
1233: *idx = idx_p = mat->rowindices;
1234: if (imark > -1) {
1235: for (i = 0; i < imark; i++) idx_p[i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1236: } else {
1237: for (i = 0; i < nzB; i++) {
1238: if (cmap[cworkB[i] / bs] < cstart) idx_p[i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1239: else break;
1240: }
1241: imark = i;
1242: }
1243: for (i = 0; i < nzA; i++) idx_p[imark + i] = cstart * bs + cworkA[i];
1244: for (i = imark; i < nzB; i++) idx_p[nzA + i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1245: }
1246: } else {
1247: if (idx) *idx = NULL;
1248: if (v) *v = NULL;
1249: }
1250: }
1251: *nz = nztot;
1252: PetscCall((*mat->A->ops->restorerow)(mat->A, lrow, &nzA, pcA, pvA));
1253: PetscCall((*mat->B->ops->restorerow)(mat->B, lrow, &nzB, pcB, pvB));
1254: PetscFunctionReturn(PETSC_SUCCESS);
1255: }
1257: static PetscErrorCode MatRestoreRow_MPISBAIJ(Mat mat, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1258: {
1259: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
1261: PetscFunctionBegin;
1262: PetscCheck(baij->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "MatGetRow() must be called first");
1263: baij->getrowactive = PETSC_FALSE;
1264: PetscFunctionReturn(PETSC_SUCCESS);
1265: }
1267: static PetscErrorCode MatGetRowUpperTriangular_MPISBAIJ(Mat A)
1268: {
1269: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1270: Mat_SeqSBAIJ *aA = (Mat_SeqSBAIJ *)a->A->data;
1272: PetscFunctionBegin;
1273: aA->getrow_utriangular = PETSC_TRUE;
1274: PetscFunctionReturn(PETSC_SUCCESS);
1275: }
1276: static PetscErrorCode MatRestoreRowUpperTriangular_MPISBAIJ(Mat A)
1277: {
1278: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1279: Mat_SeqSBAIJ *aA = (Mat_SeqSBAIJ *)a->A->data;
1281: PetscFunctionBegin;
1282: aA->getrow_utriangular = PETSC_FALSE;
1283: PetscFunctionReturn(PETSC_SUCCESS);
1284: }
1286: static PetscErrorCode MatConjugate_MPISBAIJ(Mat mat)
1287: {
1288: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)mat->data;
1290: PetscFunctionBegin;
1291: PetscCall(MatConjugate(a->A));
1292: PetscCall(MatConjugate(a->B));
1293: PetscFunctionReturn(PETSC_SUCCESS);
1294: }
1296: static PetscErrorCode MatRealPart_MPISBAIJ(Mat A)
1297: {
1298: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1300: PetscFunctionBegin;
1301: PetscCall(MatRealPart(a->A));
1302: PetscCall(MatRealPart(a->B));
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: static PetscErrorCode MatImaginaryPart_MPISBAIJ(Mat A)
1307: {
1308: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1310: PetscFunctionBegin;
1311: PetscCall(MatImaginaryPart(a->A));
1312: PetscCall(MatImaginaryPart(a->B));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: /* Check if isrow is a subset of iscol_local, called by MatCreateSubMatrix_MPISBAIJ()
1317: Input: isrow - distributed(parallel),
1318: iscol_local - locally owned (seq)
1319: */
1320: static PetscErrorCode ISEqual_private(IS isrow, IS iscol_local, PetscBool *flg)
1321: {
1322: PetscInt sz1, sz2, *a1, *a2, i, j, k, nmatch;
1323: const PetscInt *ptr1, *ptr2;
1325: PetscFunctionBegin;
1326: *flg = PETSC_FALSE;
1327: PetscCall(ISGetLocalSize(isrow, &sz1));
1328: PetscCall(ISGetLocalSize(iscol_local, &sz2));
1329: if (sz1 > sz2) PetscFunctionReturn(PETSC_SUCCESS);
1331: PetscCall(ISGetIndices(isrow, &ptr1));
1332: PetscCall(ISGetIndices(iscol_local, &ptr2));
1334: PetscCall(PetscMalloc1(sz1, &a1));
1335: PetscCall(PetscMalloc1(sz2, &a2));
1336: PetscCall(PetscArraycpy(a1, ptr1, sz1));
1337: PetscCall(PetscArraycpy(a2, ptr2, sz2));
1338: PetscCall(PetscSortInt(sz1, a1));
1339: PetscCall(PetscSortInt(sz2, a2));
1341: nmatch = 0;
1342: k = 0;
1343: for (i = 0; i < sz1; i++) {
1344: for (j = k; j < sz2; j++) {
1345: if (a1[i] == a2[j]) {
1346: k = j;
1347: nmatch++;
1348: break;
1349: }
1350: }
1351: }
1352: PetscCall(ISRestoreIndices(isrow, &ptr1));
1353: PetscCall(ISRestoreIndices(iscol_local, &ptr2));
1354: PetscCall(PetscFree(a1));
1355: PetscCall(PetscFree(a2));
1356: if (nmatch < sz1) {
1357: *flg = PETSC_FALSE;
1358: } else {
1359: *flg = PETSC_TRUE;
1360: }
1361: PetscFunctionReturn(PETSC_SUCCESS);
1362: }
1364: static PetscErrorCode MatCreateSubMatrix_MPISBAIJ(Mat mat, IS isrow, IS iscol, MatReuse call, Mat *newmat)
1365: {
1366: Mat C[2];
1367: IS iscol_local, isrow_local;
1368: PetscInt csize, csize_local, rsize;
1369: PetscBool isequal, issorted, isidentity = PETSC_FALSE;
1371: PetscFunctionBegin;
1372: PetscCall(ISGetLocalSize(iscol, &csize));
1373: PetscCall(ISGetLocalSize(isrow, &rsize));
1374: if (call == MAT_REUSE_MATRIX) {
1375: PetscCall(PetscObjectQuery((PetscObject)*newmat, "ISAllGather", (PetscObject *)&iscol_local));
1376: PetscCheck(iscol_local, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Submatrix passed in was not used before, cannot reuse");
1377: } else {
1378: PetscCall(ISAllGather(iscol, &iscol_local));
1379: PetscCall(ISSorted(iscol_local, &issorted));
1380: PetscCheck(issorted, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "For symmetric format, iscol must be sorted");
1381: }
1382: PetscCall(ISEqual_private(isrow, iscol_local, &isequal));
1383: if (!isequal) {
1384: PetscCall(ISGetLocalSize(iscol_local, &csize_local));
1385: isidentity = (PetscBool)(mat->cmap->N == csize_local);
1386: if (!isidentity) {
1387: if (call == MAT_REUSE_MATRIX) {
1388: PetscCall(PetscObjectQuery((PetscObject)*newmat, "ISAllGather_other", (PetscObject *)&isrow_local));
1389: PetscCheck(isrow_local, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Submatrix passed in was not used before, cannot reuse");
1390: } else {
1391: PetscCall(ISAllGather(isrow, &isrow_local));
1392: PetscCall(ISSorted(isrow_local, &issorted));
1393: PetscCheck(issorted, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "For symmetric format, isrow must be sorted");
1394: }
1395: }
1396: }
1397: /* now call MatCreateSubMatrix_MPIBAIJ() */
1398: PetscCall(MatCreateSubMatrix_MPIBAIJ_Private(mat, isrow, iscol_local, csize, isequal || isidentity ? call : MAT_INITIAL_MATRIX, isequal || isidentity ? newmat : C, (PetscBool)(isequal || isidentity)));
1399: if (!isequal && !isidentity) {
1400: if (call == MAT_INITIAL_MATRIX) {
1401: IS intersect;
1402: PetscInt ni;
1404: PetscCall(ISIntersect(isrow_local, iscol_local, &intersect));
1405: PetscCall(ISGetLocalSize(intersect, &ni));
1406: PetscCall(ISDestroy(&intersect));
1407: PetscCheck(ni == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot create such a submatrix: for symmetric format, when requesting an off-diagonal submatrix, isrow and iscol should have an empty intersection (number of common indices is %" PetscInt_FMT ")", ni);
1408: }
1409: PetscCall(MatCreateSubMatrix_MPIBAIJ_Private(mat, iscol, isrow_local, rsize, MAT_INITIAL_MATRIX, C + 1, PETSC_FALSE));
1410: PetscCall(MatTranspose(C[1], MAT_INPLACE_MATRIX, C + 1));
1411: PetscCall(MatAXPY(C[0], 1.0, C[1], DIFFERENT_NONZERO_PATTERN));
1412: if (call == MAT_REUSE_MATRIX) PetscCall(MatCopy(C[0], *newmat, SAME_NONZERO_PATTERN));
1413: else if (mat->rmap->bs == 1) PetscCall(MatConvert(C[0], MATAIJ, MAT_INITIAL_MATRIX, newmat));
1414: else PetscCall(MatCopy(C[0], *newmat, SAME_NONZERO_PATTERN));
1415: PetscCall(MatDestroy(C));
1416: PetscCall(MatDestroy(C + 1));
1417: }
1418: if (call == MAT_INITIAL_MATRIX) {
1419: if (!isequal && !isidentity) {
1420: PetscCall(PetscObjectCompose((PetscObject)*newmat, "ISAllGather_other", (PetscObject)isrow_local));
1421: PetscCall(ISDestroy(&isrow_local));
1422: }
1423: PetscCall(PetscObjectCompose((PetscObject)*newmat, "ISAllGather", (PetscObject)iscol_local));
1424: PetscCall(ISDestroy(&iscol_local));
1425: }
1426: PetscFunctionReturn(PETSC_SUCCESS);
1427: }
1429: static PetscErrorCode MatZeroEntries_MPISBAIJ(Mat A)
1430: {
1431: Mat_MPISBAIJ *l = (Mat_MPISBAIJ *)A->data;
1433: PetscFunctionBegin;
1434: PetscCall(MatZeroEntries(l->A));
1435: PetscCall(MatZeroEntries(l->B));
1436: PetscFunctionReturn(PETSC_SUCCESS);
1437: }
1439: static PetscErrorCode MatGetInfo_MPISBAIJ(Mat matin, MatInfoType flag, MatInfo *info)
1440: {
1441: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)matin->data;
1442: Mat A = a->A, B = a->B;
1443: PetscLogDouble isend[5], irecv[5];
1445: PetscFunctionBegin;
1446: info->block_size = (PetscReal)matin->rmap->bs;
1448: PetscCall(MatGetInfo(A, MAT_LOCAL, info));
1450: isend[0] = info->nz_used;
1451: isend[1] = info->nz_allocated;
1452: isend[2] = info->nz_unneeded;
1453: isend[3] = info->memory;
1454: isend[4] = info->mallocs;
1456: PetscCall(MatGetInfo(B, MAT_LOCAL, info));
1458: isend[0] += info->nz_used;
1459: isend[1] += info->nz_allocated;
1460: isend[2] += info->nz_unneeded;
1461: isend[3] += info->memory;
1462: isend[4] += info->mallocs;
1463: if (flag == MAT_LOCAL) {
1464: info->nz_used = isend[0];
1465: info->nz_allocated = isend[1];
1466: info->nz_unneeded = isend[2];
1467: info->memory = isend[3];
1468: info->mallocs = isend[4];
1469: } else if (flag == MAT_GLOBAL_MAX) {
1470: PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)matin)));
1472: info->nz_used = irecv[0];
1473: info->nz_allocated = irecv[1];
1474: info->nz_unneeded = irecv[2];
1475: info->memory = irecv[3];
1476: info->mallocs = irecv[4];
1477: } else if (flag == MAT_GLOBAL_SUM) {
1478: PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)matin)));
1480: info->nz_used = irecv[0];
1481: info->nz_allocated = irecv[1];
1482: info->nz_unneeded = irecv[2];
1483: info->memory = irecv[3];
1484: info->mallocs = irecv[4];
1485: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatInfoType argument %d", (int)flag);
1486: info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */
1487: info->fill_ratio_needed = 0;
1488: info->factor_mallocs = 0;
1489: PetscFunctionReturn(PETSC_SUCCESS);
1490: }
1492: static PetscErrorCode MatSetOption_MPISBAIJ(Mat A, MatOption op, PetscBool flg)
1493: {
1494: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1495: Mat_SeqSBAIJ *aA = (Mat_SeqSBAIJ *)a->A->data;
1497: PetscFunctionBegin;
1498: switch (op) {
1499: case MAT_NEW_NONZERO_LOCATIONS:
1500: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1501: case MAT_UNUSED_NONZERO_LOCATION_ERR:
1502: case MAT_KEEP_NONZERO_PATTERN:
1503: case MAT_NEW_NONZERO_LOCATION_ERR:
1504: MatCheckPreallocated(A, 1);
1505: PetscCall(MatSetOption(a->A, op, flg));
1506: PetscCall(MatSetOption(a->B, op, flg));
1507: break;
1508: case MAT_ROW_ORIENTED:
1509: MatCheckPreallocated(A, 1);
1510: a->roworiented = flg;
1512: PetscCall(MatSetOption(a->A, op, flg));
1513: PetscCall(MatSetOption(a->B, op, flg));
1514: break;
1515: case MAT_IGNORE_OFF_PROC_ENTRIES:
1516: a->donotstash = flg;
1517: break;
1518: case MAT_USE_HASH_TABLE:
1519: a->ht_flag = flg;
1520: break;
1521: case MAT_HERMITIAN:
1522: if (a->A && A->rmap->n == A->cmap->n) PetscCall(MatSetOption(a->A, op, flg));
1523: #if defined(PETSC_USE_COMPLEX)
1524: if (flg) { /* need different mat-vec ops */
1525: A->ops->mult = MatMult_MPISBAIJ_Hermitian;
1526: A->ops->multadd = MatMultAdd_MPISBAIJ_Hermitian;
1527: A->ops->multtranspose = NULL;
1528: A->ops->multtransposeadd = NULL;
1529: }
1530: #endif
1531: break;
1532: case MAT_SPD:
1533: case MAT_SYMMETRIC:
1534: if (a->A && A->rmap->n == A->cmap->n) PetscCall(MatSetOption(a->A, op, flg));
1535: #if defined(PETSC_USE_COMPLEX)
1536: if (flg) { /* restore to use default mat-vec ops */
1537: A->ops->mult = MatMult_MPISBAIJ;
1538: A->ops->multadd = MatMultAdd_MPISBAIJ;
1539: A->ops->multtranspose = MatMult_MPISBAIJ;
1540: A->ops->multtransposeadd = MatMultAdd_MPISBAIJ;
1541: }
1542: #endif
1543: break;
1544: case MAT_STRUCTURALLY_SYMMETRIC:
1545: if (a->A && A->rmap->n == A->cmap->n) PetscCall(MatSetOption(a->A, op, flg));
1546: break;
1547: case MAT_IGNORE_LOWER_TRIANGULAR:
1548: case MAT_ERROR_LOWER_TRIANGULAR:
1549: aA->ignore_ltriangular = flg;
1550: break;
1551: case MAT_GETROW_UPPERTRIANGULAR:
1552: aA->getrow_utriangular = flg;
1553: break;
1554: default:
1555: break;
1556: }
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: static PetscErrorCode MatTranspose_MPISBAIJ(Mat A, MatReuse reuse, Mat *B)
1561: {
1562: PetscFunctionBegin;
1563: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1564: if (reuse == MAT_INITIAL_MATRIX) {
1565: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, B));
1566: } else if (reuse == MAT_REUSE_MATRIX) {
1567: PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN));
1568: }
1569: PetscFunctionReturn(PETSC_SUCCESS);
1570: }
1572: static PetscErrorCode MatDiagonalScale_MPISBAIJ(Mat mat, Vec ll, Vec rr)
1573: {
1574: Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
1575: Mat a = baij->A, b = baij->B;
1576: PetscInt nv, m, n;
1578: PetscFunctionBegin;
1579: if (!ll) PetscFunctionReturn(PETSC_SUCCESS);
1581: PetscCall(MatGetLocalSize(mat, &m, &n));
1582: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "For symmetric format, local size %" PetscInt_FMT " %" PetscInt_FMT " must be same", m, n);
1584: PetscCall(VecGetLocalSize(rr, &nv));
1585: PetscCheck(nv == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left and right vector non-conforming local size");
1587: PetscCall(VecScatterBegin(baij->Mvctx, rr, baij->lvec, INSERT_VALUES, SCATTER_FORWARD));
1589: /* left diagonalscale the off-diagonal part */
1590: PetscUseTypeMethod(b, diagonalscale, ll, NULL);
1592: /* scale the diagonal part */
1593: PetscUseTypeMethod(a, diagonalscale, ll, rr);
1595: /* right diagonalscale the off-diagonal part */
1596: PetscCall(VecScatterEnd(baij->Mvctx, rr, baij->lvec, INSERT_VALUES, SCATTER_FORWARD));
1597: PetscUseTypeMethod(b, diagonalscale, NULL, baij->lvec);
1598: PetscFunctionReturn(PETSC_SUCCESS);
1599: }
1601: static PetscErrorCode MatSetUnfactored_MPISBAIJ(Mat A)
1602: {
1603: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1605: PetscFunctionBegin;
1606: PetscCall(MatSetUnfactored(a->A));
1607: PetscFunctionReturn(PETSC_SUCCESS);
1608: }
1610: static PetscErrorCode MatDuplicate_MPISBAIJ(Mat, MatDuplicateOption, Mat *);
1612: static PetscErrorCode MatEqual_MPISBAIJ(Mat A, Mat B, PetscBool *flag)
1613: {
1614: Mat_MPISBAIJ *matB = (Mat_MPISBAIJ *)B->data, *matA = (Mat_MPISBAIJ *)A->data;
1615: Mat a, b, c, d;
1616: PetscBool flg;
1618: PetscFunctionBegin;
1619: a = matA->A;
1620: b = matA->B;
1621: c = matB->A;
1622: d = matB->B;
1624: PetscCall(MatEqual(a, c, &flg));
1625: if (flg) PetscCall(MatEqual(b, d, &flg));
1626: PetscCallMPI(MPIU_Allreduce(&flg, flag, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
1627: PetscFunctionReturn(PETSC_SUCCESS);
1628: }
1630: static PetscErrorCode MatCopy_MPISBAIJ(Mat A, Mat B, MatStructure str)
1631: {
1632: PetscBool isbaij;
1634: PetscFunctionBegin;
1635: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &isbaij, MATSEQSBAIJ, MATMPISBAIJ, ""));
1636: PetscCheck(isbaij, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)B)->type_name);
1637: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1638: if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
1639: PetscCall(MatGetRowUpperTriangular(A));
1640: PetscCall(MatCopy_Basic(A, B, str));
1641: PetscCall(MatRestoreRowUpperTriangular(A));
1642: } else {
1643: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1644: Mat_MPISBAIJ *b = (Mat_MPISBAIJ *)B->data;
1646: PetscCall(MatCopy(a->A, b->A, str));
1647: PetscCall(MatCopy(a->B, b->B, str));
1648: }
1649: PetscCall(PetscObjectStateIncrease((PetscObject)B));
1650: PetscFunctionReturn(PETSC_SUCCESS);
1651: }
1653: static PetscErrorCode MatAXPY_MPISBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str)
1654: {
1655: Mat_MPISBAIJ *xx = (Mat_MPISBAIJ *)X->data, *yy = (Mat_MPISBAIJ *)Y->data;
1656: PetscBLASInt bnz, one = 1;
1657: Mat_SeqSBAIJ *xa, *ya;
1658: Mat_SeqBAIJ *xb, *yb;
1660: PetscFunctionBegin;
1661: if (str == SAME_NONZERO_PATTERN) {
1662: PetscScalar alpha = a;
1663: xa = (Mat_SeqSBAIJ *)xx->A->data;
1664: ya = (Mat_SeqSBAIJ *)yy->A->data;
1665: PetscCall(PetscBLASIntCast(xa->nz, &bnz));
1666: PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, xa->a, &one, ya->a, &one));
1667: xb = (Mat_SeqBAIJ *)xx->B->data;
1668: yb = (Mat_SeqBAIJ *)yy->B->data;
1669: PetscCall(PetscBLASIntCast(xb->nz, &bnz));
1670: PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, xb->a, &one, yb->a, &one));
1671: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
1672: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1673: PetscCall(MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE));
1674: PetscCall(MatAXPY_Basic(Y, a, X, str));
1675: PetscCall(MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE));
1676: } else {
1677: Mat B;
1678: PetscInt *nnz_d, *nnz_o, bs = Y->rmap->bs;
1679: PetscCheck(bs == X->rmap->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrices must have same block size");
1680: PetscCall(MatGetRowUpperTriangular(X));
1681: PetscCall(MatGetRowUpperTriangular(Y));
1682: PetscCall(PetscMalloc1(yy->A->rmap->N, &nnz_d));
1683: PetscCall(PetscMalloc1(yy->B->rmap->N, &nnz_o));
1684: PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B));
1685: PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name));
1686: PetscCall(MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N));
1687: PetscCall(MatSetBlockSizesFromMats(B, Y, Y));
1688: PetscCall(MatSetType(B, MATMPISBAIJ));
1689: PetscCall(MatAXPYGetPreallocation_SeqSBAIJ(yy->A, xx->A, nnz_d));
1690: PetscCall(MatAXPYGetPreallocation_MPIBAIJ(yy->B, yy->garray, xx->B, xx->garray, nnz_o));
1691: PetscCall(MatMPISBAIJSetPreallocation(B, bs, 0, nnz_d, 0, nnz_o));
1692: PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
1693: PetscCall(MatHeaderMerge(Y, &B));
1694: PetscCall(PetscFree(nnz_d));
1695: PetscCall(PetscFree(nnz_o));
1696: PetscCall(MatRestoreRowUpperTriangular(X));
1697: PetscCall(MatRestoreRowUpperTriangular(Y));
1698: }
1699: PetscFunctionReturn(PETSC_SUCCESS);
1700: }
1702: static PetscErrorCode MatCreateSubMatrices_MPISBAIJ(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
1703: {
1704: PetscInt i;
1705: PetscBool flg;
1707: PetscFunctionBegin;
1708: PetscCall(MatCreateSubMatrices_MPIBAIJ(A, n, irow, icol, scall, B)); /* B[] are sbaij matrices */
1709: for (i = 0; i < n; i++) {
1710: PetscCall(ISEqual(irow[i], icol[i], &flg));
1711: if (!flg) PetscCall(MatSeqSBAIJZeroOps_Private(*B[i]));
1712: }
1713: PetscFunctionReturn(PETSC_SUCCESS);
1714: }
1716: static PetscErrorCode MatShift_MPISBAIJ(Mat Y, PetscScalar a)
1717: {
1718: Mat_MPISBAIJ *maij = (Mat_MPISBAIJ *)Y->data;
1719: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)maij->A->data;
1721: PetscFunctionBegin;
1722: if (!Y->preallocated) {
1723: PetscCall(MatMPISBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL, 0, NULL));
1724: } else if (!aij->nz) {
1725: PetscInt nonew = aij->nonew;
1726: PetscCall(MatSeqSBAIJSetPreallocation(maij->A, Y->rmap->bs, 1, NULL));
1727: aij->nonew = nonew;
1728: }
1729: PetscCall(MatShift_Basic(Y, a));
1730: PetscFunctionReturn(PETSC_SUCCESS);
1731: }
1733: static PetscErrorCode MatZeroRowsColumns_MPISBAIJ(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1734: {
1735: Mat_MPISBAIJ *l = (Mat_MPISBAIJ *)A->data;
1736: PetscMPIInt n, p = 0;
1737: PetscInt i, j, k, r, len = 0, row, col, count;
1738: PetscInt *lrows, *owners = A->rmap->range;
1739: PetscSFNode *rrows;
1740: PetscSF sf;
1741: const PetscScalar *xx;
1742: PetscScalar *bb, *mask;
1743: Vec xmask, lmask, lvec_contrib = NULL;
1744: Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)l->B->data;
1745: PetscInt bs = A->rmap->bs, bs2 = baij->bs2;
1746: PetscScalar *aa;
1748: PetscFunctionBegin;
1749: PetscCall(PetscMPIIntCast(A->rmap->n, &n));
1750: /* create PetscSF where leaves are input rows and roots are owned rows */
1751: PetscCall(PetscMalloc1(n, &lrows));
1752: for (r = 0; r < n; ++r) lrows[r] = -1;
1753: PetscCall(PetscMalloc1(N, &rrows));
1754: for (r = 0; r < N; ++r) {
1755: const PetscInt idx = rows[r];
1756: 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);
1757: if (idx < owners[p] || owners[p + 1] <= idx) { /* short-circuit the search if the last p owns this row too */
1758: PetscCall(PetscLayoutFindOwner(A->rmap, idx, &p));
1759: }
1760: rrows[r].rank = p;
1761: rrows[r].index = rows[r] - owners[p];
1762: }
1763: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
1764: PetscCall(PetscSFSetGraph(sf, n, N, NULL, PETSC_OWN_POINTER, rrows, PETSC_OWN_POINTER));
1765: /* collect flags for rows to be zeroed */
1766: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, (PetscInt *)rows, lrows, MPI_LOR));
1767: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, (PetscInt *)rows, lrows, MPI_LOR));
1768: PetscCall(PetscSFDestroy(&sf));
1769: /* compress and put in row numbers */
1770: for (r = 0; r < n; ++r) {
1771: if (lrows[r] >= 0) lrows[len++] = r;
1772: }
1773: /* zero diagonal part of matrix */
1774: PetscCall(MatZeroRowsColumns(l->A, len, lrows, diag, x, b));
1775: /* handle off-diagonal part of matrix */
1776: PetscCall(MatCreateVecs(A, &xmask, NULL));
1777: PetscCall(VecDuplicate(l->lvec, &lmask));
1778: PetscCall(VecGetArray(xmask, &bb));
1779: for (i = 0; i < len; i++) bb[lrows[i]] = 1;
1780: PetscCall(VecRestoreArray(xmask, &bb));
1781: PetscCall(VecScatterBegin(l->Mvctx, xmask, lmask, ADD_VALUES, SCATTER_FORWARD));
1782: PetscCall(VecScatterEnd(l->Mvctx, xmask, lmask, ADD_VALUES, SCATTER_FORWARD));
1783: PetscCall(VecDestroy(&xmask));
1784: if (x) {
1785: PetscCall(VecScatterBegin(l->Mvctx, x, l->lvec, INSERT_VALUES, SCATTER_FORWARD));
1786: PetscCall(VecScatterEnd(l->Mvctx, x, l->lvec, INSERT_VALUES, SCATTER_FORWARD));
1787: PetscCall(VecGetArrayRead(l->lvec, &xx));
1788: PetscCall(VecGetArray(b, &bb));
1789: }
1790: PetscCall(VecGetArray(lmask, &mask));
1791: /* MPISBAIJ stores only the upper off-diagonal in l->B; for each zeroed local row r and
1792: non-zeroed off-process column c in that row, accumulate -A[r,c] * x[r] into lvec_contrib.
1793: A SCATTER_REVERSE below sends these contributions to b[c] on the owning (higher-rank)
1794: process, the missing symmetric lower-triangular update. We skip entries where c is
1795: also a zeroed row (mask[col] != 0) since b[c] = diag * x[c] is handled separately. */
1796: if (x) {
1797: const PetscScalar *x_vals;
1798: PetscScalar *c_vals;
1800: PetscCall(VecDuplicate(l->lvec, &lvec_contrib));
1801: PetscCall(VecGetArray(lvec_contrib, &c_vals));
1802: PetscCall(VecGetArrayRead(x, &x_vals));
1803: /* Only accumulate b[c] -= A[r,c] * x[r] when off-process col c is not also a zeroed row
1804: (mask[c] non-zero means col c is zeroed, so b[c] = diag * x[c] is already set).
1805: This mirrors the MatSeqSBAIJ pattern: if (zeroed[r] && !zeroed[c]) bb[c] -= A[r,c] * x[r].
1806: c_vals is indexed by the local B column index. */
1807: for (i = 0; i < len; ++i) {
1808: row = lrows[i];
1809: for (j = baij->i[row / bs]; j < baij->i[row / bs + 1]; ++j) {
1810: for (k = 0; k < bs; ++k) {
1811: col = baij->j[j] * bs + k;
1812: if (!PetscAbsScalar(mask[col])) {
1813: aa = baij->a + j * bs2 + (row % bs) + bs * k;
1814: c_vals[col] -= aa[0] * x_vals[row];
1815: }
1816: }
1817: }
1818: }
1819: PetscCall(VecRestoreArrayRead(x, &x_vals));
1820: PetscCall(VecRestoreArray(lvec_contrib, &c_vals));
1821: }
1822: /* remove zeroed rows of off-diagonal matrix */
1823: for (i = 0; i < len; ++i) {
1824: row = lrows[i];
1825: count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
1826: aa = baij->a + baij->i[row / bs] * bs2 + (row % bs);
1827: for (k = 0; k < count; ++k) {
1828: aa[0] = 0.0;
1829: aa += bs;
1830: }
1831: }
1832: /* loop over all elements of off process part of matrix zeroing removed columns */
1833: for (i = 0; i < l->B->rmap->N; ++i) {
1834: row = i / bs;
1835: for (j = baij->i[row]; j < baij->i[row + 1]; ++j) {
1836: for (k = 0; k < bs; ++k) {
1837: col = bs * baij->j[j] + k;
1838: if (PetscAbsScalar(mask[col])) {
1839: aa = baij->a + j * bs2 + (i % bs) + bs * k;
1840: if (x) bb[i] -= aa[0] * xx[col];
1841: aa[0] = 0.0;
1842: }
1843: }
1844: }
1845: }
1846: if (x) {
1847: PetscCall(VecRestoreArray(b, &bb));
1848: PetscCall(VecRestoreArrayRead(l->lvec, &xx));
1849: /* scatter the accumulated contributions to b[c] on higher-rank processes owning column c */
1850: PetscCall(VecScatterBegin(l->Mvctx, lvec_contrib, b, ADD_VALUES, SCATTER_REVERSE));
1851: PetscCall(VecScatterEnd(l->Mvctx, lvec_contrib, b, ADD_VALUES, SCATTER_REVERSE));
1852: PetscCall(VecDestroy(&lvec_contrib));
1853: }
1854: PetscCall(VecRestoreArray(lmask, &mask));
1855: PetscCall(VecDestroy(&lmask));
1856: PetscCall(PetscFree(lrows));
1858: /* only change matrix nonzero state if pattern was allowed to be changed */
1859: if (!((Mat_SeqSBAIJ *)l->A->data)->nonew) {
1860: PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate;
1861: PetscCallMPI(MPIU_Allreduce(&state, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A)));
1862: }
1863: PetscFunctionReturn(PETSC_SUCCESS);
1864: }
1866: static PetscErrorCode MatGetDiagonalBlock_MPISBAIJ(Mat A, Mat *a)
1867: {
1868: PetscFunctionBegin;
1869: *a = ((Mat_MPISBAIJ *)A->data)->A;
1870: PetscFunctionReturn(PETSC_SUCCESS);
1871: }
1873: static PetscErrorCode MatEliminateZeros_MPISBAIJ(Mat A, PetscBool keep)
1874: {
1875: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
1877: PetscFunctionBegin;
1878: PetscCall(MatEliminateZeros_SeqSBAIJ(a->A, keep)); // possibly keep zero diagonal coefficients
1879: PetscCall(MatEliminateZeros_SeqBAIJ(a->B, PETSC_FALSE)); // never keep zero diagonal coefficients
1880: PetscFunctionReturn(PETSC_SUCCESS);
1881: }
1883: static PetscErrorCode MatLoad_MPISBAIJ(Mat, PetscViewer);
1884: static PetscErrorCode MatGetRowMaxAbs_MPISBAIJ(Mat, Vec, PetscInt[]);
1885: static PetscErrorCode MatSOR_MPISBAIJ(Mat, Vec, PetscReal, MatSORType, PetscReal, PetscInt, PetscInt, Vec);
1887: static struct _MatOps MatOps_Values = {MatSetValues_MPISBAIJ,
1888: MatGetRow_MPISBAIJ,
1889: MatRestoreRow_MPISBAIJ,
1890: MatMult_MPISBAIJ,
1891: /* 4*/ MatMultAdd_MPISBAIJ,
1892: MatMult_MPISBAIJ, /* transpose versions are same as non-transpose */
1893: MatMultAdd_MPISBAIJ,
1894: NULL,
1895: NULL,
1896: NULL,
1897: /* 10*/ NULL,
1898: NULL,
1899: NULL,
1900: MatSOR_MPISBAIJ,
1901: MatTranspose_MPISBAIJ,
1902: /* 15*/ MatGetInfo_MPISBAIJ,
1903: MatEqual_MPISBAIJ,
1904: MatGetDiagonal_MPISBAIJ,
1905: MatDiagonalScale_MPISBAIJ,
1906: MatNorm_MPISBAIJ,
1907: /* 20*/ MatAssemblyBegin_MPISBAIJ,
1908: MatAssemblyEnd_MPISBAIJ,
1909: MatSetOption_MPISBAIJ,
1910: MatZeroEntries_MPISBAIJ,
1911: /* 24*/ NULL,
1912: NULL,
1913: NULL,
1914: NULL,
1915: NULL,
1916: /* 29*/ MatSetUp_MPI_Hash,
1917: NULL,
1918: NULL,
1919: MatGetDiagonalBlock_MPISBAIJ,
1920: NULL,
1921: /* 34*/ MatDuplicate_MPISBAIJ,
1922: NULL,
1923: NULL,
1924: NULL,
1925: NULL,
1926: /* 39*/ MatAXPY_MPISBAIJ,
1927: MatCreateSubMatrices_MPISBAIJ,
1928: MatIncreaseOverlap_MPISBAIJ,
1929: MatGetValues_MPISBAIJ,
1930: MatCopy_MPISBAIJ,
1931: /* 44*/ NULL,
1932: MatScale_MPISBAIJ,
1933: MatShift_MPISBAIJ,
1934: NULL,
1935: MatZeroRowsColumns_MPISBAIJ,
1936: /* 49*/ NULL,
1937: NULL,
1938: NULL,
1939: NULL,
1940: NULL,
1941: /* 54*/ NULL,
1942: NULL,
1943: MatSetUnfactored_MPISBAIJ,
1944: NULL,
1945: MatSetValuesBlocked_MPISBAIJ,
1946: /* 59*/ MatCreateSubMatrix_MPISBAIJ,
1947: NULL,
1948: NULL,
1949: NULL,
1950: NULL,
1951: /* 64*/ NULL,
1952: NULL,
1953: NULL,
1954: NULL,
1955: MatGetRowMaxAbs_MPISBAIJ,
1956: /* 69*/ NULL,
1957: MatConvert_MPISBAIJ_Basic,
1958: NULL,
1959: NULL,
1960: NULL,
1961: NULL,
1962: NULL,
1963: NULL,
1964: NULL,
1965: MatLoad_MPISBAIJ,
1966: /* 79*/ NULL,
1967: NULL,
1968: NULL,
1969: NULL,
1970: NULL,
1971: /* 84*/ NULL,
1972: NULL,
1973: NULL,
1974: NULL,
1975: NULL,
1976: /* 89*/ NULL,
1977: NULL,
1978: NULL,
1979: NULL,
1980: MatConjugate_MPISBAIJ,
1981: /* 94*/ NULL,
1982: NULL,
1983: MatRealPart_MPISBAIJ,
1984: MatImaginaryPart_MPISBAIJ,
1985: MatGetRowUpperTriangular_MPISBAIJ,
1986: /* 99*/ MatRestoreRowUpperTriangular_MPISBAIJ,
1987: NULL,
1988: NULL,
1989: NULL,
1990: NULL,
1991: /*104*/ NULL,
1992: NULL,
1993: NULL,
1994: NULL,
1995: NULL,
1996: /*109*/ NULL,
1997: NULL,
1998: NULL,
1999: NULL,
2000: NULL,
2001: /*114*/ NULL,
2002: NULL,
2003: NULL,
2004: NULL,
2005: NULL,
2006: /*119*/ NULL,
2007: NULL,
2008: NULL,
2009: NULL,
2010: NULL,
2011: /*124*/ NULL,
2012: NULL,
2013: MatSetBlockSizes_Default,
2014: NULL,
2015: NULL,
2016: /*129*/ NULL,
2017: MatCreateMPIMatConcatenateSeqMat_MPISBAIJ,
2018: NULL,
2019: NULL,
2020: NULL,
2021: /*134*/ NULL,
2022: NULL,
2023: MatEliminateZeros_MPISBAIJ,
2024: NULL,
2025: NULL,
2026: /*139*/ NULL,
2027: NULL,
2028: MatCopyHashToXAIJ_MPI_Hash,
2029: NULL,
2030: NULL,
2031: /*144*/ MatADot_Default,
2032: MatANorm_Default,
2033: NULL,
2034: NULL};
2036: static PetscErrorCode MatMPISBAIJSetPreallocation_MPISBAIJ(Mat B, PetscInt bs, PetscInt d_nz, const PetscInt *d_nnz, PetscInt o_nz, const PetscInt *o_nnz)
2037: {
2038: Mat_MPISBAIJ *b = (Mat_MPISBAIJ *)B->data;
2039: PetscInt i, mbs, Mbs;
2040: PetscMPIInt size;
2042: PetscFunctionBegin;
2043: if (B->hash_active) {
2044: B->ops[0] = b->cops;
2045: B->hash_active = PETSC_FALSE;
2046: }
2047: if (!B->preallocated) PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), bs, &B->bstash));
2048: PetscCall(MatSetBlockSize(B, bs));
2049: PetscCall(PetscLayoutSetUp(B->rmap));
2050: PetscCall(PetscLayoutSetUp(B->cmap));
2051: PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
2052: PetscCheck(B->rmap->N <= B->cmap->N, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "MPISBAIJ matrix cannot have more rows %" PetscInt_FMT " than columns %" PetscInt_FMT, B->rmap->N, B->cmap->N);
2053: PetscCheck(B->rmap->n <= B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "MPISBAIJ matrix cannot have more local rows %" PetscInt_FMT " than columns %" PetscInt_FMT, B->rmap->n, B->cmap->n);
2055: mbs = B->rmap->n / bs;
2056: Mbs = B->rmap->N / bs;
2057: PetscCheck(mbs * bs == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "No of local rows %" PetscInt_FMT " must be divisible by blocksize %" PetscInt_FMT, B->rmap->N, bs);
2059: B->rmap->bs = bs;
2060: b->bs2 = bs * bs;
2061: b->mbs = mbs;
2062: b->Mbs = Mbs;
2063: b->nbs = B->cmap->n / bs;
2064: b->Nbs = B->cmap->N / bs;
2066: for (i = 0; i <= b->size; i++) b->rangebs[i] = B->rmap->range[i] / bs;
2067: b->rstartbs = B->rmap->rstart / bs;
2068: b->rendbs = B->rmap->rend / bs;
2070: b->cstartbs = B->cmap->rstart / bs;
2071: b->cendbs = B->cmap->rend / bs;
2073: #if defined(PETSC_USE_CTABLE)
2074: PetscCall(PetscHMapIDestroy(&b->colmap));
2075: #else
2076: PetscCall(PetscFree(b->colmap));
2077: #endif
2078: PetscCall(PetscFree(b->garray));
2079: PetscCall(VecDestroy(&b->lvec));
2080: PetscCall(VecScatterDestroy(&b->Mvctx));
2081: PetscCall(VecDestroy(&b->slvec0));
2082: PetscCall(VecDestroy(&b->slvec0b));
2083: PetscCall(VecDestroy(&b->slvec1));
2084: PetscCall(VecDestroy(&b->slvec1a));
2085: PetscCall(VecDestroy(&b->slvec1b));
2086: PetscCall(VecScatterDestroy(&b->sMvctx));
2088: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
2090: MatSeqXAIJGetOptions_Private(b->B);
2091: PetscCall(MatDestroy(&b->B));
2092: PetscCall(MatCreate(PETSC_COMM_SELF, &b->B));
2093: PetscCall(MatSetSizes(b->B, B->rmap->n, size > 1 ? B->cmap->N : 0, B->rmap->n, size > 1 ? B->cmap->N : 0));
2094: PetscCall(MatSetType(b->B, MATSEQBAIJ));
2095: MatSeqXAIJRestoreOptions_Private(b->B);
2097: MatSeqXAIJGetOptions_Private(b->A);
2098: PetscCall(MatDestroy(&b->A));
2099: PetscCall(MatCreate(PETSC_COMM_SELF, &b->A));
2100: PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n));
2101: PetscCall(MatSetType(b->A, MATSEQSBAIJ));
2102: MatSeqXAIJRestoreOptions_Private(b->A);
2104: PetscCall(MatSeqSBAIJSetPreallocation(b->A, bs, d_nz, d_nnz));
2105: PetscCall(MatSeqBAIJSetPreallocation(b->B, bs, o_nz, o_nnz));
2107: B->preallocated = PETSC_TRUE;
2108: B->was_assembled = PETSC_FALSE;
2109: B->assembled = PETSC_FALSE;
2110: PetscFunctionReturn(PETSC_SUCCESS);
2111: }
2113: static PetscErrorCode MatMPISBAIJSetPreallocationCSR_MPISBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[])
2114: {
2115: PetscInt m, rstart, cend;
2116: PetscInt i, j, d, nz, bd, nz_max = 0, *d_nnz = NULL, *o_nnz = NULL;
2117: const PetscInt *JJ = NULL;
2118: PetscScalar *values = NULL;
2119: PetscBool roworiented = ((Mat_MPISBAIJ *)B->data)->roworiented;
2120: PetscBool nooffprocentries;
2122: PetscFunctionBegin;
2123: PetscCheck(bs >= 1, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "Invalid block size specified, must be positive but it is %" PetscInt_FMT, bs);
2124: PetscCall(PetscLayoutSetBlockSize(B->rmap, bs));
2125: PetscCall(PetscLayoutSetBlockSize(B->cmap, bs));
2126: PetscCall(PetscLayoutSetUp(B->rmap));
2127: PetscCall(PetscLayoutSetUp(B->cmap));
2128: PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
2129: m = B->rmap->n / bs;
2130: rstart = B->rmap->rstart / bs;
2131: cend = B->cmap->rend / bs;
2133: PetscCheck(!ii[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]);
2134: PetscCall(PetscMalloc2(m, &d_nnz, m, &o_nnz));
2135: for (i = 0; i < m; i++) {
2136: nz = ii[i + 1] - ii[i];
2137: PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz);
2138: /* count the ones on the diagonal and above, split into diagonal and off-diagonal portions. */
2139: JJ = jj + ii[i];
2140: bd = 0;
2141: for (j = 0; j < nz; j++) {
2142: if (*JJ >= i + rstart) break;
2143: JJ++;
2144: bd++;
2145: }
2146: d = 0;
2147: for (; j < nz; j++) {
2148: if (*JJ++ >= cend) break;
2149: d++;
2150: }
2151: d_nnz[i] = d;
2152: o_nnz[i] = nz - d - bd;
2153: nz = nz - bd;
2154: nz_max = PetscMax(nz_max, nz);
2155: }
2156: PetscCall(MatMPISBAIJSetPreallocation(B, bs, 0, d_nnz, 0, o_nnz));
2157: PetscCall(MatSetOption(B, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
2158: PetscCall(PetscFree2(d_nnz, o_nnz));
2160: values = (PetscScalar *)V;
2161: if (!values) PetscCall(PetscCalloc1(bs * bs * nz_max, &values));
2162: for (i = 0; i < m; i++) {
2163: PetscInt row = i + rstart;
2164: PetscInt ncols = ii[i + 1] - ii[i];
2165: const PetscInt *icols = jj + ii[i];
2166: if (bs == 1 || !roworiented) { /* block ordering matches the non-nested layout of MatSetValues so we can insert entire rows */
2167: const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0);
2168: PetscCall(MatSetValuesBlocked_MPISBAIJ(B, 1, &row, ncols, icols, svals, INSERT_VALUES));
2169: } else { /* block ordering does not match so we can only insert one block at a time. */
2170: PetscInt j;
2171: for (j = 0; j < ncols; j++) {
2172: const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0);
2173: PetscCall(MatSetValuesBlocked_MPISBAIJ(B, 1, &row, 1, &icols[j], svals, INSERT_VALUES));
2174: }
2175: }
2176: }
2178: if (!V) PetscCall(PetscFree(values));
2179: nooffprocentries = B->nooffprocentries;
2180: B->nooffprocentries = PETSC_TRUE;
2181: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2182: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2183: B->nooffprocentries = nooffprocentries;
2185: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
2186: PetscFunctionReturn(PETSC_SUCCESS);
2187: }
2189: /*MC
2190: MATMPISBAIJ - MATMPISBAIJ = "mpisbaij" - A matrix type to be used for distributed symmetric sparse block matrices,
2191: based on block compressed sparse row format. Only the upper triangular portion of the "diagonal" portion of
2192: the matrix is stored.
2194: For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
2195: can call `MatSetOption`(`Mat`, `MAT_HERMITIAN`);
2197: Options Database Key:
2198: . -mat_type mpisbaij - sets the matrix type to "mpisbaij" during a call to `MatSetFromOptions()`
2200: Level: beginner
2202: Note:
2203: The number of rows in the matrix must be less than or equal to the number of columns. Similarly the number of rows in the
2204: diagonal portion of the matrix of each process has to less than or equal the number of columns.
2206: .seealso: [](ch_matrices), `Mat`, `MATSBAIJ`, `MATBAIJ`, `MatCreateBAIJ()`, `MATSEQSBAIJ`, `MatType`
2207: M*/
2209: PETSC_EXTERN PetscErrorCode MatCreate_MPISBAIJ(Mat B)
2210: {
2211: Mat_MPISBAIJ *b;
2212: PetscBool flg = PETSC_FALSE;
2214: PetscFunctionBegin;
2215: PetscCall(PetscNew(&b));
2216: B->data = (void *)b;
2217: B->ops[0] = MatOps_Values;
2219: B->ops->destroy = MatDestroy_MPISBAIJ;
2220: B->ops->view = MatView_MPISBAIJ;
2221: B->assembled = PETSC_FALSE;
2222: B->insertmode = NOT_SET_VALUES;
2224: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank));
2225: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &b->size));
2227: /* build local table of row and column ownerships */
2228: PetscCall(PetscMalloc1(b->size + 2, &b->rangebs));
2230: /* build cache for off array entries formed */
2231: PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash));
2233: b->donotstash = PETSC_FALSE;
2234: b->colmap = NULL;
2235: b->garray = NULL;
2236: b->roworiented = PETSC_TRUE;
2238: /* stuff used in block assembly */
2239: b->barray = NULL;
2241: /* stuff used for matrix vector multiply */
2242: b->lvec = NULL;
2243: b->Mvctx = NULL;
2244: b->slvec0 = NULL;
2245: b->slvec0b = NULL;
2246: b->slvec1 = NULL;
2247: b->slvec1a = NULL;
2248: b->slvec1b = NULL;
2249: b->sMvctx = NULL;
2251: /* stuff for MatGetRow() */
2252: b->rowindices = NULL;
2253: b->rowvalues = NULL;
2254: b->getrowactive = PETSC_FALSE;
2256: /* hash table stuff */
2257: b->ht = NULL;
2258: b->hd = NULL;
2259: b->ht_size = 0;
2260: b->ht_flag = PETSC_FALSE;
2261: b->ht_fact = 0;
2262: b->ht_total_ct = 0;
2263: b->ht_insert_ct = 0;
2265: /* stuff for MatCreateSubMatrices_MPIBAIJ_local() */
2266: b->ijonly = PETSC_FALSE;
2268: b->in_loc = NULL;
2269: b->v_loc = NULL;
2270: b->n_loc = 0;
2272: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPISBAIJ));
2273: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPISBAIJ));
2274: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISBAIJSetPreallocation_C", MatMPISBAIJSetPreallocation_MPISBAIJ));
2275: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISBAIJSetPreallocationCSR_C", MatMPISBAIJSetPreallocationCSR_MPISBAIJ));
2276: #if defined(PETSC_HAVE_ELEMENTAL)
2277: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisbaij_elemental_C", MatConvert_MPISBAIJ_Elemental));
2278: #endif
2279: #if defined(PETSC_HAVE_SCALAPACK) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
2280: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisbaij_scalapack_C", MatConvert_SBAIJ_ScaLAPACK));
2281: #endif
2282: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisbaij_mpiaij_C", MatConvert_MPISBAIJ_Basic));
2283: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisbaij_mpibaij_C", MatConvert_MPISBAIJ_Basic));
2285: B->symmetric = PETSC_BOOL3_TRUE;
2286: B->structurally_symmetric = PETSC_BOOL3_TRUE;
2287: B->symmetry_eternal = PETSC_TRUE;
2288: B->structural_symmetry_eternal = PETSC_TRUE;
2289: #if !defined(PETSC_USE_COMPLEX)
2290: B->hermitian = PETSC_BOOL3_TRUE;
2291: #endif
2293: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPISBAIJ));
2294: PetscOptionsBegin(PetscObjectComm((PetscObject)B), NULL, "Options for loading MPISBAIJ matrix 1", "Mat");
2295: PetscCall(PetscOptionsBool("-mat_use_hash_table", "Use hash table to save memory in constructing matrix", "MatSetOption", flg, &flg, NULL));
2296: if (flg) {
2297: PetscReal fact = 1.39;
2298: PetscCall(MatSetOption(B, MAT_USE_HASH_TABLE, PETSC_TRUE));
2299: PetscCall(PetscOptionsReal("-mat_use_hash_table", "Use hash table factor", "MatMPIBAIJSetHashTableFactor", fact, &fact, NULL));
2300: if (fact <= 1.0) fact = 1.39;
2301: PetscCall(MatMPIBAIJSetHashTableFactor(B, fact));
2302: PetscCall(PetscInfo(B, "Hash table Factor used %5.2g\n", (double)fact));
2303: }
2304: PetscOptionsEnd();
2305: PetscFunctionReturn(PETSC_SUCCESS);
2306: }
2308: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
2309: /*MC
2310: MATSBAIJ - MATSBAIJ = "sbaij" - A matrix type to be used for symmetric block sparse matrices.
2312: This matrix type is identical to `MATSEQSBAIJ` when constructed with a single process communicator,
2313: and `MATMPISBAIJ` otherwise.
2315: Options Database Key:
2316: . -mat_type sbaij - sets the matrix type to `MATSBAIJ` during a call to `MatSetFromOptions()`
2318: Level: beginner
2320: .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MATMPISBAIJ`, `MatCreateSBAIJ()`, `MATSEQBAIJ`, `MATMPIBAIJ`
2321: M*/
2323: /*@
2324: MatMPISBAIJSetPreallocation - For good matrix assembly performance
2325: the user should preallocate the matrix storage by setting the parameters
2326: d_nz (or d_nnz) and o_nz (or o_nnz). By setting these parameters accurately,
2327: performance can be increased by more than a factor of 50.
2329: Collective
2331: Input Parameters:
2332: + B - the matrix
2333: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2334: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2335: . d_nz - number of block nonzeros per block row in diagonal portion of local
2336: submatrix (same for all local rows)
2337: . d_nnz - array containing the number of block nonzeros in the various block rows
2338: in the upper triangular and diagonal part of the in diagonal portion of the local
2339: (possibly different for each block row) or `NULL`. If you plan to factor the matrix you must leave room
2340: for the diagonal entry and set a value even if it is zero.
2341: . o_nz - number of block nonzeros per block row in the off-diagonal portion of local
2342: submatrix (same for all local rows).
2343: - o_nnz - array containing the number of nonzeros in the various block rows of the
2344: off-diagonal portion of the local submatrix that is right of the diagonal
2345: (possibly different for each block row) or `NULL`.
2347: Options Database Keys:
2348: + -mat_no_unroll - uses code that does not unroll the loops in the
2349: block calculations (much slower)
2350: - -mat_block_size - size of the blocks to use
2352: Level: intermediate
2354: Notes:
2356: If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one processor
2357: than it must be used on all processors that share the object for that argument.
2359: If the *_nnz parameter is given then the *_nz parameter is ignored
2361: Storage Information:
2362: For a square global matrix we define each processor's diagonal portion
2363: to be its local rows and the corresponding columns (a square submatrix);
2364: each processor's off-diagonal portion encompasses the remainder of the
2365: local matrix (a rectangular submatrix).
2367: The user can specify preallocated storage for the diagonal part of
2368: the local submatrix with either `d_nz` or `d_nnz` (not both). Set
2369: `d_nz` = `PETSC_DEFAULT` and `d_nnz` = `NULL` for PETSc to control dynamic
2370: memory allocation. Likewise, specify preallocated storage for the
2371: off-diagonal part of the local submatrix with `o_nz` or `o_nnz` (not both).
2373: You can call `MatGetInfo()` to get information on how effective the preallocation was;
2374: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
2375: You can also run with the option `-info` and look for messages with the string
2376: malloc in them to see if additional memory allocation was needed.
2378: Consider a processor that owns rows 3, 4 and 5 of a parallel matrix. In
2379: the figure below we depict these three local rows and all columns (0-11).
2381: .vb
2382: 0 1 2 3 4 5 6 7 8 9 10 11
2383: --------------------------
2384: row 3 |. . . d d d o o o o o o
2385: row 4 |. . . d d d o o o o o o
2386: row 5 |. . . d d d o o o o o o
2387: --------------------------
2388: .ve
2390: Thus, any entries in the d locations are stored in the d (diagonal)
2391: submatrix, and any entries in the o locations are stored in the
2392: o (off-diagonal) submatrix. Note that the d matrix is stored in
2393: `MATSEQSBAIJ` format and the o submatrix in `MATSEQBAIJ` format.
2395: Now `d_nz` should indicate the number of block nonzeros per row in the upper triangular
2396: plus the diagonal part of the d matrix,
2397: and `o_nz` should indicate the number of block nonzeros per row in the o matrix
2399: In general, for PDE problems in which most nonzeros are near the diagonal,
2400: one expects `d_nz` >> `o_nz`.
2402: .seealso: [](ch_matrices), `Mat`, `MATMPISBAIJ`, `MATSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`, `PetscSplitOwnership()`
2403: @*/
2404: PetscErrorCode MatMPISBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
2405: {
2406: PetscFunctionBegin;
2410: PetscTryMethod(B, "MatMPISBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, bs, d_nz, d_nnz, o_nz, o_nnz));
2411: PetscFunctionReturn(PETSC_SUCCESS);
2412: }
2414: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
2415: /*@
2416: MatCreateSBAIJ - Creates a sparse parallel matrix in symmetric block AIJ format, `MATSBAIJ`,
2417: (block compressed row). For good matrix assembly performance
2418: the user should preallocate the matrix storage by setting the parameters
2419: `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`).
2421: Collective
2423: Input Parameters:
2424: + comm - MPI communicator
2425: . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
2426: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
2427: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
2428: This value should be the same as the local size used in creating the
2429: y vector for the matrix-vector product y = Ax.
2430: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
2431: This value should be the same as the local size used in creating the
2432: x vector for the matrix-vector product y = Ax.
2433: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
2434: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
2435: . d_nz - number of block nonzeros per block row in diagonal portion of local
2436: submatrix (same for all local rows)
2437: . d_nnz - array containing the number of block nonzeros in the various block rows
2438: in the upper triangular portion of the in diagonal portion of the local
2439: (possibly different for each block block row) or `NULL`.
2440: If you plan to factor the matrix you must leave room for the diagonal entry and
2441: set its value even if it is zero.
2442: . o_nz - number of block nonzeros per block row in the off-diagonal portion of local
2443: submatrix (same for all local rows).
2444: - o_nnz - array containing the number of nonzeros in the various block rows of the
2445: off-diagonal portion of the local submatrix (possibly different for
2446: each block row) or `NULL`.
2448: Output Parameter:
2449: . A - the matrix
2451: Options Database Keys:
2452: + -mat_no_unroll - uses code that does not unroll the loops in the
2453: block calculations (much slower)
2454: . -mat_block_size - size of the blocks to use
2455: - -mat_mpi - use the parallel matrix data structures even on one processor
2456: (defaults to using SeqBAIJ format on one processor)
2458: Level: intermediate
2460: Notes:
2461: It is recommended that one use `MatCreateFromOptions()` or the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
2462: MatXXXXSetPreallocation() paradigm instead of this routine directly.
2463: [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`]
2465: The number of rows and columns must be divisible by blocksize.
2466: This matrix type does not support complex Hermitian operation.
2468: The user MUST specify either the local or global matrix dimensions
2469: (possibly both).
2471: If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one processor
2472: than it must be used on all processors that share the object for that argument.
2474: If `m` and `n` are not `PETSC_DECIDE`, then the values determines the `PetscLayout` of the matrix and the ranges returned by
2475: `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `MatGetOwnershipRangeColumn()`, and `MatGetOwnershipRangesColumn()`.
2477: If the *_nnz parameter is given then the *_nz parameter is ignored
2479: Storage Information:
2480: For a square global matrix we define each processor's diagonal portion
2481: to be its local rows and the corresponding columns (a square submatrix);
2482: each processor's off-diagonal portion encompasses the remainder of the
2483: local matrix (a rectangular submatrix).
2485: The user can specify preallocated storage for the diagonal part of
2486: the local submatrix with either `d_nz` or `d_nnz` (not both). Set
2487: `d_nz` = `PETSC_DEFAULT` and `d_nnz` = `NULL` for PETSc to control dynamic
2488: memory allocation. Likewise, specify preallocated storage for the
2489: off-diagonal part of the local submatrix with `o_nz` or `o_nnz` (not both).
2491: Consider a processor that owns rows 3, 4 and 5 of a parallel matrix. In
2492: the figure below we depict these three local rows and all columns (0-11).
2494: .vb
2495: 0 1 2 3 4 5 6 7 8 9 10 11
2496: --------------------------
2497: row 3 |. . . d d d o o o o o o
2498: row 4 |. . . d d d o o o o o o
2499: row 5 |. . . d d d o o o o o o
2500: --------------------------
2501: .ve
2503: Thus, any entries in the d locations are stored in the d (diagonal)
2504: submatrix, and any entries in the o locations are stored in the
2505: o (off-diagonal) submatrix. Note that the d matrix is stored in
2506: `MATSEQSBAIJ` format and the o submatrix in `MATSEQBAIJ` format.
2508: Now `d_nz` should indicate the number of block nonzeros per row in the upper triangular
2509: plus the diagonal part of the d matrix,
2510: and `o_nz` should indicate the number of block nonzeros per row in the o matrix.
2511: In general, for PDE problems in which most nonzeros are near the diagonal,
2512: one expects `d_nz` >> `o_nz`.
2514: .seealso: [](ch_matrices), `Mat`, `MATSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`,
2515: `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `MatGetOwnershipRangeColumn()`, `MatGetOwnershipRangesColumn()`, `PetscLayout`
2516: @*/
2517: PetscErrorCode MatCreateSBAIJ(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)
2518: {
2519: PetscMPIInt size;
2521: PetscFunctionBegin;
2522: PetscCall(MatCreate(comm, A));
2523: PetscCall(MatSetSizes(*A, m, n, M, N));
2524: PetscCallMPI(MPI_Comm_size(comm, &size));
2525: if (size > 1) {
2526: PetscCall(MatSetType(*A, MATMPISBAIJ));
2527: PetscCall(MatMPISBAIJSetPreallocation(*A, bs, d_nz, d_nnz, o_nz, o_nnz));
2528: } else {
2529: PetscCall(MatSetType(*A, MATSEQSBAIJ));
2530: PetscCall(MatSeqSBAIJSetPreallocation(*A, bs, d_nz, d_nnz));
2531: }
2532: PetscFunctionReturn(PETSC_SUCCESS);
2533: }
2535: static PetscErrorCode MatDuplicate_MPISBAIJ(Mat matin, MatDuplicateOption cpvalues, Mat *newmat)
2536: {
2537: Mat mat;
2538: Mat_MPISBAIJ *a, *oldmat = (Mat_MPISBAIJ *)matin->data;
2539: PetscInt len = 0, nt, bs = matin->rmap->bs, mbs = oldmat->mbs;
2540: PetscScalar *array;
2542: PetscFunctionBegin;
2543: *newmat = NULL;
2545: PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat));
2546: PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N));
2547: PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name));
2548: PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap));
2549: PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap));
2551: if (matin->hash_active) PetscCall(MatSetUp(mat));
2552: else {
2553: mat->factortype = matin->factortype;
2554: mat->preallocated = PETSC_TRUE;
2555: mat->assembled = PETSC_TRUE;
2556: mat->insertmode = NOT_SET_VALUES;
2558: a = (Mat_MPISBAIJ *)mat->data;
2559: a->bs2 = oldmat->bs2;
2560: a->mbs = oldmat->mbs;
2561: a->nbs = oldmat->nbs;
2562: a->Mbs = oldmat->Mbs;
2563: a->Nbs = oldmat->Nbs;
2565: a->size = oldmat->size;
2566: a->rank = oldmat->rank;
2567: a->donotstash = oldmat->donotstash;
2568: a->roworiented = oldmat->roworiented;
2569: a->rowindices = NULL;
2570: a->rowvalues = NULL;
2571: a->getrowactive = PETSC_FALSE;
2572: a->barray = NULL;
2573: a->rstartbs = oldmat->rstartbs;
2574: a->rendbs = oldmat->rendbs;
2575: a->cstartbs = oldmat->cstartbs;
2576: a->cendbs = oldmat->cendbs;
2578: /* hash table stuff */
2579: a->ht = NULL;
2580: a->hd = NULL;
2581: a->ht_size = 0;
2582: a->ht_flag = oldmat->ht_flag;
2583: a->ht_fact = oldmat->ht_fact;
2584: a->ht_total_ct = 0;
2585: a->ht_insert_ct = 0;
2587: PetscCall(PetscArraycpy(a->rangebs, oldmat->rangebs, a->size + 2));
2588: if (oldmat->colmap) {
2589: #if defined(PETSC_USE_CTABLE)
2590: PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap));
2591: #else
2592: PetscCall(PetscMalloc1(a->Nbs, &a->colmap));
2593: PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, a->Nbs));
2594: #endif
2595: } else a->colmap = NULL;
2597: if (oldmat->garray && (len = ((Mat_SeqBAIJ *)oldmat->B->data)->nbs)) {
2598: PetscCall(PetscMalloc1(len, &a->garray));
2599: PetscCall(PetscArraycpy(a->garray, oldmat->garray, len));
2600: } else a->garray = NULL;
2602: PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)matin), matin->rmap->bs, &mat->bstash));
2603: PetscCall(VecDuplicate(oldmat->lvec, &a->lvec));
2604: PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx));
2606: PetscCall(VecDuplicate(oldmat->slvec0, &a->slvec0));
2607: PetscCall(VecDuplicate(oldmat->slvec1, &a->slvec1));
2609: PetscCall(VecGetLocalSize(a->slvec1, &nt));
2610: PetscCall(VecGetArray(a->slvec1, &array));
2611: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, bs * mbs, array, &a->slvec1a));
2612: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, nt - bs * mbs, array + bs * mbs, &a->slvec1b));
2613: PetscCall(VecRestoreArray(a->slvec1, &array));
2614: PetscCall(VecGetArray(a->slvec0, &array));
2615: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, nt - bs * mbs, array + bs * mbs, &a->slvec0b));
2616: PetscCall(VecRestoreArray(a->slvec0, &array));
2618: /* ierr = VecScatterCopy(oldmat->sMvctx,&a->sMvctx); - not written yet, replaced by the lazy trick: */
2619: PetscCall(PetscObjectReference((PetscObject)oldmat->sMvctx));
2620: a->sMvctx = oldmat->sMvctx;
2622: PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A));
2623: PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B));
2624: }
2625: PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist));
2626: *newmat = mat;
2627: PetscFunctionReturn(PETSC_SUCCESS);
2628: }
2630: /* Used for both MPIBAIJ and MPISBAIJ matrices */
2631: #define MatLoad_MPISBAIJ_Binary MatLoad_MPIBAIJ_Binary
2633: static PetscErrorCode MatLoad_MPISBAIJ(Mat mat, PetscViewer viewer)
2634: {
2635: PetscBool isbinary;
2637: PetscFunctionBegin;
2638: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2639: 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);
2640: PetscCall(MatLoad_MPISBAIJ_Binary(mat, viewer));
2641: PetscFunctionReturn(PETSC_SUCCESS);
2642: }
2644: static PetscErrorCode MatGetRowMaxAbs_MPISBAIJ(Mat A, Vec v, PetscInt idx[])
2645: {
2646: Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
2647: Mat_SeqBAIJ *b = (Mat_SeqBAIJ *)a->B->data;
2648: PetscReal atmp;
2649: PetscReal *work, *svalues, *rvalues;
2650: PetscInt i, bs, mbs, *bi, *bj, brow, j, ncols, krow, kcol, col, row, Mbs, bcol;
2651: PetscMPIInt rank, size;
2652: PetscInt *rowners_bs, count, source;
2653: PetscScalar *va;
2654: MatScalar *ba;
2655: MPI_Status stat;
2657: PetscFunctionBegin;
2658: PetscCheck(!idx, PETSC_COMM_SELF, PETSC_ERR_SUP, "Send email to petsc-maint@mcs.anl.gov");
2659: PetscCall(MatGetRowMaxAbs(a->A, v, NULL));
2660: PetscCall(VecGetArray(v, &va));
2662: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
2663: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
2665: bs = A->rmap->bs;
2666: mbs = a->mbs;
2667: Mbs = a->Mbs;
2668: ba = b->a;
2669: bi = b->i;
2670: bj = b->j;
2672: /* find ownerships */
2673: rowners_bs = A->rmap->range;
2675: /* each proc creates an array to be distributed */
2676: PetscCall(PetscCalloc1(bs * Mbs, &work));
2678: /* row_max for B */
2679: if (rank != size - 1) {
2680: for (i = 0; i < mbs; i++) {
2681: ncols = bi[1] - bi[0];
2682: bi++;
2683: brow = bs * i;
2684: for (j = 0; j < ncols; j++) {
2685: bcol = bs * (*bj);
2686: for (kcol = 0; kcol < bs; kcol++) {
2687: col = bcol + kcol; /* local col index */
2688: col += rowners_bs[rank + 1]; /* global col index */
2689: for (krow = 0; krow < bs; krow++) {
2690: atmp = PetscAbsScalar(*ba);
2691: ba++;
2692: row = brow + krow; /* local row index */
2693: if (PetscRealPart(va[row]) < atmp) va[row] = atmp;
2694: if (work[col] < atmp) work[col] = atmp;
2695: }
2696: }
2697: bj++;
2698: }
2699: }
2701: /* send values to its owners */
2702: for (PetscMPIInt dest = rank + 1; dest < size; dest++) {
2703: svalues = work + rowners_bs[dest];
2704: count = rowners_bs[dest + 1] - rowners_bs[dest];
2705: PetscCallMPI(MPIU_Send(svalues, count, MPIU_REAL, dest, rank, PetscObjectComm((PetscObject)A)));
2706: }
2707: }
2709: /* receive values */
2710: if (rank) {
2711: rvalues = work;
2712: count = rowners_bs[rank + 1] - rowners_bs[rank];
2713: for (source = 0; source < rank; source++) {
2714: PetscCallMPI(MPIU_Recv(rvalues, count, MPIU_REAL, MPI_ANY_SOURCE, MPI_ANY_TAG, PetscObjectComm((PetscObject)A), &stat));
2715: /* process values */
2716: for (i = 0; i < count; i++) {
2717: if (PetscRealPart(va[i]) < rvalues[i]) va[i] = rvalues[i];
2718: }
2719: }
2720: }
2722: PetscCall(VecRestoreArray(v, &va));
2723: PetscCall(PetscFree(work));
2724: PetscFunctionReturn(PETSC_SUCCESS);
2725: }
2727: static PetscErrorCode MatSOR_MPISBAIJ(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
2728: {
2729: Mat_MPISBAIJ *mat = (Mat_MPISBAIJ *)matin->data;
2730: PetscInt mbs = mat->mbs, bs = matin->rmap->bs;
2731: PetscScalar *x, *ptr, *from;
2732: Vec bb1;
2733: const PetscScalar *b;
2735: PetscFunctionBegin;
2736: PetscCheck(its > 0 && lits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Relaxation requires global its %" PetscInt_FMT " and local its %" PetscInt_FMT " both positive", its, lits);
2737: PetscCheck(bs <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "SSOR for block size > 1 is not yet implemented");
2739: if (flag == SOR_APPLY_UPPER) {
2740: PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
2741: PetscFunctionReturn(PETSC_SUCCESS);
2742: }
2744: if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
2745: if (flag & SOR_ZERO_INITIAL_GUESS) {
2746: PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, lits, xx));
2747: its--;
2748: }
2750: PetscCall(VecDuplicate(bb, &bb1));
2751: while (its--) {
2752: /* lower triangular part: slvec0b = - B^T*xx */
2753: PetscCall((*mat->B->ops->multtranspose)(mat->B, xx, mat->slvec0b));
2755: /* copy xx into slvec0a */
2756: PetscCall(VecGetArray(mat->slvec0, &ptr));
2757: PetscCall(VecGetArray(xx, &x));
2758: PetscCall(PetscArraycpy(ptr, x, bs * mbs));
2759: PetscCall(VecRestoreArray(mat->slvec0, &ptr));
2761: PetscCall(VecScale(mat->slvec0, -1.0));
2763: /* copy bb into slvec1a */
2764: PetscCall(VecGetArray(mat->slvec1, &ptr));
2765: PetscCall(VecGetArrayRead(bb, &b));
2766: PetscCall(PetscArraycpy(ptr, b, bs * mbs));
2767: PetscCall(VecRestoreArray(mat->slvec1, &ptr));
2769: /* set slvec1b = 0 */
2770: PetscCall(PetscObjectStateIncrease((PetscObject)mat->slvec1b));
2771: PetscCall(VecZeroEntries(mat->slvec1b));
2773: PetscCall(VecScatterBegin(mat->sMvctx, mat->slvec0, mat->slvec1, ADD_VALUES, SCATTER_FORWARD));
2774: PetscCall(VecRestoreArray(xx, &x));
2775: PetscCall(VecRestoreArrayRead(bb, &b));
2776: PetscCall(VecScatterEnd(mat->sMvctx, mat->slvec0, mat->slvec1, ADD_VALUES, SCATTER_FORWARD));
2778: /* upper triangular part: bb1 = bb1 - B*x */
2779: PetscCall((*mat->B->ops->multadd)(mat->B, mat->slvec1b, mat->slvec1a, bb1));
2781: /* local diagonal sweep */
2782: PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, lits, xx));
2783: }
2784: PetscCall(VecDestroy(&bb1));
2785: } else if ((flag & SOR_LOCAL_FORWARD_SWEEP) && (its == 1) && (flag & SOR_ZERO_INITIAL_GUESS)) {
2786: PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
2787: } else if ((flag & SOR_LOCAL_BACKWARD_SWEEP) && (its == 1) && (flag & SOR_ZERO_INITIAL_GUESS)) {
2788: PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
2789: } else if (flag & SOR_EISENSTAT) {
2790: Vec xx1;
2791: PetscBool hasop;
2792: const PetscScalar *diag;
2793: PetscScalar *sl, scale = (omega - 2.0) / omega;
2794: PetscInt i, n;
2796: if (!mat->xx1) {
2797: PetscCall(VecDuplicate(bb, &mat->xx1));
2798: PetscCall(VecDuplicate(bb, &mat->bb1));
2799: }
2800: xx1 = mat->xx1;
2801: bb1 = mat->bb1;
2803: PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP), fshift, lits, 1, xx));
2805: if (!mat->diag) {
2806: /* this is wrong for same matrix with new nonzero values */
2807: PetscCall(MatCreateVecs(matin, &mat->diag, NULL));
2808: PetscCall(MatGetDiagonal(matin, mat->diag));
2809: }
2810: PetscCall(MatHasOperation(matin, MATOP_MULT_DIAGONAL_BLOCK, &hasop));
2812: if (hasop) {
2813: PetscCall(MatMultDiagonalBlock(matin, xx, bb1));
2814: PetscCall(VecAYPX(mat->slvec1a, scale, bb));
2815: } else {
2816: /*
2817: These two lines are replaced by code that may be a bit faster for a good compiler
2818: PetscCall(VecPointwiseMult(mat->slvec1a,mat->diag,xx));
2819: PetscCall(VecAYPX(mat->slvec1a,scale,bb));
2820: */
2821: PetscCall(VecGetArray(mat->slvec1a, &sl));
2822: PetscCall(VecGetArrayRead(mat->diag, &diag));
2823: PetscCall(VecGetArrayRead(bb, &b));
2824: PetscCall(VecGetArray(xx, &x));
2825: PetscCall(VecGetLocalSize(xx, &n));
2826: if (omega == 1.0) {
2827: for (i = 0; i < n; i++) sl[i] = b[i] - diag[i] * x[i];
2828: PetscCall(PetscLogFlops(2.0 * n));
2829: } else {
2830: for (i = 0; i < n; i++) sl[i] = b[i] + scale * diag[i] * x[i];
2831: PetscCall(PetscLogFlops(3.0 * n));
2832: }
2833: PetscCall(VecRestoreArray(mat->slvec1a, &sl));
2834: PetscCall(VecRestoreArrayRead(mat->diag, &diag));
2835: PetscCall(VecRestoreArrayRead(bb, &b));
2836: PetscCall(VecRestoreArray(xx, &x));
2837: }
2839: /* multiply off-diagonal portion of matrix */
2840: PetscCall(PetscObjectStateIncrease((PetscObject)mat->slvec1b));
2841: PetscCall(VecZeroEntries(mat->slvec1b));
2842: PetscCall((*mat->B->ops->multtranspose)(mat->B, xx, mat->slvec0b));
2843: PetscCall(VecGetArray(mat->slvec0, &from));
2844: PetscCall(VecGetArray(xx, &x));
2845: PetscCall(PetscArraycpy(from, x, bs * mbs));
2846: PetscCall(VecRestoreArray(mat->slvec0, &from));
2847: PetscCall(VecRestoreArray(xx, &x));
2848: PetscCall(VecScatterBegin(mat->sMvctx, mat->slvec0, mat->slvec1, ADD_VALUES, SCATTER_FORWARD));
2849: PetscCall(VecScatterEnd(mat->sMvctx, mat->slvec0, mat->slvec1, ADD_VALUES, SCATTER_FORWARD));
2850: PetscCall((*mat->B->ops->multadd)(mat->B, mat->slvec1b, mat->slvec1a, mat->slvec1a));
2852: /* local sweep */
2853: PetscCall((*mat->A->ops->sor)(mat->A, mat->slvec1a, omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP), fshift, lits, 1, xx1));
2854: PetscCall(VecAXPY(xx, 1.0, xx1));
2855: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatSORType is not supported for SBAIJ matrix format");
2856: PetscFunctionReturn(PETSC_SUCCESS);
2857: }
2859: /*@
2860: MatCreateMPISBAIJWithArrays - creates a `MATMPISBAIJ` matrix using arrays that contain in standard CSR format for the local rows.
2862: Collective
2864: Input Parameters:
2865: + comm - MPI communicator
2866: . bs - the block size, only a block size of 1 is supported
2867: . m - number of local rows (Cannot be `PETSC_DECIDE`)
2868: . n - This value should be the same as the local size used in creating the
2869: x vector for the matrix-vector product $ y = Ax $. (or `PETSC_DECIDE` to have
2870: calculated if `N` is given) For square matrices `n` is almost always `m`.
2871: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
2872: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
2873: . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that row block row of the matrix
2874: . j - column indices
2875: - a - matrix values
2877: Output Parameter:
2878: . mat - the matrix
2880: Level: intermediate
2882: Notes:
2883: The `i`, `j`, and `a` arrays ARE copied by this routine into the internal format used by PETSc;
2884: thus you CANNOT change the matrix entries by changing the values of `a` after you have
2885: called this routine. Use `MatCreateMPIAIJWithSplitArrays()` to avoid needing to copy the arrays.
2887: The `i` and `j` indices are 0 based, and `i` indices are indices corresponding to the local `j` array.
2889: .seealso: [](ch_matrices), `Mat`, `MATMPISBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatMPIAIJSetPreallocation()`, `MatMPIAIJSetPreallocationCSR()`,
2890: `MATMPIAIJ`, `MatCreateAIJ()`, `MatCreateMPIAIJWithSplitArrays()`, `MatMPISBAIJSetPreallocationCSR()`
2891: @*/
2892: PetscErrorCode MatCreateMPISBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, const PetscInt i[], const PetscInt j[], const PetscScalar a[], Mat *mat)
2893: {
2894: PetscFunctionBegin;
2895: PetscCheck(!i[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");
2896: PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "local number of rows (m) cannot be PETSC_DECIDE, or negative");
2897: PetscCall(MatCreate(comm, mat));
2898: PetscCall(MatSetSizes(*mat, m, n, M, N));
2899: PetscCall(MatSetType(*mat, MATMPISBAIJ));
2900: PetscCall(MatMPISBAIJSetPreallocationCSR(*mat, bs, i, j, a));
2901: PetscFunctionReturn(PETSC_SUCCESS);
2902: }
2904: /*@
2905: MatMPISBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATMPISBAIJ` format using the given nonzero structure and (optional) numerical values
2907: Collective
2909: Input Parameters:
2910: + B - the matrix
2911: . bs - the block size
2912: . i - the indices into `j` for the start of each local row (indices start with zero)
2913: . j - the column indices for each local row (indices start with zero) these must be sorted for each row
2914: - v - optional values in the matrix, pass `NULL` if not provided
2916: Level: advanced
2918: Notes:
2919: The `i`, `j`, and `v` arrays ARE copied by this routine into the internal format used by PETSc;
2920: thus you CANNOT change the matrix entries by changing the values of `v` after you have
2921: called this routine.
2923: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries
2924: and usually the numerical values as well
2926: Any entries passed in that are below the diagonal are ignored
2928: .seealso: [](ch_matrices), `Mat`, `MATMPISBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatMPIBAIJSetPreallocation()`, `MatCreateAIJ()`, `MATMPIAIJ`,
2929: `MatCreateMPISBAIJWithArrays()`
2930: @*/
2931: PetscErrorCode MatMPISBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[])
2932: {
2933: PetscFunctionBegin;
2934: PetscTryMethod(B, "MatMPISBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v));
2935: PetscFunctionReturn(PETSC_SUCCESS);
2936: }
2938: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
2939: {
2940: PetscInt m, N, i, rstart, nnz, Ii, bs, cbs;
2941: PetscInt *indx;
2942: PetscScalar *values;
2944: PetscFunctionBegin;
2945: PetscCall(MatGetSize(inmat, &m, &N));
2946: if (scall == MAT_INITIAL_MATRIX) { /* symbolic phase */
2947: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)inmat->data;
2948: PetscInt *dnz, *onz, mbs, Nbs, nbs;
2949: PetscInt *bindx, rmax = a->rmax, j;
2950: PetscMPIInt rank, size;
2952: PetscCall(MatGetBlockSizes(inmat, &bs, &cbs));
2953: mbs = m / bs;
2954: Nbs = N / cbs;
2955: if (n == PETSC_DECIDE) PetscCall(PetscSplitOwnershipBlock(comm, cbs, &n, &N));
2956: nbs = n / cbs;
2958: PetscCall(PetscMalloc1(rmax, &bindx));
2959: MatPreallocateBegin(comm, mbs, nbs, dnz, onz); /* inline function, output __end and __rstart are used below */
2961: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2962: PetscCallMPI(MPI_Comm_size(comm, &size));
2963: if (rank == size - 1) {
2964: /* Check sum(nbs) = Nbs */
2965: PetscCheck(__end == Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Sum of local block columns %" PetscInt_FMT " != global block columns %" PetscInt_FMT, __end, Nbs);
2966: }
2968: rstart = __rstart; /* block rstart of *outmat; see inline function MatPreallocateBegin */
2969: PetscCall(MatSetOption(inmat, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE));
2970: for (i = 0; i < mbs; i++) {
2971: PetscCall(MatGetRow_SeqSBAIJ(inmat, i * bs, &nnz, &indx, NULL)); /* non-blocked nnz and indx */
2972: nnz = nnz / bs;
2973: for (j = 0; j < nnz; j++) bindx[j] = indx[j * bs] / bs;
2974: PetscCall(MatPreallocateSet(i + rstart, nnz, bindx, dnz, onz));
2975: PetscCall(MatRestoreRow_SeqSBAIJ(inmat, i * bs, &nnz, &indx, NULL));
2976: }
2977: PetscCall(MatSetOption(inmat, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE));
2978: PetscCall(PetscFree(bindx));
2980: PetscCall(MatCreate(comm, outmat));
2981: PetscCall(MatSetSizes(*outmat, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
2982: PetscCall(MatSetBlockSizes(*outmat, bs, cbs));
2983: PetscCall(MatSetType(*outmat, MATSBAIJ));
2984: PetscCall(MatSeqSBAIJSetPreallocation(*outmat, bs, 0, dnz));
2985: PetscCall(MatMPISBAIJSetPreallocation(*outmat, bs, 0, dnz, 0, onz));
2986: MatPreallocateEnd(dnz, onz);
2987: }
2989: /* numeric phase */
2990: PetscCall(MatGetBlockSizes(inmat, &bs, &cbs));
2991: PetscCall(MatGetOwnershipRange(*outmat, &rstart, NULL));
2993: PetscCall(MatSetOption(inmat, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE));
2994: for (i = 0; i < m; i++) {
2995: PetscCall(MatGetRow_SeqSBAIJ(inmat, i, &nnz, &indx, &values));
2996: Ii = i + rstart;
2997: PetscCall(MatSetValues(*outmat, 1, &Ii, nnz, indx, values, INSERT_VALUES));
2998: PetscCall(MatRestoreRow_SeqSBAIJ(inmat, i, &nnz, &indx, &values));
2999: }
3000: PetscCall(MatSetOption(inmat, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE));
3001: PetscCall(MatAssemblyBegin(*outmat, MAT_FINAL_ASSEMBLY));
3002: PetscCall(MatAssemblyEnd(*outmat, MAT_FINAL_ASSEMBLY));
3003: PetscFunctionReturn(PETSC_SUCCESS);
3004: }