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