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