Actual source code: sbaij.c
1: /*
2: Defines the basic matrix operations for the SBAIJ (compressed row)
3: matrix storage format.
4: */
5: #include <../src/mat/impls/baij/seq/baij.h>
6: #include <../src/mat/impls/sbaij/seq/sbaij.h>
7: #include <petsc/private/kernels/blocktranspose.h>
8: #include <petscblaslapack.h>
10: #include <../src/mat/impls/sbaij/seq/relax.h>
11: #define USESHORT
12: #include <../src/mat/impls/sbaij/seq/relax.h>
14: /* defines MatSetValues_Seq_Hash(), MatAssemblyEnd_Seq_Hash(), MatSetUp_Seq_Hash() */
15: #define TYPE SBAIJ
16: #define TYPE_SBAIJ
17: #define TYPE_BS
18: #include "../src/mat/impls/aij/seq/seqhashmatsetvalues.h"
19: #undef TYPE_BS
20: #define TYPE_BS _BS
21: #define TYPE_BS_ON
22: #include "../src/mat/impls/aij/seq/seqhashmatsetvalues.h"
23: #undef TYPE_BS
24: #undef TYPE_SBAIJ
25: #include "../src/mat/impls/aij/seq/seqhashmat.h"
26: #undef TYPE
27: #undef TYPE_BS_ON
29: #if PetscDefined(HAVE_ELEMENTAL)
30: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat, MatType, MatReuse, Mat *);
31: #endif
32: #if PetscDefined(HAVE_SCALAPACK) && (PetscDefined(USE_REAL_SINGLE) || PetscDefined(USE_REAL_DOUBLE))
33: PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *);
34: #endif
35: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat, MatType, MatReuse, Mat *);
37: MatGetDiagonalMarkers(SeqSBAIJ, A->rmap->bs)
39: static PetscErrorCode MatGetColumnReductions_SeqSBAIJ(Mat A, PetscInt type, PetscReal *reductions)
40: {
41: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)A->data;
42: PetscInt m, n, bs = A->rmap->bs, bs2 = aij->bs2;
43: PetscBool implicit, hermitian;
44: const PetscInt *ai = aij->i, *aj = aij->j;
45: const MatScalar *aa = aij->a;
47: PetscFunctionBegin;
48: PetscCall(MatGetSize(A, &m, &n));
49: PetscCheck(type == NORM_2 || type == NORM_1 || type == NORM_INFINITY || type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART || type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
50: PetscCall(PetscArrayzero(reductions, n));
51: implicit = (PetscBool)(m == n && (A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE));
52: hermitian = (PetscBool)(implicit && PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE);
53: for (PetscInt i = 0; i < aij->mbs; i++) {
54: for (PetscInt k = ai[i]; k < ai[i + 1]; k++) {
55: const PetscBool offdiag = (PetscBool)(implicit && aj[k] != i);
57: for (PetscInt jb = 0; jb < bs; jb++) {
58: for (PetscInt ib = 0; ib < bs; ib++) {
59: const MatScalar value = aa[k * bs2 + jb * bs + ib];
60: const PetscInt col = aj[k] * bs + jb, row = i * bs + ib;
61: PetscReal reduction;
63: if (type == NORM_2) reduction = PetscAbsScalar(value * value);
64: else if (type == NORM_1 || type == NORM_INFINITY) reduction = PetscAbsScalar(value);
65: else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) reduction = PetscRealPart(value);
66: else reduction = PetscImaginaryPart(value);
67: if (type == NORM_INFINITY) reductions[col] = PetscMax(reduction, reductions[col]);
68: else reductions[col] += reduction;
69: if (offdiag) {
70: if (hermitian && (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART)) reduction = -reduction;
71: if (type == NORM_INFINITY) reductions[row] = PetscMax(reduction, reductions[row]);
72: else reductions[row] += reduction;
73: }
74: }
75: }
76: }
77: }
78: if (type == NORM_2) {
79: for (PetscInt i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
80: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
81: for (PetscInt i = 0; i < n; i++) reductions[i] /= m;
82: }
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done)
87: {
88: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
89: PetscInt i, j, n = a->mbs, nz = a->i[n], *tia, *tja, bs = A->rmap->bs, k, l, cnt;
90: PetscInt **ia = (PetscInt **)inia, **ja = (PetscInt **)inja;
92: PetscFunctionBegin;
93: *nn = n;
94: if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
95: if (symmetric) {
96: PetscCall(MatToSymmetricIJ_SeqAIJ(n, a->i, a->j, PETSC_FALSE, 0, 0, &tia, &tja));
97: nz = tia[n];
98: } else {
99: tia = a->i;
100: tja = a->j;
101: }
103: if (!blockcompressed && bs > 1) {
104: (*nn) *= bs;
105: /* malloc & create the natural set of indices */
106: PetscCall(PetscMalloc1((n + 1) * bs, ia));
107: if (n) {
108: (*ia)[0] = oshift;
109: for (j = 1; j < bs; j++) (*ia)[j] = (tia[1] - tia[0]) * bs + (*ia)[j - 1];
110: }
112: for (i = 1; i < n; i++) {
113: (*ia)[i * bs] = (tia[i] - tia[i - 1]) * bs + (*ia)[i * bs - 1];
114: for (j = 1; j < bs; j++) (*ia)[i * bs + j] = (tia[i + 1] - tia[i]) * bs + (*ia)[i * bs + j - 1];
115: }
116: if (n) (*ia)[n * bs] = (tia[n] - tia[n - 1]) * bs + (*ia)[n * bs - 1];
118: if (inja) {
119: PetscCall(PetscMalloc1(nz * bs * bs, ja));
120: cnt = 0;
121: for (i = 0; i < n; i++) {
122: for (j = 0; j < bs; j++) {
123: for (k = tia[i]; k < tia[i + 1]; k++) {
124: for (l = 0; l < bs; l++) (*ja)[cnt++] = bs * tja[k] + l;
125: }
126: }
127: }
128: }
130: if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
131: PetscCall(PetscFree(tia));
132: PetscCall(PetscFree(tja));
133: }
134: } else if (oshift == 1) {
135: if (symmetric) {
136: nz = tia[A->rmap->n / bs];
137: /* add 1 to i and j indices */
138: for (i = 0; i < A->rmap->n / bs + 1; i++) tia[i] = tia[i] + 1;
139: *ia = tia;
140: if (ja) {
141: for (i = 0; i < nz; i++) tja[i] = tja[i] + 1;
142: *ja = tja;
143: }
144: } else {
145: nz = a->i[A->rmap->n / bs];
146: /* malloc space and add 1 to i and j indices */
147: PetscCall(PetscMalloc1(A->rmap->n / bs + 1, ia));
148: for (i = 0; i < A->rmap->n / bs + 1; i++) (*ia)[i] = a->i[i] + 1;
149: if (ja) {
150: PetscCall(PetscMalloc1(nz, ja));
151: for (i = 0; i < nz; i++) (*ja)[i] = a->j[i] + 1;
152: }
153: }
154: } else {
155: *ia = tia;
156: if (ja) *ja = tja;
157: }
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done)
162: {
163: PetscFunctionBegin;
164: if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
165: if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
166: PetscCall(PetscFree(*ia));
167: if (ja) PetscCall(PetscFree(*ja));
168: }
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: PetscErrorCode MatDestroy_SeqSBAIJ(Mat A)
173: {
174: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
176: PetscFunctionBegin;
177: if (A->hash_active) {
178: PetscInt bs;
179: A->ops[0] = a->cops;
180: PetscCall(PetscHMapIJVDestroy(&a->ht));
181: PetscCall(MatGetBlockSize(A, &bs));
182: if (bs > 1) PetscCall(PetscHSetIJDestroy(&a->bht));
183: PetscCall(PetscFree(a->dnz));
184: PetscCall(PetscFree(a->bdnz));
185: A->hash_active = PETSC_FALSE;
186: }
187: PetscCall(PetscLogObjectState((PetscObject)A, "Rows=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, A->rmap->N, a->nz));
188: PetscCall(MatSeqXAIJFreeAIJ(A, &a->a, &a->j, &a->i));
189: PetscCall(PetscFree(a->diag));
190: PetscCall(ISDestroy(&a->row));
191: PetscCall(ISDestroy(&a->col));
192: PetscCall(ISDestroy(&a->icol));
193: PetscCall(PetscFree(a->idiag));
194: PetscCall(PetscFree(a->inode.size_csr));
195: if (a->free_imax_ilen) PetscCall(PetscFree2(a->imax, a->ilen));
196: PetscCall(PetscFree(a->solve_work));
197: PetscCall(PetscFree(a->sor_work));
198: PetscCall(PetscFree(a->solves_work));
199: PetscCall(PetscFree(a->mult_work));
200: PetscCall(PetscFree(a->saved_values));
201: if (a->free_jshort) PetscCall(PetscFree(a->jshort));
202: PetscCall(PetscFree(a->inew));
203: PetscCall(MatDestroy(&a->parent));
204: PetscCall(PetscFree(A->data));
206: PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
207: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJGetArray_C", NULL));
208: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJRestoreArray_C", NULL));
209: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatStoreValues_C", NULL));
210: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatRetrieveValues_C", NULL));
211: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetColumnIndices_C", NULL));
212: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_seqaij_C", NULL));
213: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_seqbaij_C", NULL));
214: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetPreallocation_C", NULL));
215: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetPreallocationCSR_C", NULL));
216: #if PetscDefined(HAVE_ELEMENTAL)
217: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_elemental_C", NULL));
218: #endif
219: #if PetscDefined(HAVE_SCALAPACK) && (PetscDefined(USE_REAL_SINGLE) || PetscDefined(USE_REAL_DOUBLE))
220: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_scalapack_C", NULL));
221: #endif
222: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode MatSetOption_SeqSBAIJ(Mat A, MatOption op, PetscBool flg)
227: {
228: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
230: PetscFunctionBegin;
231: switch (op) {
232: case MAT_ROW_ORIENTED:
233: a->roworiented = flg;
234: break;
235: case MAT_KEEP_NONZERO_PATTERN:
236: a->keepnonzeropattern = flg;
237: break;
238: case MAT_NEW_NONZERO_LOCATIONS:
239: a->nonew = (flg ? 0 : 1);
240: break;
241: case MAT_NEW_NONZERO_LOCATION_ERR:
242: a->nonew = (flg ? -1 : 0);
243: break;
244: case MAT_NEW_NONZERO_ALLOCATION_ERR:
245: a->nonew = (flg ? -2 : 0);
246: break;
247: case MAT_UNUSED_NONZERO_LOCATION_ERR:
248: a->nounused = (flg ? -1 : 0);
249: break;
250: case MAT_HERMITIAN:
251: if (PetscDefined(USE_COMPLEX) && flg) { /* disable transpose ops */
252: PetscInt bs;
254: PetscCall(MatGetBlockSize(A, &bs));
255: PetscCheck(bs <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for Hermitian with block size greater than 1");
256: A->ops->multtranspose = NULL;
257: A->ops->multtransposeadd = NULL;
258: }
259: break;
260: case MAT_SYMMETRIC:
261: case MAT_SPD:
262: if (PetscDefined(USE_COMPLEX) && flg) { /* An Hermitian and symmetric matrix has zero imaginary part (restore back transpose ops) */
263: A->ops->multtranspose = A->ops->mult;
264: A->ops->multtransposeadd = A->ops->multadd;
265: }
266: break;
267: case MAT_IGNORE_LOWER_TRIANGULAR:
268: case MAT_ERROR_LOWER_TRIANGULAR:
269: a->ignore_ltriangular = flg;
270: break;
271: case MAT_GETROW_UPPERTRIANGULAR:
272: a->getrow_utriangular = flg;
273: break;
274: case MAT_STRUCTURE_ONLY:
275: if (flg) {
276: PetscCall(MatXAIJDeallocatea(A, &a->a));
277: a->a = NULL;
278: }
279: break;
280: default:
281: break;
282: }
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: PetscErrorCode MatGetRow_SeqSBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
287: {
288: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
290: PetscFunctionBegin;
291: PetscCheck(!A || a->getrow_utriangular, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE) or MatGetRowUpperTriangular()");
293: /* Get the upper triangular part of the row */
294: PetscCall(MatGetRow_SeqBAIJ_private(A, row, nz, idx, v, a->i, a->j, a->a));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
299: {
300: PetscFunctionBegin;
301: if (idx) PetscCall(PetscFree(*idx));
302: if (v) PetscCall(PetscFree(*v));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: static PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A)
307: {
308: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
310: PetscFunctionBegin;
311: a->getrow_utriangular = PETSC_TRUE;
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: static PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A)
316: {
317: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
319: PetscFunctionBegin;
320: a->getrow_utriangular = PETSC_FALSE;
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: static PetscErrorCode MatTranspose_SeqSBAIJ(Mat A, MatReuse reuse, Mat *B)
325: {
326: PetscFunctionBegin;
327: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
328: if (reuse == MAT_INITIAL_MATRIX) {
329: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, B));
330: } else if (reuse == MAT_REUSE_MATRIX) {
331: PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN));
332: }
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: static PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A, PetscViewer viewer)
337: {
338: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
339: PetscInt i, j, bs = A->rmap->bs, k, l, bs2 = a->bs2;
340: PetscViewerFormat format;
341: const PetscInt *diag;
342: const char *matname;
344: PetscFunctionBegin;
345: if (A->structure_only) {
346: PetscCall(MatView_SeqBAIJ(A, viewer));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
349: PetscCall(PetscViewerGetFormat(viewer, &format));
350: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
351: } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
352: Mat aij;
354: if (A->factortype && bs > 1) {
355: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n"));
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
358: PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &aij));
359: if (((PetscObject)A)->name) PetscCall(PetscObjectGetName((PetscObject)A, &matname));
360: if (((PetscObject)A)->name) PetscCall(PetscObjectSetName((PetscObject)aij, matname));
361: PetscCall(MatView_SeqAIJ(aij, viewer));
362: PetscCall(MatDestroy(&aij));
363: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
364: Mat B;
366: PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &B));
367: if (((PetscObject)A)->name) PetscCall(PetscObjectGetName((PetscObject)A, &matname));
368: if (((PetscObject)A)->name) PetscCall(PetscObjectSetName((PetscObject)B, matname));
369: PetscCall(MatView_SeqAIJ(B, viewer));
370: PetscCall(MatDestroy(&B));
371: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
372: PetscFunctionReturn(PETSC_SUCCESS);
373: } else {
374: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
375: if (A->factortype) { /* for factored matrix */
376: PetscCheck(bs <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "matrix is factored with bs>1. Not implemented yet");
377: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, &diag, NULL));
378: for (i = 0; i < a->mbs; i++) { /* for row block i */
379: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
380: /* diagonal entry */
381: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[diag[i]]) > 0.0) {
382: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]]), (double)PetscImaginaryPart(1.0 / a->a[diag[i]])));
383: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[diag[i]]) < 0.0) {
384: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]]), -(double)PetscImaginaryPart(1.0 / a->a[diag[i]])));
385: } else {
386: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]])));
387: }
388: /* off-diagonal entries */
389: for (k = a->i[i]; k < a->i[i + 1] - 1; k++) {
390: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[k]) > 0.0) {
391: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", bs * a->j[k], (double)PetscRealPart(a->a[k]), (double)PetscImaginaryPart(a->a[k])));
392: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[k]) < 0.0) {
393: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", bs * a->j[k], (double)PetscRealPart(a->a[k]), -(double)PetscImaginaryPart(a->a[k])));
394: } else {
395: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k], (double)PetscRealPart(a->a[k])));
396: }
397: }
398: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
399: }
401: } else { /* for non-factored matrix */
402: for (i = 0; i < a->mbs; i++) { /* for row block i */
403: for (j = 0; j < bs; j++) { /* for row bs*i + j */
404: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j));
405: for (k = a->i[i]; k < a->i[i + 1]; k++) { /* for column block */
406: for (l = 0; l < bs; l++) { /* for column */
407: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0) {
408: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), (double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j])));
409: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0) {
410: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), -(double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j])));
411: } else {
412: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j])));
413: }
414: }
415: }
416: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
417: }
418: }
419: }
420: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
421: }
422: PetscCall(PetscViewerFlush(viewer));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: #include <petscdraw.h>
427: static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw, void *Aa)
428: {
429: Mat A = (Mat)Aa;
430: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
431: PetscInt row, i, j, k, l, mbs = a->mbs, bs = A->rmap->bs, bs2 = a->bs2;
432: PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r;
433: MatScalar *aa;
434: PetscViewer viewer;
435: int color;
437: PetscFunctionBegin;
438: PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
439: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
441: /* loop over matrix elements drawing boxes */
443: PetscDrawCollectiveBegin(draw);
444: PetscCall(PetscDrawString(draw, .3 * (xl + xr), .3 * (yl + yr), PETSC_DRAW_BLACK, "symmetric"));
445: /* Blue for negative, Cyan for zero and Red for positive */
446: color = PETSC_DRAW_BLUE;
447: for (i = 0, row = 0; i < mbs; i++, row += bs) {
448: for (j = a->i[i]; j < a->i[i + 1]; j++) {
449: y_l = A->rmap->N - row - 1.0;
450: y_r = y_l + 1.0;
451: x_l = a->j[j] * bs;
452: x_r = x_l + 1.0;
453: aa = a->a + j * bs2;
454: for (k = 0; k < bs; k++) {
455: for (l = 0; l < bs; l++) {
456: if (PetscRealPart(*aa++) >= 0.) continue;
457: PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
458: }
459: }
460: }
461: }
462: color = PETSC_DRAW_CYAN;
463: for (i = 0, row = 0; i < mbs; i++, row += bs) {
464: for (j = a->i[i]; j < a->i[i + 1]; j++) {
465: y_l = A->rmap->N - row - 1.0;
466: y_r = y_l + 1.0;
467: x_l = a->j[j] * bs;
468: x_r = x_l + 1.0;
469: aa = a->a + j * bs2;
470: for (k = 0; k < bs; k++) {
471: for (l = 0; l < bs; l++) {
472: if (PetscRealPart(*aa++) != 0.) continue;
473: PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
474: }
475: }
476: }
477: }
478: color = PETSC_DRAW_RED;
479: for (i = 0, row = 0; i < mbs; i++, row += bs) {
480: for (j = a->i[i]; j < a->i[i + 1]; j++) {
481: y_l = A->rmap->N - row - 1.0;
482: y_r = y_l + 1.0;
483: x_l = a->j[j] * bs;
484: x_r = x_l + 1.0;
485: aa = a->a + j * bs2;
486: for (k = 0; k < bs; k++) {
487: for (l = 0; l < bs; l++) {
488: if (PetscRealPart(*aa++) <= 0.) continue;
489: PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
490: }
491: }
492: }
493: }
494: PetscDrawCollectiveEnd(draw);
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A, PetscViewer viewer)
499: {
500: PetscReal xl, yl, xr, yr, w, h;
501: PetscDraw draw;
502: PetscBool isnull;
504: PetscFunctionBegin;
505: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
506: PetscCall(PetscDrawIsNull(draw, &isnull));
507: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
509: xr = A->rmap->N;
510: yr = A->rmap->N;
511: h = yr / 10.0;
512: w = xr / 10.0;
513: xr += w;
514: yr += h;
515: xl = -w;
516: yl = -h;
517: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
518: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
519: PetscCall(PetscDrawZoom(draw, MatView_SeqSBAIJ_Draw_Zoom, A));
520: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
521: PetscCall(PetscDrawSave(draw));
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: /* Used for both MPIBAIJ and MPISBAIJ matrices */
526: #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary
528: PetscErrorCode MatView_SeqSBAIJ(Mat A, PetscViewer viewer)
529: {
530: PetscBool isascii, isbinary, isdraw;
532: PetscFunctionBegin;
533: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
534: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
535: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
536: if (isascii) {
537: PetscCall(MatView_SeqSBAIJ_ASCII(A, viewer));
538: } else if (isbinary) {
539: PetscCall(MatView_SeqSBAIJ_Binary(A, viewer));
540: } else if (isdraw) {
541: PetscCall(MatView_SeqSBAIJ_Draw(A, viewer));
542: } else {
543: Mat B;
544: const char *matname;
545: PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &B));
546: if (((PetscObject)A)->name) PetscCall(PetscObjectGetName((PetscObject)A, &matname));
547: if (((PetscObject)A)->name) PetscCall(PetscObjectSetName((PetscObject)B, matname));
548: PetscCall(MatView(B, viewer));
549: PetscCall(MatDestroy(&B));
550: }
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: PetscErrorCode MatGetValues_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], PetscScalar v[])
555: {
556: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
557: PetscInt *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j;
558: PetscInt *ai = a->i, *ailen = a->ilen;
559: PetscInt brow, bcol, ridx, cidx, bs = A->rmap->bs, bs2 = a->bs2;
560: MatScalar *ap, *aa = a->a;
561: PetscBool roworiented = a->roworiented;
562: PetscScalar *value;
564: PetscFunctionBegin;
565: for (k = 0; k < m; k++) { /* loop over rows */
566: row = im[k];
567: if (row < 0) continue; /* negative row */
568: PetscCheck(row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, A->rmap->N - 1);
569: brow = row / bs;
570: rp = aj + ai[brow];
571: ap = aa + bs2 * ai[brow];
572: nrow = ailen[brow];
573: for (l = 0; l < n; l++) { /* loop over columns */
574: if (in[l] < 0) continue; /* negative column */
575: PetscCheck(in[l] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[l], A->cmap->n - 1);
576: value = roworiented ? &v[l + k * n] : &v[k + l * m];
577: col = in[l];
578: bcol = col / bs;
579: cidx = col % bs;
580: ridx = row % bs;
581: high = nrow;
582: low = 0; /* assume unsorted */
583: while (high - low > 5) {
584: t = (low + high) / 2;
585: if (rp[t] > bcol) high = t;
586: else low = t;
587: }
588: for (i = low; i < high; i++) {
589: if (rp[i] > bcol) break;
590: if (rp[i] == bcol) {
591: *value = ap[bs2 * i + bs * cidx + ridx];
592: goto finished;
593: }
594: }
595: *value = 0.0;
596: finished:;
597: }
598: }
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: static PetscErrorCode MatPermute_SeqSBAIJ(Mat A, IS rowp, IS colp, Mat *B)
603: {
604: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = NULL;
605: Mat_SeqBAIJ *c = NULL;
606: Mat C;
607: IS browp, bcolp;
608: const PetscInt *row, *col;
609: PetscInt *irow, *icol, *lens, *bi, *bj, *bilen;
610: MatScalar *ba, *work;
611: PetscInt mbs, nbs, bs = A->rmap->bs, bs2 = a->bs2;
612: PetscBool same = (PetscBool)(rowp == colp), implicit = (PetscBool)(A->rmap->N == A->cmap->N && (A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE)), hermitian;
614: PetscFunctionBegin;
615: if (same == PETSC_FALSE) PetscCall(ISEqualUnsorted(rowp, colp, &same));
616: hermitian = (PetscBool)(implicit == PETSC_TRUE && PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE);
617: same = (PetscBool)(implicit == PETSC_TRUE && same == PETSC_TRUE);
618: PetscCall(ISCompressIndicesGeneral(A->rmap->N, A->rmap->n, bs, 1, &rowp, &browp));
619: if (rowp == colp) {
620: bcolp = browp;
621: PetscCall(PetscObjectReference((PetscObject)bcolp));
622: } else PetscCall(ISCompressIndicesGeneral(A->cmap->N, A->cmap->n, bs, 1, &colp, &bcolp));
623: PetscCall(ISGetLocalSize(browp, &mbs));
624: PetscCall(ISGetLocalSize(bcolp, &nbs));
625: PetscCheck(mbs == a->mbs && nbs == a->nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Row and column index sets must be block permutations");
626: PetscCall(ISGetIndices(browp, &row));
627: PetscCall(ISGetIndices(bcolp, &col));
628: PetscCall(PetscMalloc2(mbs, &irow, nbs, &icol));
629: for (PetscInt i = 0; i < mbs; i++) {
630: PetscCheck(row[i] >= 0 && row[i] < mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Row index set is not a block permutation");
631: irow[row[i]] = i;
632: }
633: for (PetscInt i = 0; i < nbs; i++) {
634: PetscCheck(col[i] >= 0 && col[i] < nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Column index set is not a block permutation");
635: icol[col[i]] = i;
636: }
637: PetscCall(PetscCalloc1(mbs, &lens));
638: for (PetscInt i = 0; i < a->mbs; i++) {
639: for (PetscInt k = a->i[i]; k < a->i[i + 1]; k++) {
640: const PetscInt j = a->j[k];
642: if (same == PETSC_TRUE) lens[PetscMin(irow[i], icol[j])]++;
643: else {
644: lens[irow[i]]++;
645: if (implicit == PETSC_TRUE && i != j) lens[irow[j]]++;
646: }
647: }
648: }
649: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
650: PetscCall(MatSetSizes(C, mbs * bs, nbs * bs, mbs * bs, nbs * bs));
651: PetscCall(MatSetType(C, same == PETSC_TRUE ? ((PetscObject)A)->type_name : MATSEQBAIJ));
652: if (same == PETSC_TRUE) {
653: PetscCall(MatSeqSBAIJSetPreallocation(C, bs, 0, lens));
654: b = (Mat_SeqSBAIJ *)C->data;
655: bi = b->i;
656: bj = b->j;
657: ba = b->a;
658: bilen = b->ilen;
659: } else {
660: PetscCall(MatSeqBAIJSetPreallocation(C, bs, 0, lens));
661: c = (Mat_SeqBAIJ *)C->data;
662: bi = c->i;
663: bj = c->j;
664: ba = c->a;
665: bilen = c->ilen;
666: }
667: PetscCall(PetscFree(lens));
668: for (PetscInt i = 0; i < a->mbs; i++) {
669: for (PetscInt k = a->i[i]; k < a->i[i + 1]; k++) {
670: const PetscInt j = a->j[k];
671: PetscInt r = irow[i], colidx = icol[j], pos;
672: MatScalar *block;
674: if (same == PETSC_TRUE && r > colidx) {
675: pos = r;
676: r = colidx;
677: colidx = pos;
678: }
679: pos = bi[r] + bilen[r]++;
680: bj[pos] = colidx;
681: block = ba + pos * bs2;
682: PetscCall(PetscArraycpy(block, a->a + k * bs2, bs2));
683: if (same == PETSC_TRUE && irow[i] > icol[j]) {
684: PetscCall(PetscKernel_A_gets_transpose_A_N(block, bs));
685: if (hermitian == PETSC_TRUE) {
686: for (PetscInt q = 0; q < bs2; q++) block[q] = PetscConj(block[q]);
687: }
688: }
689: if (same == PETSC_FALSE && implicit == PETSC_TRUE && i != j) {
690: r = irow[j];
691: colidx = icol[i];
692: pos = bi[r] + bilen[r]++;
693: bj[pos] = colidx;
694: block = ba + pos * bs2;
695: PetscCall(PetscArraycpy(block, a->a + k * bs2, bs2));
696: PetscCall(PetscKernel_A_gets_transpose_A_N(block, bs));
697: if (hermitian == PETSC_TRUE) {
698: for (PetscInt q = 0; q < bs2; q++) block[q] = PetscConj(block[q]);
699: }
700: }
701: }
702: }
703: PetscCall(PetscMalloc1(bs2, &work));
704: for (PetscInt i = 0; i < mbs; i++) PetscCall(PetscSortIntWithDataArray(bilen[i], PetscSafePointerPlusOffset(bj, bi[i]), PetscSafePointerPlusOffset(ba, bi[i] * bs2), bs2 * sizeof(MatScalar), work));
705: PetscCall(PetscFree(work));
706: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
707: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
708: if (same == PETSC_TRUE) PetscCall(MatPropagateSymmetryOptions(A, C));
709: PetscCall(PetscFree2(irow, icol));
710: PetscCall(ISRestoreIndices(browp, &row));
711: PetscCall(ISRestoreIndices(bcolp, &col));
712: PetscCall(ISDestroy(&browp));
713: PetscCall(ISDestroy(&bcolp));
714: *B = C;
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
718: PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
719: {
720: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
721: PetscInt *rp, k, low, high, t, ii, jj, row, nrow, i, col, l, rmax, N, lastcol = -1;
722: PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen;
723: PetscInt *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs, stepval;
724: PetscBool roworiented = a->roworiented;
725: const PetscScalar *value = v;
726: MatScalar *ap = NULL, *aa = a->a, *bap;
728: PetscFunctionBegin;
729: if (roworiented) stepval = (n - 1) * bs;
730: else stepval = (m - 1) * bs;
731: for (k = 0; k < m; k++) { /* loop over added rows */
732: row = im[k];
733: if (row < 0) continue;
734: PetscCheck(row < a->mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block index row too large %" PetscInt_FMT " max %" PetscInt_FMT, row, a->mbs - 1);
735: rp = aj + ai[row];
736: ap = PetscSafePointerPlusOffset(aa, bs2 * ai[row]);
737: rmax = imax[row];
738: nrow = ailen[row];
739: low = 0;
740: high = nrow;
741: for (l = 0; l < n; l++) { /* loop over added columns */
742: if (in[l] < 0) continue;
743: col = in[l];
744: PetscCheck(col < a->nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block index column too large %" PetscInt_FMT " max %" PetscInt_FMT, col, a->nbs - 1);
745: if (col < row) {
746: 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)");
747: continue; /* ignore lower triangular block */
748: }
749: if (!A->structure_only) {
750: if (roworiented) value = v + k * (stepval + bs) * bs + l * bs;
751: else value = v + l * (stepval + bs) * bs + k * bs;
752: }
754: if (col <= lastcol) low = 0;
755: else high = nrow;
757: lastcol = col;
758: while (high - low > 7) {
759: t = (low + high) / 2;
760: if (rp[t] > col) high = t;
761: else low = t;
762: }
763: for (i = low; i < high; i++) {
764: if (rp[i] > col) break;
765: if (rp[i] == col) {
766: if (A->structure_only) goto noinsert2;
767: bap = ap + bs2 * i;
768: if (roworiented) {
769: if (is == ADD_VALUES) {
770: for (ii = 0; ii < bs; ii++, value += stepval) {
771: for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
772: }
773: } else {
774: for (ii = 0; ii < bs; ii++, value += stepval) {
775: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
776: }
777: }
778: } else {
779: if (is == ADD_VALUES) {
780: for (ii = 0; ii < bs; ii++, value += stepval) {
781: for (jj = 0; jj < bs; jj++) *bap++ += *value++;
782: }
783: } else {
784: for (ii = 0; ii < bs; ii++, value += stepval) {
785: for (jj = 0; jj < bs; jj++) *bap++ = *value++;
786: }
787: }
788: }
789: goto noinsert2;
790: }
791: }
792: if (nonew == 1) goto noinsert2;
793: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new block index nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col);
794: if (A->structure_only) MatSeqXAIJReallocateAIJ_structure_only(A, a->mbs, bs2, nrow, row, col, rmax, ai, aj, rp, imax, nonew, MatScalar);
795: else MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
796: N = nrow++ - 1;
797: high++;
798: /* shift up all the later entries in this row */
799: PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
800: rp[i] = col;
801: if (!A->structure_only) {
802: PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
803: PetscCall(PetscArrayzero(ap + bs2 * i, bs2));
804: bap = ap + bs2 * i;
805: if (roworiented) {
806: for (ii = 0; ii < bs; ii++, value += stepval) {
807: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
808: }
809: } else {
810: for (ii = 0; ii < bs; ii++, value += stepval) {
811: for (jj = 0; jj < bs; jj++) *bap++ = *value++;
812: }
813: }
814: }
815: noinsert2:;
816: low = i;
817: }
818: ailen[row] = nrow;
819: }
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }
823: static PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A, MatAssemblyType mode)
824: {
825: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
826: PetscInt fshift = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax;
827: PetscInt m = A->rmap->N, *ip, N, *ailen = a->ilen;
828: PetscInt mbs = a->mbs, bs2 = a->bs2, rmax = 0;
829: MatScalar *aa = a->a, *ap;
831: PetscFunctionBegin;
832: if (mode == MAT_FLUSH_ASSEMBLY || (A->was_assembled && A->ass_nonzerostate == A->nonzerostate)) PetscFunctionReturn(PETSC_SUCCESS);
834: if (m) rmax = ailen[0];
835: for (i = 1; i < mbs; i++) {
836: /* move each row back by the amount of empty slots (fshift) before it*/
837: fshift += imax[i - 1] - ailen[i - 1];
838: rmax = PetscMax(rmax, ailen[i]);
839: if (fshift) {
840: ip = aj + ai[i];
841: N = ailen[i];
842: PetscCall(PetscArraymove(ip - fshift, ip, N));
843: if (!A->structure_only) {
844: ap = aa + bs2 * ai[i];
845: PetscCall(PetscArraymove(ap - bs2 * fshift, ap, bs2 * N));
846: }
847: }
848: ai[i] = ai[i - 1] + ailen[i - 1];
849: }
850: if (mbs) {
851: fshift += imax[mbs - 1] - ailen[mbs - 1];
852: ai[mbs] = ai[mbs - 1] + ailen[mbs - 1];
853: }
854: /* reset ilen and imax for each row */
855: for (i = 0; i < mbs; i++) ailen[i] = imax[i] = ai[i + 1] - ai[i];
856: a->nz = ai[mbs];
858: PetscCheck(!fshift || a->nounused != -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT " block size %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, A->rmap->bs, fshift * bs2);
860: PetscCall(PetscInfo(A, "Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT ", block size %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded, %" PetscInt_FMT " used\n", m, A->rmap->N, A->rmap->bs, fshift * bs2, a->nz * bs2));
861: PetscCall(PetscInfo(A, "Number of mallocs during MatSetValues is %" PetscInt_FMT "\n", a->reallocs));
862: PetscCall(PetscInfo(A, "Most nonzeros blocks in any row is %" PetscInt_FMT "\n", rmax));
864: A->info.mallocs += a->reallocs;
865: a->reallocs = 0;
866: A->info.nz_unneeded = (PetscReal)fshift * bs2;
867: a->rmax = rmax;
869: if (!A->structure_only && A->cmap->n < 65536 && A->cmap->bs == 1) {
870: if (a->jshort && a->free_jshort) {
871: /* when matrix data structure is changed, previous jshort must be replaced */
872: PetscCall(PetscFree(a->jshort));
873: }
874: PetscCall(PetscMalloc1(a->i[A->rmap->n], &a->jshort));
875: for (i = 0; i < a->i[A->rmap->n]; i++) a->jshort[i] = (short)a->j[i];
876: A->ops->mult = MatMult_SeqSBAIJ_1_ushort;
877: A->ops->sor = MatSOR_SeqSBAIJ_ushort;
878: a->free_jshort = PETSC_TRUE;
879: }
880: PetscFunctionReturn(PETSC_SUCCESS);
881: }
883: /* Only add/insert a(i,j) with i<=j (blocks).
884: Any a(i,j) with i>j input by user is ignored.
885: */
887: PetscErrorCode MatSetValues_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
888: {
889: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
890: PetscInt *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N, lastcol = -1;
891: PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen, roworiented = a->roworiented;
892: PetscInt *aj = a->j, nonew = a->nonew, bs = A->rmap->bs, brow, bcol;
893: PetscInt ridx, cidx, bs2 = a->bs2;
894: MatScalar *ap, value = 0.0, *aa = a->a, *bap;
896: PetscFunctionBegin;
897: for (k = 0; k < m; k++) { /* loop over added rows */
898: row = im[k]; /* row number */
899: brow = row / bs; /* block row number */
900: if (row < 0) continue;
901: PetscCheck(row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, A->rmap->N - 1);
902: rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
903: ap = PetscSafePointerPlusOffset(aa, bs2 * ai[brow]); /*ptr to beginning of element value of the row block*/
904: rmax = imax[brow]; /* maximum space allocated for this row */
905: nrow = ailen[brow]; /* actual length of this row */
906: low = 0;
907: high = nrow;
908: for (l = 0; l < n; l++) { /* loop over added columns */
909: if (in[l] < 0) continue;
910: PetscCheck(in[l] < A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[l], A->cmap->N - 1);
911: col = in[l];
912: bcol = col / bs; /* block col number */
914: if (brow > bcol) {
915: 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)");
916: continue; /* ignore lower triangular values */
917: }
919: ridx = row % bs;
920: cidx = col % bs; /*row and col index inside the block */
921: if ((brow == bcol && ridx <= cidx) || (brow < bcol)) {
922: /* element value a(k,l) */
923: if (!A->structure_only) {
924: if (roworiented) value = v[l + k * n];
925: else value = v[k + l * m];
926: }
928: /* move pointer bap to a(k,l) quickly and add/insert value */
929: if (col <= lastcol) low = 0;
930: else high = nrow;
932: lastcol = col;
933: while (high - low > 7) {
934: t = (low + high) / 2;
935: if (rp[t] > bcol) high = t;
936: else low = t;
937: }
938: for (i = low; i < high; i++) {
939: if (rp[i] > bcol) break;
940: if (rp[i] == bcol) {
941: if (A->structure_only) goto noinsert1;
942: bap = ap + bs2 * i + bs * cidx + ridx;
943: if (is == ADD_VALUES) *bap += value;
944: else *bap = value;
945: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
946: if (brow == bcol && ridx < cidx) {
947: bap = ap + bs2 * i + bs * ridx + cidx;
948: if (is == ADD_VALUES) *bap += value;
949: else *bap = value;
950: }
951: goto noinsert1;
952: }
953: }
955: if (nonew == 1) goto noinsert1;
956: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col);
957: if (A->structure_only) MatSeqXAIJReallocateAIJ_structure_only(A, a->mbs, bs2, nrow, brow, bcol, rmax, ai, aj, rp, imax, nonew, MatScalar);
958: else MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
960: N = nrow++ - 1;
961: high++;
962: /* shift up all the later entries in this row */
963: PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
964: rp[i] = bcol;
965: if (!A->structure_only) {
966: PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
967: PetscCall(PetscArrayzero(ap + bs2 * i, bs2));
968: ap[bs2 * i + bs * cidx + ridx] = value;
969: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
970: if (brow == bcol && ridx < cidx) ap[bs2 * i + bs * ridx + cidx] = value;
971: }
972: noinsert1:;
973: low = i;
974: }
975: } /* end of loop over added columns */
976: ailen[brow] = nrow;
977: } /* end of loop over added rows */
978: PetscFunctionReturn(PETSC_SUCCESS);
979: }
981: static PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA, IS row, const MatFactorInfo *info)
982: {
983: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)inA->data;
984: Mat outA;
985: PetscBool row_identity;
987: PetscFunctionBegin;
988: PetscCheck(info->levels == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only levels=0 is supported for in-place icc");
989: PetscCall(ISIdentity(row, &row_identity));
990: PetscCheck(row_identity, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported");
991: PetscCheck(inA->rmap->bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix block size %" PetscInt_FMT " is not supported", inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */
993: outA = inA;
994: PetscCall(PetscFree(inA->solvertype));
995: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &inA->solvertype));
997: inA->factortype = MAT_FACTOR_ICC;
998: PetscCall(MatSeqSBAIJSetNumericFactorization_inplace(inA, row_identity));
1000: PetscCall(PetscObjectReference((PetscObject)row));
1001: PetscCall(ISDestroy(&a->row));
1002: a->row = row;
1003: PetscCall(PetscObjectReference((PetscObject)row));
1004: PetscCall(ISDestroy(&a->col));
1005: a->col = row;
1007: /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
1008: if (a->icol) PetscCall(ISInvertPermutation(row, PETSC_DECIDE, &a->icol));
1010: if (!a->solve_work) PetscCall(PetscMalloc1(inA->rmap->N + inA->rmap->bs, &a->solve_work));
1012: PetscCall(MatCholeskyFactorNumeric(outA, inA, info));
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: static PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat, PetscInt *indices)
1017: {
1018: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data;
1019: PetscInt i, nz, n;
1021: PetscFunctionBegin;
1022: nz = baij->maxnz;
1023: n = mat->cmap->n;
1024: for (i = 0; i < nz; i++) baij->j[i] = indices[i];
1026: baij->nz = nz;
1027: for (i = 0; i < n; i++) baij->ilen[i] = baij->imax[i];
1029: PetscCall(MatSetOption(mat, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: /*@
1034: MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1035: in a `MATSEQSBAIJ` matrix.
1037: Input Parameters:
1038: + mat - the `MATSEQSBAIJ` matrix
1039: - indices - the column indices
1041: Level: advanced
1043: Notes:
1044: This can be called if you have precomputed the nonzero structure of the
1045: matrix and want to provide it to the matrix object to improve the performance
1046: of the `MatSetValues()` operation.
1048: You MUST have set the correct numbers of nonzeros per row in the call to
1049: `MatCreateSeqSBAIJ()`, and the columns indices MUST be sorted.
1051: MUST be called before any calls to `MatSetValues()`
1053: .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatCreateSeqSBAIJ`
1054: @*/
1055: PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat, PetscInt *indices)
1056: {
1057: PetscFunctionBegin;
1059: PetscAssertPointer(indices, 2);
1060: PetscUseMethod(mat, "MatSeqSBAIJSetColumnIndices_C", (Mat, PetscInt *), (mat, indices));
1061: PetscFunctionReturn(PETSC_SUCCESS);
1062: }
1064: static PetscErrorCode MatCopy_SeqSBAIJ(Mat A, Mat B, MatStructure str)
1065: {
1066: PetscBool isbaij;
1068: PetscFunctionBegin;
1069: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &isbaij, MATSEQSBAIJ, MATMPISBAIJ, ""));
1070: PetscCheck(isbaij, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)B)->type_name);
1071: /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */
1072: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1073: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1074: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data;
1076: PetscCheck(a->i[a->mbs] == b->i[b->mbs], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of nonzeros in two matrices are different");
1077: PetscCheck(a->mbs == b->mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of rows in two matrices are different");
1078: PetscCheck(a->bs2 == b->bs2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different block size");
1079: PetscCall(PetscArraycpy(b->a, a->a, a->bs2 * a->i[a->mbs]));
1080: PetscCall(PetscObjectStateIncrease((PetscObject)B));
1081: } else {
1082: PetscCall(MatGetRowUpperTriangular(A));
1083: PetscCall(MatCopy_Basic(A, B, str));
1084: PetscCall(MatRestoreRowUpperTriangular(A));
1085: }
1086: PetscFunctionReturn(PETSC_SUCCESS);
1087: }
1089: static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A, PetscScalar *array[])
1090: {
1091: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1093: PetscFunctionBegin;
1094: *array = a->a;
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A, PetscScalar *array[])
1099: {
1100: PetscFunctionBegin;
1101: *array = NULL;
1102: PetscFunctionReturn(PETSC_SUCCESS);
1103: }
1105: PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y, Mat X, PetscInt *nnz)
1106: {
1107: PetscInt bs = Y->rmap->bs, mbs = Y->rmap->N / bs;
1108: Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ *)X->data;
1109: Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ *)Y->data;
1111: PetscFunctionBegin;
1112: /* Set the number of nonzeros in the new matrix */
1113: PetscCall(MatAXPYGetPreallocation_SeqX_private(mbs, x->i, x->j, y->i, y->j, nnz));
1114: PetscFunctionReturn(PETSC_SUCCESS);
1115: }
1117: static PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str)
1118: {
1119: Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ *)X->data, *y = (Mat_SeqSBAIJ *)Y->data;
1120: PetscInt bs = Y->rmap->bs, bs2 = bs * bs;
1121: PetscBLASInt one = 1;
1123: PetscFunctionBegin;
1124: if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) {
1125: PetscBool e = x->nz == y->nz && x->mbs == y->mbs ? PETSC_TRUE : PETSC_FALSE;
1126: if (e) {
1127: PetscCall(PetscArraycmp(x->i, y->i, x->mbs + 1, &e));
1128: if (e) {
1129: PetscCall(PetscArraycmp(x->j, y->j, x->i[x->mbs], &e));
1130: if (e) str = SAME_NONZERO_PATTERN;
1131: }
1132: }
1133: if (!e) PetscCheck(str != SAME_NONZERO_PATTERN, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatStructure is not SAME_NONZERO_PATTERN");
1134: }
1135: if (str == SAME_NONZERO_PATTERN) {
1136: PetscScalar alpha = a;
1137: PetscBLASInt bnz;
1138: PetscCall(PetscBLASIntCast(x->nz * bs2, &bnz));
1139: PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one));
1140: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
1141: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1142: PetscCall(MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE));
1143: PetscCall(MatAXPY_Basic(Y, a, X, str));
1144: PetscCall(MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE));
1145: } else {
1146: Mat B;
1147: PetscInt *nnz;
1148: PetscCheck(bs == X->rmap->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrices must have same block size");
1149: PetscCall(MatGetRowUpperTriangular(X));
1150: PetscCall(MatGetRowUpperTriangular(Y));
1151: PetscCall(PetscMalloc1(Y->rmap->N, &nnz));
1152: PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B));
1153: PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name));
1154: PetscCall(MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N));
1155: PetscCall(MatSetBlockSizesFromMats(B, Y, Y));
1156: PetscCall(MatSetType(B, ((PetscObject)Y)->type_name));
1157: PetscCall(MatAXPYGetPreallocation_SeqSBAIJ(Y, X, nnz));
1158: PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, nnz));
1160: PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
1162: PetscCall(MatHeaderMerge(Y, &B));
1163: PetscCall(PetscFree(nnz));
1164: PetscCall(MatRestoreRowUpperTriangular(X));
1165: PetscCall(MatRestoreRowUpperTriangular(Y));
1166: }
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: static PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A, PetscBool *flg)
1171: {
1172: PetscFunctionBegin;
1173: *flg = PETSC_TRUE;
1174: PetscFunctionReturn(PETSC_SUCCESS);
1175: }
1177: static PetscErrorCode MatConjugate_SeqSBAIJ(Mat A)
1178: {
1179: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1180: PetscInt i, nz = a->bs2 * a->i[a->mbs];
1181: MatScalar *aa = a->a;
1183: PetscFunctionBegin;
1184: for (i = 0; i < nz; i++) aa[i] = PetscConj(aa[i]);
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: static PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1189: {
1190: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1191: PetscInt i, nz = a->bs2 * a->i[a->mbs];
1192: MatScalar *aa = a->a;
1194: PetscFunctionBegin;
1195: for (i = 0; i < nz; i++) aa[i] = PetscRealPart(aa[i]);
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: static PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1200: {
1201: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1202: PetscInt i, nz = a->bs2 * a->i[a->mbs];
1203: MatScalar *aa = a->a;
1205: PetscFunctionBegin;
1206: for (i = 0; i < nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: static PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A, PetscInt is_n, const PetscInt is_idx[], PetscScalar diag, Vec x, Vec b)
1211: {
1212: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)A->data;
1213: PetscInt i, j, k, count;
1214: PetscInt bs = A->rmap->bs, bs2 = baij->bs2, row, col;
1215: PetscScalar zero = 0.0;
1216: MatScalar *aa;
1217: const PetscScalar *xx;
1218: PetscScalar *bb;
1219: PetscBool *zeroed, vecs = PETSC_FALSE;
1221: PetscFunctionBegin;
1222: /* fix right-hand side if needed */
1223: if (x && b) {
1224: PetscCall(VecGetArrayRead(x, &xx));
1225: PetscCall(VecGetArray(b, &bb));
1226: vecs = PETSC_TRUE;
1227: }
1229: /* zero the columns */
1230: PetscCall(PetscCalloc1(A->rmap->n, &zeroed));
1231: for (i = 0; i < is_n; i++) {
1232: PetscCheck(is_idx[i] >= 0 && is_idx[i] < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", is_idx[i]);
1233: zeroed[is_idx[i]] = PETSC_TRUE;
1234: }
1235: if (vecs) {
1236: for (i = 0; i < A->rmap->N; i++) {
1237: row = i / bs;
1238: for (j = baij->i[row]; j < baij->i[row + 1]; j++) {
1239: for (k = 0; k < bs; k++) {
1240: col = bs * baij->j[j] + k;
1241: if (col <= i) continue;
1242: aa = baij->a + j * bs2 + (i % bs) + bs * k;
1243: if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0] * xx[col];
1244: if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0] * xx[i];
1245: }
1246: }
1247: }
1248: for (i = 0; i < is_n; i++) bb[is_idx[i]] = diag * xx[is_idx[i]];
1249: }
1251: for (i = 0; i < A->rmap->N; i++) {
1252: if (!zeroed[i]) {
1253: row = i / bs;
1254: for (j = baij->i[row]; j < baij->i[row + 1]; j++) {
1255: for (k = 0; k < bs; k++) {
1256: col = bs * baij->j[j] + k;
1257: if (zeroed[col]) {
1258: aa = baij->a + j * bs2 + (i % bs) + bs * k;
1259: aa[0] = 0.0;
1260: }
1261: }
1262: }
1263: }
1264: }
1265: PetscCall(PetscFree(zeroed));
1266: if (vecs) {
1267: PetscCall(VecRestoreArrayRead(x, &xx));
1268: PetscCall(VecRestoreArray(b, &bb));
1269: }
1271: /* zero the rows */
1272: for (i = 0; i < is_n; i++) {
1273: row = is_idx[i];
1274: count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
1275: aa = PetscSafePointerPlusOffset(baij->a, baij->i[row / bs] * bs2 + (row % bs));
1276: for (k = 0; k < count; k++) {
1277: aa[0] = zero;
1278: aa += bs;
1279: }
1280: if (diag != 0.0) PetscUseTypeMethod(A, setvalues, 1, &row, 1, &row, &diag, INSERT_VALUES);
1281: }
1282: PetscCall(MatAssemblyEnd_SeqSBAIJ(A, MAT_FINAL_ASSEMBLY));
1283: PetscFunctionReturn(PETSC_SUCCESS);
1284: }
1286: static PetscErrorCode MatShift_SeqSBAIJ(Mat Y, PetscScalar a)
1287: {
1288: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)Y->data;
1290: PetscFunctionBegin;
1291: if (!Y->preallocated || !aij->nz) PetscCall(MatSeqSBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL));
1292: PetscCall(MatShift_Basic(Y, a));
1293: PetscFunctionReturn(PETSC_SUCCESS);
1294: }
1296: PetscErrorCode MatEliminateZeros_SeqSBAIJ(Mat A, PetscBool keep)
1297: {
1298: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1299: PetscInt fshift = 0, fshift_prev = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax, j, k;
1300: PetscInt m = A->rmap->N, *ailen = a->ilen;
1301: PetscInt mbs = a->mbs, bs2 = a->bs2, rmax = 0;
1302: MatScalar *aa = a->a, *ap;
1303: PetscBool zero;
1305: PetscFunctionBegin;
1306: PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot eliminate zeros for unassembled matrix");
1307: if (m) rmax = ailen[0];
1308: for (i = 1, a->nonzerorowcnt = 0; i <= mbs; i++) {
1309: for (k = ai[i - 1]; k < ai[i]; k++) {
1310: zero = PETSC_TRUE;
1311: ap = aa + bs2 * k;
1312: for (j = 0; j < bs2 && zero; j++) {
1313: if (ap[j] != 0.0) zero = PETSC_FALSE;
1314: }
1315: if (zero && (aj[k] != i - 1 || !keep)) fshift++;
1316: else {
1317: if (zero && aj[k] == i - 1) PetscCall(PetscInfo(A, "Keep the diagonal block at row %" PetscInt_FMT "\n", i - 1));
1318: aj[k - fshift] = aj[k];
1319: PetscCall(PetscArraymove(ap - bs2 * fshift, ap, bs2));
1320: }
1321: }
1322: ai[i - 1] -= fshift_prev;
1323: fshift_prev = fshift;
1324: ailen[i - 1] = imax[i - 1] = ai[i] - fshift - ai[i - 1];
1325: a->nonzerorowcnt += ((ai[i] - fshift - ai[i - 1]) > 0);
1326: rmax = PetscMax(rmax, ailen[i - 1]);
1327: }
1328: if (fshift) {
1329: if (mbs) {
1330: ai[mbs] -= fshift;
1331: a->nz = ai[mbs];
1332: }
1333: PetscCall(PetscInfo(A, "Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; zeros eliminated: %" PetscInt_FMT "; nonzeros left: %" PetscInt_FMT "\n", m, A->cmap->n, fshift, a->nz));
1334: A->nonzerostate++;
1335: A->info.nz_unneeded += (PetscReal)fshift;
1336: a->rmax = rmax;
1337: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1338: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1339: }
1340: PetscFunctionReturn(PETSC_SUCCESS);
1341: }
1343: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1344: MatGetRow_SeqSBAIJ,
1345: MatRestoreRow_SeqSBAIJ,
1346: MatMult_SeqSBAIJ_N,
1347: /* 4*/ MatMultAdd_SeqSBAIJ_N,
1348: MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1349: MatMultAdd_SeqSBAIJ_N,
1350: NULL,
1351: NULL,
1352: NULL,
1353: /* 10*/ NULL,
1354: NULL,
1355: MatCholeskyFactor_SeqSBAIJ,
1356: MatSOR_SeqSBAIJ,
1357: MatTranspose_SeqSBAIJ,
1358: /* 15*/ MatGetInfo_SeqSBAIJ,
1359: MatEqual_SeqSBAIJ,
1360: MatGetDiagonal_SeqSBAIJ,
1361: MatDiagonalScale_SeqSBAIJ,
1362: MatNorm_SeqSBAIJ,
1363: /* 20*/ NULL,
1364: MatAssemblyEnd_SeqSBAIJ,
1365: MatSetOption_SeqSBAIJ,
1366: MatZeroEntries_SeqSBAIJ,
1367: /* 24*/ NULL,
1368: NULL,
1369: NULL,
1370: NULL,
1371: NULL,
1372: /* 29*/ MatSetUp_Seq_Hash,
1373: NULL,
1374: NULL,
1375: NULL,
1376: NULL,
1377: /* 34*/ MatDuplicate_SeqSBAIJ,
1378: NULL,
1379: NULL,
1380: NULL,
1381: MatICCFactor_SeqSBAIJ,
1382: /* 39*/ MatAXPY_SeqSBAIJ,
1383: MatCreateSubMatrices_SeqSBAIJ,
1384: MatIncreaseOverlap_SeqSBAIJ,
1385: MatGetValues_SeqSBAIJ,
1386: MatCopy_SeqSBAIJ,
1387: /* 44*/ NULL,
1388: MatScale_SeqSBAIJ,
1389: MatShift_SeqSBAIJ,
1390: NULL,
1391: MatZeroRowsColumns_SeqSBAIJ,
1392: /* 49*/ NULL,
1393: MatGetRowIJ_SeqSBAIJ,
1394: MatRestoreRowIJ_SeqSBAIJ,
1395: NULL,
1396: NULL,
1397: /* 54*/ NULL,
1398: NULL,
1399: NULL,
1400: MatPermute_SeqSBAIJ,
1401: MatSetValuesBlocked_SeqSBAIJ,
1402: /* 59*/ MatCreateSubMatrix_SeqSBAIJ,
1403: NULL,
1404: NULL,
1405: NULL,
1406: NULL,
1407: /* 64*/ NULL,
1408: NULL,
1409: NULL,
1410: NULL,
1411: MatGetRowMaxAbs_SeqSBAIJ,
1412: /* 69*/ NULL,
1413: MatConvert_MPISBAIJ_Basic,
1414: NULL,
1415: NULL,
1416: NULL,
1417: /* 74*/ NULL,
1418: NULL,
1419: NULL,
1420: MatGetInertia_SeqSBAIJ,
1421: MatLoad_SeqSBAIJ,
1422: /* 79*/ NULL,
1423: NULL,
1424: MatIsStructurallySymmetric_SeqSBAIJ,
1425: NULL,
1426: NULL,
1427: /* 84*/ NULL,
1428: NULL,
1429: NULL,
1430: NULL,
1431: NULL,
1432: /* 89*/ NULL,
1433: NULL,
1434: NULL,
1435: NULL,
1436: MatConjugate_SeqSBAIJ,
1437: /* 94*/ NULL,
1438: NULL,
1439: MatRealPart_SeqSBAIJ,
1440: MatImaginaryPart_SeqSBAIJ,
1441: MatGetRowUpperTriangular_SeqSBAIJ,
1442: /* 99*/ MatRestoreRowUpperTriangular_SeqSBAIJ,
1443: NULL,
1444: NULL,
1445: NULL,
1446: NULL,
1447: /*104*/ NULL,
1448: NULL,
1449: NULL,
1450: NULL,
1451: NULL,
1452: /*109*/ NULL,
1453: NULL,
1454: NULL,
1455: NULL,
1456: NULL,
1457: /*114*/ NULL,
1458: MatGetColumnReductions_SeqSBAIJ,
1459: NULL,
1460: NULL,
1461: NULL,
1462: /*119*/ NULL,
1463: NULL,
1464: NULL,
1465: NULL,
1466: NULL,
1467: /*124*/ NULL,
1468: MatSetBlockSizes_Default,
1469: NULL,
1470: NULL,
1471: NULL,
1472: /*129*/ MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ,
1473: NULL,
1474: NULL,
1475: NULL,
1476: NULL,
1477: /*134*/ NULL,
1478: MatEliminateZeros_SeqSBAIJ,
1479: NULL,
1480: NULL,
1481: NULL,
1482: /*139*/ NULL,
1483: MatCopyHashToXAIJ_Seq_Hash,
1484: NULL,
1485: NULL,
1486: NULL,
1487: /*144*/ NULL,
1488: NULL,
1489: NULL,
1490: NULL};
1492: static PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1493: {
1494: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data;
1495: PetscInt nz = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2;
1497: PetscFunctionBegin;
1498: PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1500: /* allocate space for values if not already there */
1501: if (!aij->saved_values) PetscCall(PetscMalloc1(nz + 1, &aij->saved_values));
1503: /* copy values over */
1504: PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz));
1505: PetscFunctionReturn(PETSC_SUCCESS);
1506: }
1508: static PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1509: {
1510: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data;
1511: PetscInt nz = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2;
1513: PetscFunctionBegin;
1514: PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1515: PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first");
1517: /* copy values over */
1518: PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz));
1519: PetscFunctionReturn(PETSC_SUCCESS);
1520: }
1522: static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[])
1523: {
1524: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data;
1525: PetscInt i, mbs, nbs, bs2;
1526: PetscBool skipallocation = PETSC_FALSE, flg = PETSC_FALSE, realalloc = PETSC_FALSE;
1528: PetscFunctionBegin;
1529: if (B->hash_active) {
1530: PetscInt bs;
1531: B->ops[0] = b->cops;
1532: PetscCall(PetscHMapIJVDestroy(&b->ht));
1533: PetscCall(MatGetBlockSize(B, &bs));
1534: if (bs > 1) PetscCall(PetscHSetIJDestroy(&b->bht));
1535: PetscCall(PetscFree(b->dnz));
1536: PetscCall(PetscFree(b->bdnz));
1537: B->hash_active = PETSC_FALSE;
1538: }
1539: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
1541: PetscCall(MatSetBlockSize(B, bs));
1542: PetscCall(PetscLayoutSetUp(B->rmap));
1543: PetscCall(PetscLayoutSetUp(B->cmap));
1544: PetscCheck(B->rmap->N <= B->cmap->N, PETSC_COMM_SELF, PETSC_ERR_SUP, "SEQSBAIJ matrix cannot have more rows %" PetscInt_FMT " than columns %" PetscInt_FMT, B->rmap->N, B->cmap->N);
1545: PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
1547: B->preallocated = PETSC_TRUE;
1549: mbs = B->rmap->N / bs;
1550: nbs = B->cmap->n / bs;
1551: bs2 = bs * bs;
1553: PetscCheck(mbs * bs == B->rmap->N && nbs * bs == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows, cols must be divisible by blocksize");
1555: if (nz == MAT_SKIP_ALLOCATION) {
1556: skipallocation = PETSC_TRUE;
1557: nz = 0;
1558: }
1560: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1561: PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz);
1562: if (nnz) {
1563: for (i = 0; i < mbs; i++) {
1564: PetscCheck(nnz[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, nnz[i]);
1565: PetscCheck(nnz[i] <= nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be greater than block row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " block rowlength %" PetscInt_FMT, i, nnz[i], nbs);
1566: }
1567: }
1569: B->ops->mult = MatMult_SeqSBAIJ_N;
1570: B->ops->multadd = MatMultAdd_SeqSBAIJ_N;
1571: B->ops->multtranspose = MatMult_SeqSBAIJ_N;
1572: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N;
1574: PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL));
1575: if (!flg) {
1576: switch (bs) {
1577: case 1:
1578: B->ops->mult = MatMult_SeqSBAIJ_1;
1579: B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1580: B->ops->multtranspose = MatMult_SeqSBAIJ_1;
1581: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1582: break;
1583: case 2:
1584: B->ops->mult = MatMult_SeqSBAIJ_2;
1585: B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1586: B->ops->multtranspose = MatMult_SeqSBAIJ_2;
1587: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1588: break;
1589: case 3:
1590: B->ops->mult = MatMult_SeqSBAIJ_3;
1591: B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1592: B->ops->multtranspose = MatMult_SeqSBAIJ_3;
1593: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1594: break;
1595: case 4:
1596: B->ops->mult = MatMult_SeqSBAIJ_4;
1597: B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1598: B->ops->multtranspose = MatMult_SeqSBAIJ_4;
1599: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1600: break;
1601: case 5:
1602: B->ops->mult = MatMult_SeqSBAIJ_5;
1603: B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1604: B->ops->multtranspose = MatMult_SeqSBAIJ_5;
1605: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1606: break;
1607: case 6:
1608: B->ops->mult = MatMult_SeqSBAIJ_6;
1609: B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1610: B->ops->multtranspose = MatMult_SeqSBAIJ_6;
1611: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1612: break;
1613: case 7:
1614: B->ops->mult = MatMult_SeqSBAIJ_7;
1615: B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1616: B->ops->multtranspose = MatMult_SeqSBAIJ_7;
1617: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1618: break;
1619: }
1620: }
1622: b->mbs = mbs;
1623: b->nbs = nbs;
1624: if (!skipallocation) {
1625: if (!b->imax) {
1626: PetscCall(PetscMalloc2(mbs, &b->imax, mbs, &b->ilen));
1627: b->free_imax_ilen = PETSC_TRUE;
1628: }
1629: if (!nnz) {
1630: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1631: else if (nz <= 0) nz = 1;
1632: nz = PetscMin(nbs, nz);
1633: for (i = 0; i < mbs; i++) b->imax[i] = nz;
1634: PetscCall(PetscIntMultError(nz, mbs, &nz));
1635: } else {
1636: PetscInt64 nz64 = 0;
1637: for (i = 0; i < mbs; i++) {
1638: b->imax[i] = nnz[i];
1639: nz64 += nnz[i];
1640: }
1641: PetscCall(PetscIntCast(nz64, &nz));
1642: }
1643: /* b->ilen will count nonzeros in each block row so far. */
1644: for (i = 0; i < mbs; i++) b->ilen[i] = 0;
1645: /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1647: /* allocate the matrix space */
1648: PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i));
1649: if (!B->structure_only) {
1650: PetscCall(PetscShmgetAllocateArray(bs2 * nz, sizeof(PetscScalar), (void **)&b->a));
1651: PetscCall(PetscArrayzero(b->a, nz * bs2));
1652: }
1653: PetscCall(PetscShmgetAllocateArray(nz, sizeof(PetscInt), (void **)&b->j));
1654: PetscCall(PetscShmgetAllocateArray(B->rmap->n + 1, sizeof(PetscInt), (void **)&b->i));
1655: PetscCall(PetscArrayzero(b->j, nz));
1656: b->free_a = PETSC_TRUE;
1657: b->free_ij = PETSC_TRUE;
1659: /* pointer to beginning of each row */
1660: b->i[0] = 0;
1661: for (i = 1; i < mbs + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1];
1662: } else {
1663: b->free_a = PETSC_FALSE;
1664: b->free_ij = PETSC_FALSE;
1665: }
1667: b->bs2 = bs2;
1668: b->nz = 0;
1669: b->maxnz = nz;
1670: b->inew = NULL;
1671: b->jnew = NULL;
1672: b->anew = NULL;
1673: b->a2anew = NULL;
1674: b->permute = PETSC_FALSE;
1676: B->was_assembled = PETSC_FALSE;
1677: B->assembled = PETSC_FALSE;
1678: if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: static PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[])
1683: {
1684: PetscInt i, j, m, nz, anz, nz_max = 0, *nnz;
1685: PetscScalar *values = NULL;
1686: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data;
1687: PetscBool roworiented = b->roworiented;
1688: PetscBool ilw = b->ignore_ltriangular;
1690: PetscFunctionBegin;
1691: PetscCheck(bs >= 1, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "Invalid block size specified, must be positive but it is %" PetscInt_FMT, bs);
1692: PetscCall(PetscLayoutSetBlockSize(B->rmap, bs));
1693: PetscCall(PetscLayoutSetBlockSize(B->cmap, bs));
1694: PetscCall(PetscLayoutSetUp(B->rmap));
1695: PetscCall(PetscLayoutSetUp(B->cmap));
1696: PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
1697: m = B->rmap->n / bs;
1699: PetscCheck(!ii[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]);
1700: PetscCall(PetscMalloc1(m + 1, &nnz));
1701: for (i = 0; i < m; i++) {
1702: nz = ii[i + 1] - ii[i];
1703: PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz);
1704: PetscCheckSorted(nz, jj + ii[i]);
1705: anz = 0;
1706: for (j = 0; j < nz; j++) {
1707: /* count only values on the diagonal or above */
1708: if (jj[ii[i] + j] >= i) {
1709: anz = nz - j;
1710: break;
1711: }
1712: }
1713: nz_max = PetscMax(nz_max, nz);
1714: nnz[i] = anz;
1715: }
1716: PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, nnz));
1717: PetscCall(PetscFree(nnz));
1719: values = (PetscScalar *)V;
1720: if (!values && !B->structure_only) PetscCall(PetscCalloc1(bs * bs * nz_max, &values));
1721: b->ignore_ltriangular = PETSC_TRUE;
1722: for (i = 0; i < m; i++) {
1723: PetscInt ncols = ii[i + 1] - ii[i];
1724: const PetscInt *icols = jj + ii[i];
1726: if (B->structure_only) PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, ncols, icols, NULL, INSERT_VALUES));
1727: else if (!roworiented || bs == 1) {
1728: const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0);
1730: PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, ncols, icols, svals, INSERT_VALUES));
1731: } else {
1732: for (j = 0; j < ncols; j++) {
1733: const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0);
1735: PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, 1, &icols[j], svals, INSERT_VALUES));
1736: }
1737: }
1738: }
1739: if (!V) PetscCall(PetscFree(values));
1740: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1741: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1742: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1743: b->ignore_ltriangular = ilw;
1744: PetscFunctionReturn(PETSC_SUCCESS);
1745: }
1747: /*
1748: This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization
1749: */
1750: PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B, PetscBool natural)
1751: {
1752: PetscBool flg = PETSC_FALSE;
1753: PetscInt bs = B->rmap->bs;
1755: PetscFunctionBegin;
1756: PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL));
1757: if (flg) bs = 8;
1759: if (!natural) {
1760: switch (bs) {
1761: case 1:
1762: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace;
1763: break;
1764: case 2:
1765: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1766: break;
1767: case 3:
1768: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1769: break;
1770: case 4:
1771: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1772: break;
1773: case 5:
1774: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1775: break;
1776: case 6:
1777: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1778: break;
1779: case 7:
1780: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1781: break;
1782: default:
1783: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N;
1784: break;
1785: }
1786: } else {
1787: switch (bs) {
1788: case 1:
1789: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace;
1790: break;
1791: case 2:
1792: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1793: break;
1794: case 3:
1795: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1796: break;
1797: case 4:
1798: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1799: break;
1800: case 5:
1801: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1802: break;
1803: case 6:
1804: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1805: break;
1806: case 7:
1807: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1808: break;
1809: default:
1810: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
1811: break;
1812: }
1813: }
1814: PetscFunctionReturn(PETSC_SUCCESS);
1815: }
1817: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *);
1818: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType, MatReuse, Mat *);
1819: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A, MatSolverType *type)
1820: {
1821: PetscFunctionBegin;
1822: *type = MATSOLVERPETSC;
1823: PetscFunctionReturn(PETSC_SUCCESS);
1824: }
1826: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A, MatFactorType ftype, Mat *B)
1827: {
1828: PetscInt n = A->rmap->n;
1830: PetscFunctionBegin;
1831: if (PetscDefined(USE_COMPLEX) && (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) && A->hermitian == PETSC_BOOL3_TRUE && A->symmetric != PETSC_BOOL3_TRUE) {
1832: PetscCall(PetscInfo(A, "Hermitian MAT_FACTOR_CHOLESKY or MAT_FACTOR_ICC are not supported. Use MAT_FACTOR_LU instead.\n"));
1833: *B = NULL;
1834: PetscFunctionReturn(PETSC_SUCCESS);
1835: }
1837: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
1838: PetscCall(MatSetSizes(*B, n, n, n, n));
1839: PetscCheck(ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC, PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
1840: PetscCall(MatSetType(*B, MATSEQSBAIJ));
1841: PetscCall(MatSeqSBAIJSetPreallocation(*B, A->rmap->bs, MAT_SKIP_ALLOCATION, NULL));
1843: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ;
1844: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;
1845: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]));
1846: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC]));
1848: (*B)->factortype = ftype;
1849: (*B)->canuseordering = PETSC_TRUE;
1850: PetscCall(PetscFree((*B)->solvertype));
1851: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
1852: PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_petsc));
1853: PetscFunctionReturn(PETSC_SUCCESS);
1854: }
1856: /*@
1857: MatSeqSBAIJGetArray - gives access to the array where the numerical data for a `MATSEQSBAIJ` matrix is stored
1859: Not Collective
1861: Input Parameter:
1862: . A - a `MATSEQSBAIJ` matrix
1864: Output Parameter:
1865: . array - pointer to the data
1867: Level: intermediate
1869: .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatSeqSBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()`
1870: @*/
1871: PetscErrorCode MatSeqSBAIJGetArray(Mat A, PetscScalar *array[])
1872: {
1873: PetscFunctionBegin;
1874: PetscUseMethod(A, "MatSeqSBAIJGetArray_C", (Mat, PetscScalar **), (A, array));
1875: PetscFunctionReturn(PETSC_SUCCESS);
1876: }
1878: /*@
1879: MatSeqSBAIJRestoreArray - returns access to the array where the numerical data for a `MATSEQSBAIJ` matrix is stored obtained by `MatSeqSBAIJGetArray()`
1881: Not Collective
1883: Input Parameters:
1884: + A - a `MATSEQSBAIJ` matrix
1885: - array - pointer to the data
1887: Level: intermediate
1889: .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatSeqSBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()`
1890: @*/
1891: PetscErrorCode MatSeqSBAIJRestoreArray(Mat A, PetscScalar *array[])
1892: {
1893: PetscFunctionBegin;
1894: PetscUseMethod(A, "MatSeqSBAIJRestoreArray_C", (Mat, PetscScalar **), (A, array));
1895: PetscCall(PetscObjectStateIncrease((PetscObject)A));
1896: PetscFunctionReturn(PETSC_SUCCESS);
1897: }
1899: /*MC
1900: MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1901: based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored.
1903: For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
1904: can call `MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE)`.
1906: Options Database Key:
1907: . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to `MatSetFromOptions()`
1909: Level: beginner
1911: Notes:
1912: Call `MatSetOption(A, MAT_STRUCTURE_ONLY, PETSC_TRUE)` before preallocation or `MatSetUp()` to store only the nonzero pattern.
1913: The assembled matrix has no numerical value array. Row and column indices supplied during insertion are retained, while numerical values are ignored.
1914: Such matrices can be used for structural operations, but not for numerical operations.
1916: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not
1917: stored and they are assumed to be symmetric to the upper triangular). If you call `MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_FALSE)` or use
1918: the options database `-mat_ignore_lower_triangular false` it will generate an error if you try to set a value in the lower triangular portion.
1920: The number of rows in the matrix must be less than or equal to the number of columns
1922: .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatCreateSeqSBAIJ()`, `MatType`, `MATMPISBAIJ`
1923: M*/
1924: PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1925: {
1926: Mat_SeqSBAIJ *b;
1927: PetscMPIInt size;
1928: PetscBool no_unroll = PETSC_FALSE, no_inode = PETSC_FALSE;
1930: PetscFunctionBegin;
1931: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
1932: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");
1934: PetscCall(PetscNew(&b));
1935: B->data = (void *)b;
1936: B->ops[0] = MatOps_Values;
1938: B->ops->destroy = MatDestroy_SeqSBAIJ;
1939: B->ops->view = MatView_SeqSBAIJ;
1940: b->row = NULL;
1941: b->icol = NULL;
1942: b->reallocs = 0;
1943: b->saved_values = NULL;
1944: b->inode.limit = 5;
1945: b->inode.max_limit = 5;
1947: b->roworiented = PETSC_TRUE;
1948: b->nonew = 0;
1949: b->diag = NULL;
1950: b->solve_work = NULL;
1951: b->mult_work = NULL;
1952: B->spptr = NULL;
1953: B->info.nz_unneeded = (PetscReal)b->maxnz * b->bs2;
1954: b->keepnonzeropattern = PETSC_FALSE;
1956: b->inew = NULL;
1957: b->jnew = NULL;
1958: b->anew = NULL;
1959: b->a2anew = NULL;
1960: b->permute = PETSC_FALSE;
1962: b->ignore_ltriangular = PETSC_TRUE;
1964: PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_ignore_lower_triangular", &b->ignore_ltriangular, NULL));
1966: b->getrow_utriangular = PETSC_FALSE;
1968: PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_getrow_uppertriangular", &b->getrow_utriangular, NULL));
1970: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJGetArray_C", MatSeqSBAIJGetArray_SeqSBAIJ));
1971: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJRestoreArray_C", MatSeqSBAIJRestoreArray_SeqSBAIJ));
1972: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqSBAIJ));
1973: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqSBAIJ));
1974: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetColumnIndices_C", MatSeqSBAIJSetColumnIndices_SeqSBAIJ));
1975: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqaij_C", MatConvert_SeqSBAIJ_SeqAIJ));
1976: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqbaij_C", MatConvert_SeqSBAIJ_SeqBAIJ));
1977: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocation_C", MatSeqSBAIJSetPreallocation_SeqSBAIJ));
1978: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocationCSR_C", MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ));
1979: #if PetscDefined(HAVE_ELEMENTAL)
1980: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_elemental_C", MatConvert_SeqSBAIJ_Elemental));
1981: #endif
1982: #if PetscDefined(HAVE_SCALAPACK) && (PetscDefined(USE_REAL_SINGLE) || PetscDefined(USE_REAL_DOUBLE))
1983: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_scalapack_C", MatConvert_SBAIJ_ScaLAPACK));
1984: #endif
1986: B->symmetry_eternal = PETSC_TRUE;
1987: B->structural_symmetry_eternal = PETSC_TRUE;
1988: B->symmetric = PETSC_BOOL3_TRUE;
1989: B->structurally_symmetric = PETSC_BOOL3_TRUE;
1990: #if !PetscDefined(USE_COMPLEX)
1991: B->hermitian = PETSC_BOOL3_TRUE;
1992: #endif
1994: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQSBAIJ));
1996: PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "Options for SEQSBAIJ matrix", "Mat");
1997: PetscCall(PetscOptionsBool("-mat_no_unroll", "Do not optimize for inodes (slower)", NULL, no_unroll, &no_unroll, NULL));
1998: if (no_unroll) PetscCall(PetscInfo(B, "Not using Inode routines due to -mat_no_unroll\n"));
1999: PetscCall(PetscOptionsBool("-mat_no_inode", "Do not optimize for inodes (slower)", NULL, no_inode, &no_inode, NULL));
2000: if (no_inode) PetscCall(PetscInfo(B, "Not using Inode routines due to -mat_no_inode\n"));
2001: PetscCall(PetscOptionsInt("-mat_inode_limit", "Do not use inodes larger than this value", NULL, b->inode.limit, &b->inode.limit, NULL));
2002: PetscOptionsEnd();
2003: b->inode.use = (PetscBool)(!(no_unroll || no_inode));
2004: if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit;
2005: PetscFunctionReturn(PETSC_SUCCESS);
2006: }
2008: /*@
2009: MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
2010: compressed row) `MATSEQSBAIJ` format. For good matrix assembly performance the
2011: user should preallocate the matrix storage by setting the parameter `nz`
2012: (or the array `nnz`).
2014: Collective
2016: Input Parameters:
2017: + B - the symmetric matrix
2018: . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
2019: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
2020: . nz - number of block nonzeros per block row (same for all rows)
2021: - nnz - array containing the number of block nonzeros in the upper triangular plus
2022: diagonal portion of each block (possibly different for each block row) or `NULL`
2024: Options Database Keys:
2025: + -mat_no_unroll - uses code that does not unroll the loops in the block calculations (much slower)
2026: - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in
2028: Level: intermediate
2030: Notes:
2031: Specify the preallocated storage with either `nz` or `nnz` (not both).
2032: Set `nz` = `PETSC_DEFAULT` and `nnz` = `NULL` for PETSc to control dynamic memory
2033: allocation. See [Sparse Matrices](sec_matsparse) for details.
2035: You can call `MatGetInfo()` to get information on how effective the preallocation was;
2036: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
2037: You can also run with the option `-info` and look for messages with the string
2038: malloc in them to see if additional memory allocation was needed.
2040: If the `nnz` parameter is given then the `nz` parameter is ignored
2042: .seealso: [](ch_matrices), `Mat`, [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()`
2043: @*/
2044: PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[])
2045: {
2046: PetscFunctionBegin;
2050: PetscTryMethod(B, "MatSeqSBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[]), (B, bs, nz, nnz));
2051: PetscFunctionReturn(PETSC_SUCCESS);
2052: }
2054: /*@
2055: MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATSEQSBAIJ` format using the given nonzero structure and (optional) numerical values
2057: Input Parameters:
2058: + B - the matrix
2059: . bs - size of block, the blocks are ALWAYS square.
2060: . i - the indices into `j` for the start of each local row (indices start with zero)
2061: . j - the column indices for each local row (indices start with zero) these must be sorted for each row
2062: - v - optional values in the matrix, use `NULL` if not provided
2064: Level: advanced
2066: Notes:
2067: The `i`,`j`,`v` values are COPIED with this routine; to avoid the copy use `MatCreateSeqSBAIJWithArrays()`
2069: The order of the entries in values is specified by the `MatOption` `MAT_ROW_ORIENTED`. For example, C programs
2070: may want to use the default `MAT_ROW_ORIENTED` = `PETSC_TRUE` and use an array v[nnz][bs][bs] where the second index is
2071: over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
2072: `MAT_ROW_ORIENTED` = `PETSC_FALSE` and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2073: block column and the second index is over columns within a block.
2075: Any entries provided that lie below the diagonal are ignored
2077: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries
2078: and usually the numerical values as well
2080: .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValuesBlocked()`, `MatSeqSBAIJSetPreallocation()`
2081: @*/
2082: PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[])
2083: {
2084: PetscFunctionBegin;
2088: PetscTryMethod(B, "MatSeqSBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v));
2089: PetscFunctionReturn(PETSC_SUCCESS);
2090: }
2092: /*@
2093: MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in (block
2094: compressed row) `MATSEQSBAIJ` format. For good matrix assembly performance the
2095: user should preallocate the matrix storage by setting the parameter `nz`
2096: (or the array `nnz`).
2098: Collective
2100: Input Parameters:
2101: + comm - MPI communicator, set to `PETSC_COMM_SELF`
2102: . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
2103: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2104: . m - number of rows
2105: . n - number of columns
2106: . nz - number of block nonzeros per block row (same for all rows)
2107: - nnz - array containing the number of block nonzeros in the upper triangular plus
2108: diagonal portion of each block (possibly different for each block row) or `NULL`
2110: Output Parameter:
2111: . A - the symmetric matrix
2113: Options Database Keys:
2114: + -mat_no_unroll - uses code that does not unroll the loops in the block calculations (much slower)
2115: - -mat_block_size - size of the blocks to use
2117: Level: intermediate
2119: Notes:
2120: It is recommended that one use `MatCreateFromOptions()` or the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
2121: MatXXXXSetPreallocation() paradigm instead of this routine directly.
2122: [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`]
2124: The number of rows and columns must be divisible by blocksize.
2125: This matrix type does not support complex Hermitian operation.
2127: Specify the preallocated storage with either `nz` or `nnz` (not both).
2128: Set `nz` = `PETSC_DEFAULT` and `nnz` = `NULL` for PETSc to control dynamic memory
2129: allocation. See [Sparse Matrices](sec_matsparse) for details.
2131: If the `nnz` parameter is given then the `nz` parameter is ignored
2133: .seealso: [](ch_matrices), `Mat`, [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()`
2134: @*/
2135: PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A)
2136: {
2137: PetscFunctionBegin;
2138: PetscCall(MatCreate(comm, A));
2139: PetscCall(MatSetSizes(*A, m, n, m, n));
2140: PetscCall(MatSetType(*A, MATSEQSBAIJ));
2141: PetscCall(MatSeqSBAIJSetPreallocation(*A, bs, nz, (PetscInt *)nnz));
2142: PetscFunctionReturn(PETSC_SUCCESS);
2143: }
2145: PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B)
2146: {
2147: Mat C;
2148: Mat_SeqSBAIJ *c, *a = (Mat_SeqSBAIJ *)A->data;
2149: PetscInt i, mbs = a->mbs, nz = a->nz, bs2 = a->bs2;
2151: PetscFunctionBegin;
2152: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Cannot duplicate unassembled matrix");
2153: PetscCheck(a->i[mbs] == nz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt matrix");
2155: *B = NULL;
2156: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
2157: PetscCall(MatSetSizes(C, A->rmap->N, A->cmap->n, A->rmap->N, A->cmap->n));
2158: PetscCall(MatSetBlockSizesFromMats(C, A, A));
2159: PetscCall(MatSetType(C, MATSEQSBAIJ));
2160: PetscCall(MatSetOption(C, MAT_STRUCTURE_ONLY, A->structure_only));
2161: c = (Mat_SeqSBAIJ *)C->data;
2163: C->preallocated = PETSC_TRUE;
2164: C->factortype = A->factortype;
2165: c->row = NULL;
2166: c->icol = NULL;
2167: c->saved_values = NULL;
2168: c->keepnonzeropattern = a->keepnonzeropattern;
2169: C->assembled = PETSC_TRUE;
2171: PetscCall(PetscLayoutReference(A->rmap, &C->rmap));
2172: PetscCall(PetscLayoutReference(A->cmap, &C->cmap));
2173: c->bs2 = a->bs2;
2174: c->mbs = a->mbs;
2175: c->nbs = a->nbs;
2177: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2178: c->imax = a->imax;
2179: c->ilen = a->ilen;
2180: c->free_imax_ilen = PETSC_FALSE;
2181: } else {
2182: PetscCall(PetscMalloc2(mbs + 1, &c->imax, mbs + 1, &c->ilen));
2183: for (i = 0; i < mbs; i++) {
2184: c->imax[i] = a->imax[i];
2185: c->ilen[i] = a->ilen[i];
2186: }
2187: c->free_imax_ilen = PETSC_TRUE;
2188: }
2190: /* allocate the matrix space */
2191: if (!A->structure_only) PetscCall(PetscShmgetAllocateArray(bs2 * nz, sizeof(PetscScalar), (void **)&c->a));
2192: c->free_a = PETSC_TRUE;
2193: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2194: if (!A->structure_only) PetscCall(PetscArrayzero(c->a, bs2 * nz));
2195: c->i = a->i;
2196: c->j = a->j;
2197: c->free_ij = PETSC_FALSE;
2198: c->parent = A;
2199: PetscCall(PetscObjectReference((PetscObject)A));
2200: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
2201: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
2202: } else {
2203: PetscCall(PetscShmgetAllocateArray(nz, sizeof(PetscInt), (void **)&c->j));
2204: PetscCall(PetscShmgetAllocateArray(mbs + 1, sizeof(PetscInt), (void **)&c->i));
2205: PetscCall(PetscArraycpy(c->i, a->i, mbs + 1));
2206: c->free_ij = PETSC_TRUE;
2207: }
2208: if (mbs > 0) {
2209: if (cpvalues != MAT_SHARE_NONZERO_PATTERN) PetscCall(PetscArraycpy(c->j, a->j, nz));
2210: if (!A->structure_only) {
2211: if (cpvalues == MAT_COPY_VALUES) PetscCall(PetscArraycpy(c->a, a->a, bs2 * nz));
2212: else PetscCall(PetscArrayzero(c->a, bs2 * nz));
2213: }
2214: if (a->jshort) {
2215: /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */
2216: /* if the parent matrix is reassembled, this child matrix will never notice */
2217: PetscCall(PetscMalloc1(nz, &c->jshort));
2218: PetscCall(PetscArraycpy(c->jshort, a->jshort, nz));
2220: c->free_jshort = PETSC_TRUE;
2221: }
2222: }
2224: c->roworiented = a->roworiented;
2225: c->nonew = a->nonew;
2226: c->nz = a->nz;
2227: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
2228: c->solve_work = NULL;
2229: c->mult_work = NULL;
2231: *B = C;
2232: PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist));
2233: PetscFunctionReturn(PETSC_SUCCESS);
2234: }
2236: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
2237: #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary
2239: PetscErrorCode MatLoad_SeqSBAIJ(Mat mat, PetscViewer viewer)
2240: {
2241: PetscBool isbinary;
2243: PetscFunctionBegin;
2244: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2245: 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);
2246: PetscCall(MatLoad_SeqSBAIJ_Binary(mat, viewer));
2247: PetscFunctionReturn(PETSC_SUCCESS);
2248: }
2250: /*@
2251: MatCreateSeqSBAIJWithArrays - Creates an sequential `MATSEQSBAIJ` matrix using matrix elements
2252: (upper triangular entries in CSR format) provided by the user.
2254: Collective
2256: Input Parameters:
2257: + comm - must be an MPI communicator of size 1
2258: . bs - size of block
2259: . m - number of rows
2260: . n - number of columns
2261: . 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
2262: . j - column indices
2263: - a - matrix values
2265: Output Parameter:
2266: . mat - the matrix
2268: Level: advanced
2270: Notes:
2271: The `i`, `j`, and `a` arrays are not copied by this routine, the user must free these arrays
2272: once the matrix is destroyed
2274: You cannot set new nonzero locations into this matrix, that will generate an error.
2276: The `i` and `j` indices are 0 based
2278: When block size is greater than 1 the matrix values must be stored using the `MATSBAIJ` storage format. For block size of 1
2279: it is the regular CSR format excluding the lower triangular elements.
2281: .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSBAIJ()`, `MatCreateSeqSBAIJ()`
2282: @*/
2283: PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat)
2284: {
2285: Mat_SeqSBAIJ *sbaij;
2287: PetscFunctionBegin;
2288: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "block size %" PetscInt_FMT " > 1 is not supported yet", bs);
2289: PetscCheck(m == 0 || i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");
2291: PetscCall(MatCreate(comm, mat));
2292: PetscCall(MatSetSizes(*mat, m, n, m, n));
2293: PetscCall(MatSetType(*mat, MATSEQSBAIJ));
2294: PetscCall(MatSeqSBAIJSetPreallocation(*mat, bs, MAT_SKIP_ALLOCATION, NULL));
2295: sbaij = (Mat_SeqSBAIJ *)(*mat)->data;
2296: PetscCall(PetscMalloc2(m, &sbaij->imax, m, &sbaij->ilen));
2298: sbaij->i = i;
2299: sbaij->j = j;
2300: sbaij->a = a;
2302: sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2303: sbaij->free_a = PETSC_FALSE;
2304: sbaij->free_ij = PETSC_FALSE;
2305: sbaij->free_imax_ilen = PETSC_TRUE;
2307: for (PetscInt ii = 0; ii < m; ii++) {
2308: sbaij->ilen[ii] = sbaij->imax[ii] = i[ii + 1] - i[ii];
2309: PetscCheck(i[ii + 1] >= i[ii], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT, ii, i[ii + 1] - i[ii]);
2310: }
2311: if (PetscDefined(USE_DEBUG)) {
2312: for (PetscInt ii = 0; ii < sbaij->i[m]; ii++) {
2313: PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]);
2314: PetscCheck(j[ii] < n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column index too large at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]);
2315: }
2316: }
2318: PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
2319: PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
2320: PetscFunctionReturn(PETSC_SUCCESS);
2321: }
2323: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
2324: {
2325: PetscFunctionBegin;
2326: PetscCall(MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm, inmat, n, scall, outmat));
2327: PetscFunctionReturn(PETSC_SUCCESS);
2328: }