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