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