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