Actual source code: baij.c

  1: /*
  2:     Defines the basic matrix operations for the BAIJ (compressed row)
  3:   matrix storage format.
  4: */
  5: #include <../src/mat/impls/baij/seq/baij.h>
  6: #include <petscblaslapack.h>
  7: #include <petsc/private/kernels/blockinvert.h>
  8: #include <petsc/private/kernels/blockmatmult.h>

 10: /* defines MatSetValues_Seq_Hash(), MatAssemblyEnd_Seq_Hash(), MatSetUp_Seq_Hash() */
 11: #define TYPE BAIJ
 12: #define TYPE_BS
 13: #include "../src/mat/impls/aij/seq/seqhashmatsetvalues.h"
 14: #undef TYPE_BS
 15: #define TYPE_BS _BS
 16: #define TYPE_BS_ON
 17: #include "../src/mat/impls/aij/seq/seqhashmatsetvalues.h"
 18: #undef TYPE_BS
 19: #include "../src/mat/impls/aij/seq/seqhashmat.h"
 20: #undef TYPE
 21: #undef TYPE_BS_ON

 23: #if defined(PETSC_HAVE_HYPRE)
 24: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat, MatType, MatReuse, Mat *);
 25: #endif

 27: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
 28: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqBAIJMKL(Mat, MatType, MatReuse, Mat *);
 29: #endif
 30: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat, MatType, MatReuse, Mat *);

 32: static PetscErrorCode MatGetColumnReductions_SeqBAIJ(Mat A, PetscInt type, PetscReal *reductions)
 33: {
 34:   Mat_SeqBAIJ *a_aij = (Mat_SeqBAIJ *)A->data;
 35:   PetscInt     m, n, ib, jb, bs = A->rmap->bs;
 36:   MatScalar   *a_val = a_aij->a;

 38:   PetscFunctionBegin;
 39:   PetscCall(MatGetSize(A, &m, &n));
 40:   PetscCall(PetscArrayzero(reductions, n));
 41:   if (type == NORM_2) {
 42:     for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
 43:       for (jb = 0; jb < bs; jb++) {
 44:         for (ib = 0; ib < bs; ib++) {
 45:           reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val * *a_val);
 46:           a_val++;
 47:         }
 48:       }
 49:     }
 50:   } else if (type == NORM_1) {
 51:     for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
 52:       for (jb = 0; jb < bs; jb++) {
 53:         for (ib = 0; ib < bs; ib++) {
 54:           reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val);
 55:           a_val++;
 56:         }
 57:       }
 58:     }
 59:   } else if (type == NORM_INFINITY) {
 60:     for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
 61:       for (jb = 0; jb < bs; jb++) {
 62:         for (ib = 0; ib < bs; ib++) {
 63:           int col         = A->cmap->rstart + a_aij->j[i] * bs + jb;
 64:           reductions[col] = PetscMax(PetscAbsScalar(*a_val), reductions[col]);
 65:           a_val++;
 66:         }
 67:       }
 68:     }
 69:   } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
 70:     for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
 71:       for (jb = 0; jb < bs; jb++) {
 72:         for (ib = 0; ib < bs; ib++) {
 73:           reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscRealPart(*a_val);
 74:           a_val++;
 75:         }
 76:       }
 77:     }
 78:   } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
 79:     for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
 80:       for (jb = 0; jb < bs; jb++) {
 81:         for (ib = 0; ib < bs; ib++) {
 82:           reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscImaginaryPart(*a_val);
 83:           a_val++;
 84:         }
 85:       }
 86:     }
 87:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
 88:   if (type == NORM_2) {
 89:     for (PetscInt i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
 90:   } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
 91:     for (PetscInt i = 0; i < n; i++) reductions[i] /= m;
 92:   }
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

 96: static PetscErrorCode MatInvertBlockDiagonal_SeqBAIJ(Mat A, const PetscScalar **values)
 97: {
 98:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
 99:   PetscInt    *diag_offset, i, bs = A->rmap->bs, mbs = a->mbs, ipvt[5], bs2 = bs * bs, *v_pivots;
100:   MatScalar   *v     = a->a, *odiag, *diag, work[25], *v_work;
101:   PetscReal    shift = 0.0;
102:   PetscBool    allowzeropivot, zeropivotdetected = PETSC_FALSE;

104:   PetscFunctionBegin;
105:   allowzeropivot = PetscNot(A->erroriffailure);

107:   if (a->idiagvalid) {
108:     if (values) *values = a->idiag;
109:     PetscFunctionReturn(PETSC_SUCCESS);
110:   }
111:   PetscCall(MatMarkDiagonal_SeqBAIJ(A));
112:   diag_offset = a->diag;
113:   if (!a->idiag) { PetscCall(PetscMalloc1(bs2 * mbs, &a->idiag)); }
114:   diag = a->idiag;
115:   if (values) *values = a->idiag;
116:   /* factor and invert each block */
117:   switch (bs) {
118:   case 1:
119:     for (i = 0; i < mbs; i++) {
120:       odiag   = v + 1 * diag_offset[i];
121:       diag[0] = odiag[0];

123:       if (PetscAbsScalar(diag[0] + shift) < PETSC_MACHINE_EPSILON) {
124:         if (allowzeropivot) {
125:           A->factorerrortype             = MAT_FACTOR_NUMERIC_ZEROPIVOT;
126:           A->factorerror_zeropivot_value = PetscAbsScalar(diag[0]);
127:           A->factorerror_zeropivot_row   = i;
128:           PetscCall(PetscInfo(A, "Zero pivot, row %" PetscInt_FMT "\n", i));
129:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT " pivot value %g tolerance %g", i, (double)PetscAbsScalar(diag[0]), (double)PETSC_MACHINE_EPSILON);
130:       }

132:       diag[0] = (PetscScalar)1.0 / (diag[0] + shift);
133:       diag += 1;
134:     }
135:     break;
136:   case 2:
137:     for (i = 0; i < mbs; i++) {
138:       odiag   = v + 4 * diag_offset[i];
139:       diag[0] = odiag[0];
140:       diag[1] = odiag[1];
141:       diag[2] = odiag[2];
142:       diag[3] = odiag[3];
143:       PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected));
144:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
145:       diag += 4;
146:     }
147:     break;
148:   case 3:
149:     for (i = 0; i < mbs; i++) {
150:       odiag   = v + 9 * diag_offset[i];
151:       diag[0] = odiag[0];
152:       diag[1] = odiag[1];
153:       diag[2] = odiag[2];
154:       diag[3] = odiag[3];
155:       diag[4] = odiag[4];
156:       diag[5] = odiag[5];
157:       diag[6] = odiag[6];
158:       diag[7] = odiag[7];
159:       diag[8] = odiag[8];
160:       PetscCall(PetscKernel_A_gets_inverse_A_3(diag, shift, allowzeropivot, &zeropivotdetected));
161:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
162:       diag += 9;
163:     }
164:     break;
165:   case 4:
166:     for (i = 0; i < mbs; i++) {
167:       odiag = v + 16 * diag_offset[i];
168:       PetscCall(PetscArraycpy(diag, odiag, 16));
169:       PetscCall(PetscKernel_A_gets_inverse_A_4(diag, shift, allowzeropivot, &zeropivotdetected));
170:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
171:       diag += 16;
172:     }
173:     break;
174:   case 5:
175:     for (i = 0; i < mbs; i++) {
176:       odiag = v + 25 * diag_offset[i];
177:       PetscCall(PetscArraycpy(diag, odiag, 25));
178:       PetscCall(PetscKernel_A_gets_inverse_A_5(diag, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
179:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
180:       diag += 25;
181:     }
182:     break;
183:   case 6:
184:     for (i = 0; i < mbs; i++) {
185:       odiag = v + 36 * diag_offset[i];
186:       PetscCall(PetscArraycpy(diag, odiag, 36));
187:       PetscCall(PetscKernel_A_gets_inverse_A_6(diag, shift, allowzeropivot, &zeropivotdetected));
188:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
189:       diag += 36;
190:     }
191:     break;
192:   case 7:
193:     for (i = 0; i < mbs; i++) {
194:       odiag = v + 49 * diag_offset[i];
195:       PetscCall(PetscArraycpy(diag, odiag, 49));
196:       PetscCall(PetscKernel_A_gets_inverse_A_7(diag, shift, allowzeropivot, &zeropivotdetected));
197:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
198:       diag += 49;
199:     }
200:     break;
201:   default:
202:     PetscCall(PetscMalloc2(bs, &v_work, bs, &v_pivots));
203:     for (i = 0; i < mbs; i++) {
204:       odiag = v + bs2 * diag_offset[i];
205:       PetscCall(PetscArraycpy(diag, odiag, bs2));
206:       PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
207:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
208:       diag += bs2;
209:     }
210:     PetscCall(PetscFree2(v_work, v_pivots));
211:   }
212:   a->idiagvalid = PETSC_TRUE;
213:   PetscFunctionReturn(PETSC_SUCCESS);
214: }

216: static PetscErrorCode MatSOR_SeqBAIJ(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
217: {
218:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data;
219:   PetscScalar       *x, *work, *w, *workt, *t;
220:   const MatScalar   *v, *aa = a->a, *idiag;
221:   const PetscScalar *b, *xb;
222:   PetscScalar        s[7], xw[7] = {0}; /* avoid some compilers thinking xw is uninitialized */
223:   PetscInt           m = a->mbs, i, i2, nz, bs = A->rmap->bs, bs2 = bs * bs, k, j, idx, it;
224:   const PetscInt    *diag, *ai = a->i, *aj = a->j, *vi;

226:   PetscFunctionBegin;
227:   its = its * lits;
228:   PetscCheck(!(flag & SOR_EISENSTAT), PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for Eisenstat");
229:   PetscCheck(its > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Relaxation requires global its %" PetscInt_FMT " and local its %" PetscInt_FMT " both positive", its, lits);
230:   PetscCheck(!fshift, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for diagonal shift");
231:   PetscCheck(omega == 1.0, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for non-trivial relaxation factor");
232:   PetscCheck(!(flag & SOR_APPLY_UPPER) && !(flag & SOR_APPLY_LOWER), PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for applying upper or lower triangular parts");

234:   if (!a->idiagvalid) PetscCall(MatInvertBlockDiagonal(A, NULL));

236:   if (!m) PetscFunctionReturn(PETSC_SUCCESS);
237:   diag  = a->diag;
238:   idiag = a->idiag;
239:   k     = PetscMax(A->rmap->n, A->cmap->n);
240:   if (!a->mult_work) PetscCall(PetscMalloc1(k + 1, &a->mult_work));
241:   if (!a->sor_workt) PetscCall(PetscMalloc1(k, &a->sor_workt));
242:   if (!a->sor_work) PetscCall(PetscMalloc1(bs, &a->sor_work));
243:   work = a->mult_work;
244:   t    = a->sor_workt;
245:   w    = a->sor_work;

247:   PetscCall(VecGetArray(xx, &x));
248:   PetscCall(VecGetArrayRead(bb, &b));

250:   if (flag & SOR_ZERO_INITIAL_GUESS) {
251:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
252:       switch (bs) {
253:       case 1:
254:         PetscKernel_v_gets_A_times_w_1(x, idiag, b);
255:         t[0] = b[0];
256:         i2   = 1;
257:         idiag += 1;
258:         for (i = 1; i < m; i++) {
259:           v    = aa + ai[i];
260:           vi   = aj + ai[i];
261:           nz   = diag[i] - ai[i];
262:           s[0] = b[i2];
263:           for (j = 0; j < nz; j++) {
264:             xw[0] = x[vi[j]];
265:             PetscKernel_v_gets_v_minus_A_times_w_1(s, (v + j), xw);
266:           }
267:           t[i2] = s[0];
268:           PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
269:           x[i2] = xw[0];
270:           idiag += 1;
271:           i2 += 1;
272:         }
273:         break;
274:       case 2:
275:         PetscKernel_v_gets_A_times_w_2(x, idiag, b);
276:         t[0] = b[0];
277:         t[1] = b[1];
278:         i2   = 2;
279:         idiag += 4;
280:         for (i = 1; i < m; i++) {
281:           v    = aa + 4 * ai[i];
282:           vi   = aj + ai[i];
283:           nz   = diag[i] - ai[i];
284:           s[0] = b[i2];
285:           s[1] = b[i2 + 1];
286:           for (j = 0; j < nz; j++) {
287:             idx   = 2 * vi[j];
288:             it    = 4 * j;
289:             xw[0] = x[idx];
290:             xw[1] = x[1 + idx];
291:             PetscKernel_v_gets_v_minus_A_times_w_2(s, (v + it), xw);
292:           }
293:           t[i2]     = s[0];
294:           t[i2 + 1] = s[1];
295:           PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
296:           x[i2]     = xw[0];
297:           x[i2 + 1] = xw[1];
298:           idiag += 4;
299:           i2 += 2;
300:         }
301:         break;
302:       case 3:
303:         PetscKernel_v_gets_A_times_w_3(x, idiag, b);
304:         t[0] = b[0];
305:         t[1] = b[1];
306:         t[2] = b[2];
307:         i2   = 3;
308:         idiag += 9;
309:         for (i = 1; i < m; i++) {
310:           v    = aa + 9 * ai[i];
311:           vi   = aj + ai[i];
312:           nz   = diag[i] - ai[i];
313:           s[0] = b[i2];
314:           s[1] = b[i2 + 1];
315:           s[2] = b[i2 + 2];
316:           while (nz--) {
317:             idx   = 3 * (*vi++);
318:             xw[0] = x[idx];
319:             xw[1] = x[1 + idx];
320:             xw[2] = x[2 + idx];
321:             PetscKernel_v_gets_v_minus_A_times_w_3(s, v, xw);
322:             v += 9;
323:           }
324:           t[i2]     = s[0];
325:           t[i2 + 1] = s[1];
326:           t[i2 + 2] = s[2];
327:           PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
328:           x[i2]     = xw[0];
329:           x[i2 + 1] = xw[1];
330:           x[i2 + 2] = xw[2];
331:           idiag += 9;
332:           i2 += 3;
333:         }
334:         break;
335:       case 4:
336:         PetscKernel_v_gets_A_times_w_4(x, idiag, b);
337:         t[0] = b[0];
338:         t[1] = b[1];
339:         t[2] = b[2];
340:         t[3] = b[3];
341:         i2   = 4;
342:         idiag += 16;
343:         for (i = 1; i < m; i++) {
344:           v    = aa + 16 * ai[i];
345:           vi   = aj + ai[i];
346:           nz   = diag[i] - ai[i];
347:           s[0] = b[i2];
348:           s[1] = b[i2 + 1];
349:           s[2] = b[i2 + 2];
350:           s[3] = b[i2 + 3];
351:           while (nz--) {
352:             idx   = 4 * (*vi++);
353:             xw[0] = x[idx];
354:             xw[1] = x[1 + idx];
355:             xw[2] = x[2 + idx];
356:             xw[3] = x[3 + idx];
357:             PetscKernel_v_gets_v_minus_A_times_w_4(s, v, xw);
358:             v += 16;
359:           }
360:           t[i2]     = s[0];
361:           t[i2 + 1] = s[1];
362:           t[i2 + 2] = s[2];
363:           t[i2 + 3] = s[3];
364:           PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
365:           x[i2]     = xw[0];
366:           x[i2 + 1] = xw[1];
367:           x[i2 + 2] = xw[2];
368:           x[i2 + 3] = xw[3];
369:           idiag += 16;
370:           i2 += 4;
371:         }
372:         break;
373:       case 5:
374:         PetscKernel_v_gets_A_times_w_5(x, idiag, b);
375:         t[0] = b[0];
376:         t[1] = b[1];
377:         t[2] = b[2];
378:         t[3] = b[3];
379:         t[4] = b[4];
380:         i2   = 5;
381:         idiag += 25;
382:         for (i = 1; i < m; i++) {
383:           v    = aa + 25 * ai[i];
384:           vi   = aj + ai[i];
385:           nz   = diag[i] - ai[i];
386:           s[0] = b[i2];
387:           s[1] = b[i2 + 1];
388:           s[2] = b[i2 + 2];
389:           s[3] = b[i2 + 3];
390:           s[4] = b[i2 + 4];
391:           while (nz--) {
392:             idx   = 5 * (*vi++);
393:             xw[0] = x[idx];
394:             xw[1] = x[1 + idx];
395:             xw[2] = x[2 + idx];
396:             xw[3] = x[3 + idx];
397:             xw[4] = x[4 + idx];
398:             PetscKernel_v_gets_v_minus_A_times_w_5(s, v, xw);
399:             v += 25;
400:           }
401:           t[i2]     = s[0];
402:           t[i2 + 1] = s[1];
403:           t[i2 + 2] = s[2];
404:           t[i2 + 3] = s[3];
405:           t[i2 + 4] = s[4];
406:           PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
407:           x[i2]     = xw[0];
408:           x[i2 + 1] = xw[1];
409:           x[i2 + 2] = xw[2];
410:           x[i2 + 3] = xw[3];
411:           x[i2 + 4] = xw[4];
412:           idiag += 25;
413:           i2 += 5;
414:         }
415:         break;
416:       case 6:
417:         PetscKernel_v_gets_A_times_w_6(x, idiag, b);
418:         t[0] = b[0];
419:         t[1] = b[1];
420:         t[2] = b[2];
421:         t[3] = b[3];
422:         t[4] = b[4];
423:         t[5] = b[5];
424:         i2   = 6;
425:         idiag += 36;
426:         for (i = 1; i < m; i++) {
427:           v    = aa + 36 * ai[i];
428:           vi   = aj + ai[i];
429:           nz   = diag[i] - ai[i];
430:           s[0] = b[i2];
431:           s[1] = b[i2 + 1];
432:           s[2] = b[i2 + 2];
433:           s[3] = b[i2 + 3];
434:           s[4] = b[i2 + 4];
435:           s[5] = b[i2 + 5];
436:           while (nz--) {
437:             idx   = 6 * (*vi++);
438:             xw[0] = x[idx];
439:             xw[1] = x[1 + idx];
440:             xw[2] = x[2 + idx];
441:             xw[3] = x[3 + idx];
442:             xw[4] = x[4 + idx];
443:             xw[5] = x[5 + idx];
444:             PetscKernel_v_gets_v_minus_A_times_w_6(s, v, xw);
445:             v += 36;
446:           }
447:           t[i2]     = s[0];
448:           t[i2 + 1] = s[1];
449:           t[i2 + 2] = s[2];
450:           t[i2 + 3] = s[3];
451:           t[i2 + 4] = s[4];
452:           t[i2 + 5] = s[5];
453:           PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
454:           x[i2]     = xw[0];
455:           x[i2 + 1] = xw[1];
456:           x[i2 + 2] = xw[2];
457:           x[i2 + 3] = xw[3];
458:           x[i2 + 4] = xw[4];
459:           x[i2 + 5] = xw[5];
460:           idiag += 36;
461:           i2 += 6;
462:         }
463:         break;
464:       case 7:
465:         PetscKernel_v_gets_A_times_w_7(x, idiag, b);
466:         t[0] = b[0];
467:         t[1] = b[1];
468:         t[2] = b[2];
469:         t[3] = b[3];
470:         t[4] = b[4];
471:         t[5] = b[5];
472:         t[6] = b[6];
473:         i2   = 7;
474:         idiag += 49;
475:         for (i = 1; i < m; i++) {
476:           v    = aa + 49 * ai[i];
477:           vi   = aj + ai[i];
478:           nz   = diag[i] - ai[i];
479:           s[0] = b[i2];
480:           s[1] = b[i2 + 1];
481:           s[2] = b[i2 + 2];
482:           s[3] = b[i2 + 3];
483:           s[4] = b[i2 + 4];
484:           s[5] = b[i2 + 5];
485:           s[6] = b[i2 + 6];
486:           while (nz--) {
487:             idx   = 7 * (*vi++);
488:             xw[0] = x[idx];
489:             xw[1] = x[1 + idx];
490:             xw[2] = x[2 + idx];
491:             xw[3] = x[3 + idx];
492:             xw[4] = x[4 + idx];
493:             xw[5] = x[5 + idx];
494:             xw[6] = x[6 + idx];
495:             PetscKernel_v_gets_v_minus_A_times_w_7(s, v, xw);
496:             v += 49;
497:           }
498:           t[i2]     = s[0];
499:           t[i2 + 1] = s[1];
500:           t[i2 + 2] = s[2];
501:           t[i2 + 3] = s[3];
502:           t[i2 + 4] = s[4];
503:           t[i2 + 5] = s[5];
504:           t[i2 + 6] = s[6];
505:           PetscKernel_v_gets_A_times_w_7(xw, idiag, s);
506:           x[i2]     = xw[0];
507:           x[i2 + 1] = xw[1];
508:           x[i2 + 2] = xw[2];
509:           x[i2 + 3] = xw[3];
510:           x[i2 + 4] = xw[4];
511:           x[i2 + 5] = xw[5];
512:           x[i2 + 6] = xw[6];
513:           idiag += 49;
514:           i2 += 7;
515:         }
516:         break;
517:       default:
518:         PetscKernel_w_gets_Ar_times_v(bs, bs, b, idiag, x);
519:         PetscCall(PetscArraycpy(t, b, bs));
520:         i2 = bs;
521:         idiag += bs2;
522:         for (i = 1; i < m; i++) {
523:           v  = aa + bs2 * ai[i];
524:           vi = aj + ai[i];
525:           nz = diag[i] - ai[i];

527:           PetscCall(PetscArraycpy(w, b + i2, bs));
528:           /* copy all rows of x that are needed into contiguous space */
529:           workt = work;
530:           for (j = 0; j < nz; j++) {
531:             PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
532:             workt += bs;
533:           }
534:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, v, work);
535:           PetscCall(PetscArraycpy(t + i2, w, bs));
536:           PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, x + i2);

538:           idiag += bs2;
539:           i2 += bs;
540:         }
541:         break;
542:       }
543:       /* for logging purposes assume number of nonzero in lower half is 1/2 of total */
544:       PetscCall(PetscLogFlops(1.0 * bs2 * a->nz));
545:       xb = t;
546:     } else xb = b;
547:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
548:       idiag = a->idiag + bs2 * (a->mbs - 1);
549:       i2    = bs * (m - 1);
550:       switch (bs) {
551:       case 1:
552:         s[0] = xb[i2];
553:         PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
554:         x[i2] = xw[0];
555:         i2 -= 1;
556:         for (i = m - 2; i >= 0; i--) {
557:           v    = aa + (diag[i] + 1);
558:           vi   = aj + diag[i] + 1;
559:           nz   = ai[i + 1] - diag[i] - 1;
560:           s[0] = xb[i2];
561:           for (j = 0; j < nz; j++) {
562:             xw[0] = x[vi[j]];
563:             PetscKernel_v_gets_v_minus_A_times_w_1(s, (v + j), xw);
564:           }
565:           PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
566:           x[i2] = xw[0];
567:           idiag -= 1;
568:           i2 -= 1;
569:         }
570:         break;
571:       case 2:
572:         s[0] = xb[i2];
573:         s[1] = xb[i2 + 1];
574:         PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
575:         x[i2]     = xw[0];
576:         x[i2 + 1] = xw[1];
577:         i2 -= 2;
578:         idiag -= 4;
579:         for (i = m - 2; i >= 0; i--) {
580:           v    = aa + 4 * (diag[i] + 1);
581:           vi   = aj + diag[i] + 1;
582:           nz   = ai[i + 1] - diag[i] - 1;
583:           s[0] = xb[i2];
584:           s[1] = xb[i2 + 1];
585:           for (j = 0; j < nz; j++) {
586:             idx   = 2 * vi[j];
587:             it    = 4 * j;
588:             xw[0] = x[idx];
589:             xw[1] = x[1 + idx];
590:             PetscKernel_v_gets_v_minus_A_times_w_2(s, (v + it), xw);
591:           }
592:           PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
593:           x[i2]     = xw[0];
594:           x[i2 + 1] = xw[1];
595:           idiag -= 4;
596:           i2 -= 2;
597:         }
598:         break;
599:       case 3:
600:         s[0] = xb[i2];
601:         s[1] = xb[i2 + 1];
602:         s[2] = xb[i2 + 2];
603:         PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
604:         x[i2]     = xw[0];
605:         x[i2 + 1] = xw[1];
606:         x[i2 + 2] = xw[2];
607:         i2 -= 3;
608:         idiag -= 9;
609:         for (i = m - 2; i >= 0; i--) {
610:           v    = aa + 9 * (diag[i] + 1);
611:           vi   = aj + diag[i] + 1;
612:           nz   = ai[i + 1] - diag[i] - 1;
613:           s[0] = xb[i2];
614:           s[1] = xb[i2 + 1];
615:           s[2] = xb[i2 + 2];
616:           while (nz--) {
617:             idx   = 3 * (*vi++);
618:             xw[0] = x[idx];
619:             xw[1] = x[1 + idx];
620:             xw[2] = x[2 + idx];
621:             PetscKernel_v_gets_v_minus_A_times_w_3(s, v, xw);
622:             v += 9;
623:           }
624:           PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
625:           x[i2]     = xw[0];
626:           x[i2 + 1] = xw[1];
627:           x[i2 + 2] = xw[2];
628:           idiag -= 9;
629:           i2 -= 3;
630:         }
631:         break;
632:       case 4:
633:         s[0] = xb[i2];
634:         s[1] = xb[i2 + 1];
635:         s[2] = xb[i2 + 2];
636:         s[3] = xb[i2 + 3];
637:         PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
638:         x[i2]     = xw[0];
639:         x[i2 + 1] = xw[1];
640:         x[i2 + 2] = xw[2];
641:         x[i2 + 3] = xw[3];
642:         i2 -= 4;
643:         idiag -= 16;
644:         for (i = m - 2; i >= 0; i--) {
645:           v    = aa + 16 * (diag[i] + 1);
646:           vi   = aj + diag[i] + 1;
647:           nz   = ai[i + 1] - diag[i] - 1;
648:           s[0] = xb[i2];
649:           s[1] = xb[i2 + 1];
650:           s[2] = xb[i2 + 2];
651:           s[3] = xb[i2 + 3];
652:           while (nz--) {
653:             idx   = 4 * (*vi++);
654:             xw[0] = x[idx];
655:             xw[1] = x[1 + idx];
656:             xw[2] = x[2 + idx];
657:             xw[3] = x[3 + idx];
658:             PetscKernel_v_gets_v_minus_A_times_w_4(s, v, xw);
659:             v += 16;
660:           }
661:           PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
662:           x[i2]     = xw[0];
663:           x[i2 + 1] = xw[1];
664:           x[i2 + 2] = xw[2];
665:           x[i2 + 3] = xw[3];
666:           idiag -= 16;
667:           i2 -= 4;
668:         }
669:         break;
670:       case 5:
671:         s[0] = xb[i2];
672:         s[1] = xb[i2 + 1];
673:         s[2] = xb[i2 + 2];
674:         s[3] = xb[i2 + 3];
675:         s[4] = xb[i2 + 4];
676:         PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
677:         x[i2]     = xw[0];
678:         x[i2 + 1] = xw[1];
679:         x[i2 + 2] = xw[2];
680:         x[i2 + 3] = xw[3];
681:         x[i2 + 4] = xw[4];
682:         i2 -= 5;
683:         idiag -= 25;
684:         for (i = m - 2; i >= 0; i--) {
685:           v    = aa + 25 * (diag[i] + 1);
686:           vi   = aj + diag[i] + 1;
687:           nz   = ai[i + 1] - diag[i] - 1;
688:           s[0] = xb[i2];
689:           s[1] = xb[i2 + 1];
690:           s[2] = xb[i2 + 2];
691:           s[3] = xb[i2 + 3];
692:           s[4] = xb[i2 + 4];
693:           while (nz--) {
694:             idx   = 5 * (*vi++);
695:             xw[0] = x[idx];
696:             xw[1] = x[1 + idx];
697:             xw[2] = x[2 + idx];
698:             xw[3] = x[3 + idx];
699:             xw[4] = x[4 + idx];
700:             PetscKernel_v_gets_v_minus_A_times_w_5(s, v, xw);
701:             v += 25;
702:           }
703:           PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
704:           x[i2]     = xw[0];
705:           x[i2 + 1] = xw[1];
706:           x[i2 + 2] = xw[2];
707:           x[i2 + 3] = xw[3];
708:           x[i2 + 4] = xw[4];
709:           idiag -= 25;
710:           i2 -= 5;
711:         }
712:         break;
713:       case 6:
714:         s[0] = xb[i2];
715:         s[1] = xb[i2 + 1];
716:         s[2] = xb[i2 + 2];
717:         s[3] = xb[i2 + 3];
718:         s[4] = xb[i2 + 4];
719:         s[5] = xb[i2 + 5];
720:         PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
721:         x[i2]     = xw[0];
722:         x[i2 + 1] = xw[1];
723:         x[i2 + 2] = xw[2];
724:         x[i2 + 3] = xw[3];
725:         x[i2 + 4] = xw[4];
726:         x[i2 + 5] = xw[5];
727:         i2 -= 6;
728:         idiag -= 36;
729:         for (i = m - 2; i >= 0; i--) {
730:           v    = aa + 36 * (diag[i] + 1);
731:           vi   = aj + diag[i] + 1;
732:           nz   = ai[i + 1] - diag[i] - 1;
733:           s[0] = xb[i2];
734:           s[1] = xb[i2 + 1];
735:           s[2] = xb[i2 + 2];
736:           s[3] = xb[i2 + 3];
737:           s[4] = xb[i2 + 4];
738:           s[5] = xb[i2 + 5];
739:           while (nz--) {
740:             idx   = 6 * (*vi++);
741:             xw[0] = x[idx];
742:             xw[1] = x[1 + idx];
743:             xw[2] = x[2 + idx];
744:             xw[3] = x[3 + idx];
745:             xw[4] = x[4 + idx];
746:             xw[5] = x[5 + idx];
747:             PetscKernel_v_gets_v_minus_A_times_w_6(s, v, xw);
748:             v += 36;
749:           }
750:           PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
751:           x[i2]     = xw[0];
752:           x[i2 + 1] = xw[1];
753:           x[i2 + 2] = xw[2];
754:           x[i2 + 3] = xw[3];
755:           x[i2 + 4] = xw[4];
756:           x[i2 + 5] = xw[5];
757:           idiag -= 36;
758:           i2 -= 6;
759:         }
760:         break;
761:       case 7:
762:         s[0] = xb[i2];
763:         s[1] = xb[i2 + 1];
764:         s[2] = xb[i2 + 2];
765:         s[3] = xb[i2 + 3];
766:         s[4] = xb[i2 + 4];
767:         s[5] = xb[i2 + 5];
768:         s[6] = xb[i2 + 6];
769:         PetscKernel_v_gets_A_times_w_7(x, idiag, b);
770:         x[i2]     = xw[0];
771:         x[i2 + 1] = xw[1];
772:         x[i2 + 2] = xw[2];
773:         x[i2 + 3] = xw[3];
774:         x[i2 + 4] = xw[4];
775:         x[i2 + 5] = xw[5];
776:         x[i2 + 6] = xw[6];
777:         i2 -= 7;
778:         idiag -= 49;
779:         for (i = m - 2; i >= 0; i--) {
780:           v    = aa + 49 * (diag[i] + 1);
781:           vi   = aj + diag[i] + 1;
782:           nz   = ai[i + 1] - diag[i] - 1;
783:           s[0] = xb[i2];
784:           s[1] = xb[i2 + 1];
785:           s[2] = xb[i2 + 2];
786:           s[3] = xb[i2 + 3];
787:           s[4] = xb[i2 + 4];
788:           s[5] = xb[i2 + 5];
789:           s[6] = xb[i2 + 6];
790:           while (nz--) {
791:             idx   = 7 * (*vi++);
792:             xw[0] = x[idx];
793:             xw[1] = x[1 + idx];
794:             xw[2] = x[2 + idx];
795:             xw[3] = x[3 + idx];
796:             xw[4] = x[4 + idx];
797:             xw[5] = x[5 + idx];
798:             xw[6] = x[6 + idx];
799:             PetscKernel_v_gets_v_minus_A_times_w_7(s, v, xw);
800:             v += 49;
801:           }
802:           PetscKernel_v_gets_A_times_w_7(xw, idiag, s);
803:           x[i2]     = xw[0];
804:           x[i2 + 1] = xw[1];
805:           x[i2 + 2] = xw[2];
806:           x[i2 + 3] = xw[3];
807:           x[i2 + 4] = xw[4];
808:           x[i2 + 5] = xw[5];
809:           x[i2 + 6] = xw[6];
810:           idiag -= 49;
811:           i2 -= 7;
812:         }
813:         break;
814:       default:
815:         PetscCall(PetscArraycpy(w, xb + i2, bs));
816:         PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, x + i2);
817:         i2 -= bs;
818:         idiag -= bs2;
819:         for (i = m - 2; i >= 0; i--) {
820:           v  = aa + bs2 * (diag[i] + 1);
821:           vi = aj + diag[i] + 1;
822:           nz = ai[i + 1] - diag[i] - 1;

824:           PetscCall(PetscArraycpy(w, xb + i2, bs));
825:           /* copy all rows of x that are needed into contiguous space */
826:           workt = work;
827:           for (j = 0; j < nz; j++) {
828:             PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
829:             workt += bs;
830:           }
831:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, v, work);
832:           PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, x + i2);

834:           idiag -= bs2;
835:           i2 -= bs;
836:         }
837:         break;
838:       }
839:       PetscCall(PetscLogFlops(1.0 * bs2 * (a->nz)));
840:     }
841:     its--;
842:   }
843:   while (its--) {
844:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
845:       idiag = a->idiag;
846:       i2    = 0;
847:       switch (bs) {
848:       case 1:
849:         for (i = 0; i < m; i++) {
850:           v    = aa + ai[i];
851:           vi   = aj + ai[i];
852:           nz   = ai[i + 1] - ai[i];
853:           s[0] = b[i2];
854:           for (j = 0; j < nz; j++) {
855:             xw[0] = x[vi[j]];
856:             PetscKernel_v_gets_v_minus_A_times_w_1(s, (v + j), xw);
857:           }
858:           PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
859:           x[i2] += xw[0];
860:           idiag += 1;
861:           i2 += 1;
862:         }
863:         break;
864:       case 2:
865:         for (i = 0; i < m; i++) {
866:           v    = aa + 4 * ai[i];
867:           vi   = aj + ai[i];
868:           nz   = ai[i + 1] - ai[i];
869:           s[0] = b[i2];
870:           s[1] = b[i2 + 1];
871:           for (j = 0; j < nz; j++) {
872:             idx   = 2 * vi[j];
873:             it    = 4 * j;
874:             xw[0] = x[idx];
875:             xw[1] = x[1 + idx];
876:             PetscKernel_v_gets_v_minus_A_times_w_2(s, (v + it), xw);
877:           }
878:           PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
879:           x[i2] += xw[0];
880:           x[i2 + 1] += xw[1];
881:           idiag += 4;
882:           i2 += 2;
883:         }
884:         break;
885:       case 3:
886:         for (i = 0; i < m; i++) {
887:           v    = aa + 9 * ai[i];
888:           vi   = aj + ai[i];
889:           nz   = ai[i + 1] - ai[i];
890:           s[0] = b[i2];
891:           s[1] = b[i2 + 1];
892:           s[2] = b[i2 + 2];
893:           while (nz--) {
894:             idx   = 3 * (*vi++);
895:             xw[0] = x[idx];
896:             xw[1] = x[1 + idx];
897:             xw[2] = x[2 + idx];
898:             PetscKernel_v_gets_v_minus_A_times_w_3(s, v, xw);
899:             v += 9;
900:           }
901:           PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
902:           x[i2] += xw[0];
903:           x[i2 + 1] += xw[1];
904:           x[i2 + 2] += xw[2];
905:           idiag += 9;
906:           i2 += 3;
907:         }
908:         break;
909:       case 4:
910:         for (i = 0; i < m; i++) {
911:           v    = aa + 16 * ai[i];
912:           vi   = aj + ai[i];
913:           nz   = ai[i + 1] - ai[i];
914:           s[0] = b[i2];
915:           s[1] = b[i2 + 1];
916:           s[2] = b[i2 + 2];
917:           s[3] = b[i2 + 3];
918:           while (nz--) {
919:             idx   = 4 * (*vi++);
920:             xw[0] = x[idx];
921:             xw[1] = x[1 + idx];
922:             xw[2] = x[2 + idx];
923:             xw[3] = x[3 + idx];
924:             PetscKernel_v_gets_v_minus_A_times_w_4(s, v, xw);
925:             v += 16;
926:           }
927:           PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
928:           x[i2] += xw[0];
929:           x[i2 + 1] += xw[1];
930:           x[i2 + 2] += xw[2];
931:           x[i2 + 3] += xw[3];
932:           idiag += 16;
933:           i2 += 4;
934:         }
935:         break;
936:       case 5:
937:         for (i = 0; i < m; i++) {
938:           v    = aa + 25 * ai[i];
939:           vi   = aj + ai[i];
940:           nz   = ai[i + 1] - ai[i];
941:           s[0] = b[i2];
942:           s[1] = b[i2 + 1];
943:           s[2] = b[i2 + 2];
944:           s[3] = b[i2 + 3];
945:           s[4] = b[i2 + 4];
946:           while (nz--) {
947:             idx   = 5 * (*vi++);
948:             xw[0] = x[idx];
949:             xw[1] = x[1 + idx];
950:             xw[2] = x[2 + idx];
951:             xw[3] = x[3 + idx];
952:             xw[4] = x[4 + idx];
953:             PetscKernel_v_gets_v_minus_A_times_w_5(s, v, xw);
954:             v += 25;
955:           }
956:           PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
957:           x[i2] += xw[0];
958:           x[i2 + 1] += xw[1];
959:           x[i2 + 2] += xw[2];
960:           x[i2 + 3] += xw[3];
961:           x[i2 + 4] += xw[4];
962:           idiag += 25;
963:           i2 += 5;
964:         }
965:         break;
966:       case 6:
967:         for (i = 0; i < m; i++) {
968:           v    = aa + 36 * ai[i];
969:           vi   = aj + ai[i];
970:           nz   = ai[i + 1] - ai[i];
971:           s[0] = b[i2];
972:           s[1] = b[i2 + 1];
973:           s[2] = b[i2 + 2];
974:           s[3] = b[i2 + 3];
975:           s[4] = b[i2 + 4];
976:           s[5] = b[i2 + 5];
977:           while (nz--) {
978:             idx   = 6 * (*vi++);
979:             xw[0] = x[idx];
980:             xw[1] = x[1 + idx];
981:             xw[2] = x[2 + idx];
982:             xw[3] = x[3 + idx];
983:             xw[4] = x[4 + idx];
984:             xw[5] = x[5 + idx];
985:             PetscKernel_v_gets_v_minus_A_times_w_6(s, v, xw);
986:             v += 36;
987:           }
988:           PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
989:           x[i2] += xw[0];
990:           x[i2 + 1] += xw[1];
991:           x[i2 + 2] += xw[2];
992:           x[i2 + 3] += xw[3];
993:           x[i2 + 4] += xw[4];
994:           x[i2 + 5] += xw[5];
995:           idiag += 36;
996:           i2 += 6;
997:         }
998:         break;
999:       case 7:
1000:         for (i = 0; i < m; i++) {
1001:           v    = aa + 49 * ai[i];
1002:           vi   = aj + ai[i];
1003:           nz   = ai[i + 1] - ai[i];
1004:           s[0] = b[i2];
1005:           s[1] = b[i2 + 1];
1006:           s[2] = b[i2 + 2];
1007:           s[3] = b[i2 + 3];
1008:           s[4] = b[i2 + 4];
1009:           s[5] = b[i2 + 5];
1010:           s[6] = b[i2 + 6];
1011:           while (nz--) {
1012:             idx   = 7 * (*vi++);
1013:             xw[0] = x[idx];
1014:             xw[1] = x[1 + idx];
1015:             xw[2] = x[2 + idx];
1016:             xw[3] = x[3 + idx];
1017:             xw[4] = x[4 + idx];
1018:             xw[5] = x[5 + idx];
1019:             xw[6] = x[6 + idx];
1020:             PetscKernel_v_gets_v_minus_A_times_w_7(s, v, xw);
1021:             v += 49;
1022:           }
1023:           PetscKernel_v_gets_A_times_w_7(xw, idiag, s);
1024:           x[i2] += xw[0];
1025:           x[i2 + 1] += xw[1];
1026:           x[i2 + 2] += xw[2];
1027:           x[i2 + 3] += xw[3];
1028:           x[i2 + 4] += xw[4];
1029:           x[i2 + 5] += xw[5];
1030:           x[i2 + 6] += xw[6];
1031:           idiag += 49;
1032:           i2 += 7;
1033:         }
1034:         break;
1035:       default:
1036:         for (i = 0; i < m; i++) {
1037:           v  = aa + bs2 * ai[i];
1038:           vi = aj + ai[i];
1039:           nz = ai[i + 1] - ai[i];

1041:           PetscCall(PetscArraycpy(w, b + i2, bs));
1042:           /* copy all rows of x that are needed into contiguous space */
1043:           workt = work;
1044:           for (j = 0; j < nz; j++) {
1045:             PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
1046:             workt += bs;
1047:           }
1048:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, v, work);
1049:           PetscKernel_w_gets_w_plus_Ar_times_v(bs, bs, w, idiag, x + i2);

1051:           idiag += bs2;
1052:           i2 += bs;
1053:         }
1054:         break;
1055:       }
1056:       PetscCall(PetscLogFlops(2.0 * bs2 * a->nz));
1057:     }
1058:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1059:       idiag = a->idiag + bs2 * (a->mbs - 1);
1060:       i2    = bs * (m - 1);
1061:       switch (bs) {
1062:       case 1:
1063:         for (i = m - 1; i >= 0; i--) {
1064:           v    = aa + ai[i];
1065:           vi   = aj + ai[i];
1066:           nz   = ai[i + 1] - ai[i];
1067:           s[0] = b[i2];
1068:           for (j = 0; j < nz; j++) {
1069:             xw[0] = x[vi[j]];
1070:             PetscKernel_v_gets_v_minus_A_times_w_1(s, (v + j), xw);
1071:           }
1072:           PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
1073:           x[i2] += xw[0];
1074:           idiag -= 1;
1075:           i2 -= 1;
1076:         }
1077:         break;
1078:       case 2:
1079:         for (i = m - 1; i >= 0; i--) {
1080:           v    = aa + 4 * ai[i];
1081:           vi   = aj + ai[i];
1082:           nz   = ai[i + 1] - ai[i];
1083:           s[0] = b[i2];
1084:           s[1] = b[i2 + 1];
1085:           for (j = 0; j < nz; j++) {
1086:             idx   = 2 * vi[j];
1087:             it    = 4 * j;
1088:             xw[0] = x[idx];
1089:             xw[1] = x[1 + idx];
1090:             PetscKernel_v_gets_v_minus_A_times_w_2(s, (v + it), xw);
1091:           }
1092:           PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
1093:           x[i2] += xw[0];
1094:           x[i2 + 1] += xw[1];
1095:           idiag -= 4;
1096:           i2 -= 2;
1097:         }
1098:         break;
1099:       case 3:
1100:         for (i = m - 1; i >= 0; i--) {
1101:           v    = aa + 9 * ai[i];
1102:           vi   = aj + ai[i];
1103:           nz   = ai[i + 1] - ai[i];
1104:           s[0] = b[i2];
1105:           s[1] = b[i2 + 1];
1106:           s[2] = b[i2 + 2];
1107:           while (nz--) {
1108:             idx   = 3 * (*vi++);
1109:             xw[0] = x[idx];
1110:             xw[1] = x[1 + idx];
1111:             xw[2] = x[2 + idx];
1112:             PetscKernel_v_gets_v_minus_A_times_w_3(s, v, xw);
1113:             v += 9;
1114:           }
1115:           PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
1116:           x[i2] += xw[0];
1117:           x[i2 + 1] += xw[1];
1118:           x[i2 + 2] += xw[2];
1119:           idiag -= 9;
1120:           i2 -= 3;
1121:         }
1122:         break;
1123:       case 4:
1124:         for (i = m - 1; i >= 0; i--) {
1125:           v    = aa + 16 * ai[i];
1126:           vi   = aj + ai[i];
1127:           nz   = ai[i + 1] - ai[i];
1128:           s[0] = b[i2];
1129:           s[1] = b[i2 + 1];
1130:           s[2] = b[i2 + 2];
1131:           s[3] = b[i2 + 3];
1132:           while (nz--) {
1133:             idx   = 4 * (*vi++);
1134:             xw[0] = x[idx];
1135:             xw[1] = x[1 + idx];
1136:             xw[2] = x[2 + idx];
1137:             xw[3] = x[3 + idx];
1138:             PetscKernel_v_gets_v_minus_A_times_w_4(s, v, xw);
1139:             v += 16;
1140:           }
1141:           PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
1142:           x[i2] += xw[0];
1143:           x[i2 + 1] += xw[1];
1144:           x[i2 + 2] += xw[2];
1145:           x[i2 + 3] += xw[3];
1146:           idiag -= 16;
1147:           i2 -= 4;
1148:         }
1149:         break;
1150:       case 5:
1151:         for (i = m - 1; i >= 0; i--) {
1152:           v    = aa + 25 * ai[i];
1153:           vi   = aj + ai[i];
1154:           nz   = ai[i + 1] - ai[i];
1155:           s[0] = b[i2];
1156:           s[1] = b[i2 + 1];
1157:           s[2] = b[i2 + 2];
1158:           s[3] = b[i2 + 3];
1159:           s[4] = b[i2 + 4];
1160:           while (nz--) {
1161:             idx   = 5 * (*vi++);
1162:             xw[0] = x[idx];
1163:             xw[1] = x[1 + idx];
1164:             xw[2] = x[2 + idx];
1165:             xw[3] = x[3 + idx];
1166:             xw[4] = x[4 + idx];
1167:             PetscKernel_v_gets_v_minus_A_times_w_5(s, v, xw);
1168:             v += 25;
1169:           }
1170:           PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
1171:           x[i2] += xw[0];
1172:           x[i2 + 1] += xw[1];
1173:           x[i2 + 2] += xw[2];
1174:           x[i2 + 3] += xw[3];
1175:           x[i2 + 4] += xw[4];
1176:           idiag -= 25;
1177:           i2 -= 5;
1178:         }
1179:         break;
1180:       case 6:
1181:         for (i = m - 1; i >= 0; i--) {
1182:           v    = aa + 36 * ai[i];
1183:           vi   = aj + ai[i];
1184:           nz   = ai[i + 1] - ai[i];
1185:           s[0] = b[i2];
1186:           s[1] = b[i2 + 1];
1187:           s[2] = b[i2 + 2];
1188:           s[3] = b[i2 + 3];
1189:           s[4] = b[i2 + 4];
1190:           s[5] = b[i2 + 5];
1191:           while (nz--) {
1192:             idx   = 6 * (*vi++);
1193:             xw[0] = x[idx];
1194:             xw[1] = x[1 + idx];
1195:             xw[2] = x[2 + idx];
1196:             xw[3] = x[3 + idx];
1197:             xw[4] = x[4 + idx];
1198:             xw[5] = x[5 + idx];
1199:             PetscKernel_v_gets_v_minus_A_times_w_6(s, v, xw);
1200:             v += 36;
1201:           }
1202:           PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
1203:           x[i2] += xw[0];
1204:           x[i2 + 1] += xw[1];
1205:           x[i2 + 2] += xw[2];
1206:           x[i2 + 3] += xw[3];
1207:           x[i2 + 4] += xw[4];
1208:           x[i2 + 5] += xw[5];
1209:           idiag -= 36;
1210:           i2 -= 6;
1211:         }
1212:         break;
1213:       case 7:
1214:         for (i = m - 1; i >= 0; i--) {
1215:           v    = aa + 49 * ai[i];
1216:           vi   = aj + ai[i];
1217:           nz   = ai[i + 1] - ai[i];
1218:           s[0] = b[i2];
1219:           s[1] = b[i2 + 1];
1220:           s[2] = b[i2 + 2];
1221:           s[3] = b[i2 + 3];
1222:           s[4] = b[i2 + 4];
1223:           s[5] = b[i2 + 5];
1224:           s[6] = b[i2 + 6];
1225:           while (nz--) {
1226:             idx   = 7 * (*vi++);
1227:             xw[0] = x[idx];
1228:             xw[1] = x[1 + idx];
1229:             xw[2] = x[2 + idx];
1230:             xw[3] = x[3 + idx];
1231:             xw[4] = x[4 + idx];
1232:             xw[5] = x[5 + idx];
1233:             xw[6] = x[6 + idx];
1234:             PetscKernel_v_gets_v_minus_A_times_w_7(s, v, xw);
1235:             v += 49;
1236:           }
1237:           PetscKernel_v_gets_A_times_w_7(xw, idiag, s);
1238:           x[i2] += xw[0];
1239:           x[i2 + 1] += xw[1];
1240:           x[i2 + 2] += xw[2];
1241:           x[i2 + 3] += xw[3];
1242:           x[i2 + 4] += xw[4];
1243:           x[i2 + 5] += xw[5];
1244:           x[i2 + 6] += xw[6];
1245:           idiag -= 49;
1246:           i2 -= 7;
1247:         }
1248:         break;
1249:       default:
1250:         for (i = m - 1; i >= 0; i--) {
1251:           v  = aa + bs2 * ai[i];
1252:           vi = aj + ai[i];
1253:           nz = ai[i + 1] - ai[i];

1255:           PetscCall(PetscArraycpy(w, b + i2, bs));
1256:           /* copy all rows of x that are needed into contiguous space */
1257:           workt = work;
1258:           for (j = 0; j < nz; j++) {
1259:             PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
1260:             workt += bs;
1261:           }
1262:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, v, work);
1263:           PetscKernel_w_gets_w_plus_Ar_times_v(bs, bs, w, idiag, x + i2);

1265:           idiag -= bs2;
1266:           i2 -= bs;
1267:         }
1268:         break;
1269:       }
1270:       PetscCall(PetscLogFlops(2.0 * bs2 * (a->nz)));
1271:     }
1272:   }
1273:   PetscCall(VecRestoreArray(xx, &x));
1274:   PetscCall(VecRestoreArrayRead(bb, &b));
1275:   PetscFunctionReturn(PETSC_SUCCESS);
1276: }

1278: /*
1279:     Special version for direct calls from Fortran (Used in PETSc-fun3d)
1280: */
1281: #if defined(PETSC_HAVE_FORTRAN_CAPS)
1282:   #define matsetvaluesblocked4_ MATSETVALUESBLOCKED4
1283: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
1284:   #define matsetvaluesblocked4_ matsetvaluesblocked4
1285: #endif

1287: PETSC_EXTERN void matsetvaluesblocked4_(Mat *AA, PetscInt *mm, const PetscInt im[], PetscInt *nn, const PetscInt in[], const PetscScalar v[])
1288: {
1289:   Mat                A = *AA;
1290:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data;
1291:   PetscInt          *rp, k, low, high, t, ii, jj, row, nrow, i, col, l, N, m = *mm, n = *nn;
1292:   PetscInt          *ai = a->i, *ailen = a->ilen;
1293:   PetscInt          *aj = a->j, stepval, lastcol = -1;
1294:   const PetscScalar *value = v;
1295:   MatScalar         *ap, *aa = a->a, *bap;

1297:   PetscFunctionBegin;
1298:   if (A->rmap->bs != 4) SETERRABORT(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Can only be called with a block size of 4");
1299:   stepval = (n - 1) * 4;
1300:   for (k = 0; k < m; k++) { /* loop over added rows */
1301:     row  = im[k];
1302:     rp   = aj + ai[row];
1303:     ap   = aa + 16 * ai[row];
1304:     nrow = ailen[row];
1305:     low  = 0;
1306:     high = nrow;
1307:     for (l = 0; l < n; l++) { /* loop over added columns */
1308:       col = in[l];
1309:       if (col <= lastcol) low = 0;
1310:       else high = nrow;
1311:       lastcol = col;
1312:       value   = v + k * (stepval + 4 + l) * 4;
1313:       while (high - low > 7) {
1314:         t = (low + high) / 2;
1315:         if (rp[t] > col) high = t;
1316:         else low = t;
1317:       }
1318:       for (i = low; i < high; i++) {
1319:         if (rp[i] > col) break;
1320:         if (rp[i] == col) {
1321:           bap = ap + 16 * i;
1322:           for (ii = 0; ii < 4; ii++, value += stepval) {
1323:             for (jj = ii; jj < 16; jj += 4) bap[jj] += *value++;
1324:           }
1325:           goto noinsert2;
1326:         }
1327:       }
1328:       N = nrow++ - 1;
1329:       high++; /* added new column index thus must search to one higher than before */
1330:       /* shift up all the later entries in this row */
1331:       for (ii = N; ii >= i; ii--) {
1332:         rp[ii + 1] = rp[ii];
1333:         PetscCallVoid(PetscArraycpy(ap + 16 * (ii + 1), ap + 16 * (ii), 16));
1334:       }
1335:       if (N >= i) PetscCallVoid(PetscArrayzero(ap + 16 * i, 16));
1336:       rp[i] = col;
1337:       bap   = ap + 16 * i;
1338:       for (ii = 0; ii < 4; ii++, value += stepval) {
1339:         for (jj = ii; jj < 16; jj += 4) bap[jj] = *value++;
1340:       }
1341:     noinsert2:;
1342:       low = i;
1343:     }
1344:     ailen[row] = nrow;
1345:   }
1346:   PetscFunctionReturnVoid();
1347: }

1349: #if defined(PETSC_HAVE_FORTRAN_CAPS)
1350:   #define matsetvalues4_ MATSETVALUES4
1351: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
1352:   #define matsetvalues4_ matsetvalues4
1353: #endif

1355: PETSC_EXTERN void matsetvalues4_(Mat *AA, PetscInt *mm, PetscInt *im, PetscInt *nn, PetscInt *in, PetscScalar *v)
1356: {
1357:   Mat          A = *AA;
1358:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1359:   PetscInt    *rp, k, low, high, t, row, nrow, i, col, l, N, n = *nn, m = *mm;
1360:   PetscInt    *ai = a->i, *ailen = a->ilen;
1361:   PetscInt    *aj = a->j, brow, bcol;
1362:   PetscInt     ridx, cidx, lastcol = -1;
1363:   MatScalar   *ap, value, *aa      = a->a, *bap;

1365:   PetscFunctionBegin;
1366:   for (k = 0; k < m; k++) { /* loop over added rows */
1367:     row  = im[k];
1368:     brow = row / 4;
1369:     rp   = aj + ai[brow];
1370:     ap   = aa + 16 * ai[brow];
1371:     nrow = ailen[brow];
1372:     low  = 0;
1373:     high = nrow;
1374:     for (l = 0; l < n; l++) { /* loop over added columns */
1375:       col   = in[l];
1376:       bcol  = col / 4;
1377:       ridx  = row % 4;
1378:       cidx  = col % 4;
1379:       value = v[l + k * n];
1380:       if (col <= lastcol) low = 0;
1381:       else high = nrow;
1382:       lastcol = col;
1383:       while (high - low > 7) {
1384:         t = (low + high) / 2;
1385:         if (rp[t] > bcol) high = t;
1386:         else low = t;
1387:       }
1388:       for (i = low; i < high; i++) {
1389:         if (rp[i] > bcol) break;
1390:         if (rp[i] == bcol) {
1391:           bap = ap + 16 * i + 4 * cidx + ridx;
1392:           *bap += value;
1393:           goto noinsert1;
1394:         }
1395:       }
1396:       N = nrow++ - 1;
1397:       high++; /* added new column thus must search to one higher than before */
1398:       /* shift up all the later entries in this row */
1399:       PetscCallVoid(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
1400:       PetscCallVoid(PetscArraymove(ap + 16 * i + 16, ap + 16 * i, 16 * (N - i + 1)));
1401:       PetscCallVoid(PetscArrayzero(ap + 16 * i, 16));
1402:       rp[i]                        = bcol;
1403:       ap[16 * i + 4 * cidx + ridx] = value;
1404:     noinsert1:;
1405:       low = i;
1406:     }
1407:     ailen[brow] = nrow;
1408:   }
1409:   PetscFunctionReturnVoid();
1410: }

1412: /*
1413:      Checks for missing diagonals
1414: */
1415: PetscErrorCode MatMissingDiagonal_SeqBAIJ(Mat A, PetscBool *missing, PetscInt *d)
1416: {
1417:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1418:   PetscInt    *diag, *ii = a->i, i;

1420:   PetscFunctionBegin;
1421:   PetscCall(MatMarkDiagonal_SeqBAIJ(A));
1422:   *missing = PETSC_FALSE;
1423:   if (A->rmap->n > 0 && !ii) {
1424:     *missing = PETSC_TRUE;
1425:     if (d) *d = 0;
1426:     PetscCall(PetscInfo(A, "Matrix has no entries therefore is missing diagonal\n"));
1427:   } else {
1428:     PetscInt n;
1429:     n    = PetscMin(a->mbs, a->nbs);
1430:     diag = a->diag;
1431:     for (i = 0; i < n; i++) {
1432:       if (diag[i] >= ii[i + 1]) {
1433:         *missing = PETSC_TRUE;
1434:         if (d) *d = i;
1435:         PetscCall(PetscInfo(A, "Matrix is missing block diagonal number %" PetscInt_FMT "\n", i));
1436:         break;
1437:       }
1438:     }
1439:   }
1440:   PetscFunctionReturn(PETSC_SUCCESS);
1441: }

1443: PetscErrorCode MatMarkDiagonal_SeqBAIJ(Mat A)
1444: {
1445:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1446:   PetscInt     i, j, m = a->mbs;

1448:   PetscFunctionBegin;
1449:   if (!a->diag) {
1450:     PetscCall(PetscMalloc1(m, &a->diag));
1451:     a->free_diag = PETSC_TRUE;
1452:   }
1453:   for (i = 0; i < m; i++) {
1454:     a->diag[i] = a->i[i + 1];
1455:     for (j = a->i[i]; j < a->i[i + 1]; j++) {
1456:       if (a->j[j] == i) {
1457:         a->diag[i] = j;
1458:         break;
1459:       }
1460:     }
1461:   }
1462:   PetscFunctionReturn(PETSC_SUCCESS);
1463: }

1465: static PetscErrorCode MatGetRowIJ_SeqBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done)
1466: {
1467:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1468:   PetscInt     i, j, n = a->mbs, nz = a->i[n], *tia, *tja, bs = A->rmap->bs, k, l, cnt;
1469:   PetscInt   **ia = (PetscInt **)inia, **ja = (PetscInt **)inja;

1471:   PetscFunctionBegin;
1472:   *nn = n;
1473:   if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
1474:   if (symmetric) {
1475:     PetscCall(MatToSymmetricIJ_SeqAIJ(n, a->i, a->j, PETSC_TRUE, 0, 0, &tia, &tja));
1476:     nz = tia[n];
1477:   } else {
1478:     tia = a->i;
1479:     tja = a->j;
1480:   }

1482:   if (!blockcompressed && bs > 1) {
1483:     (*nn) *= bs;
1484:     /* malloc & create the natural set of indices */
1485:     PetscCall(PetscMalloc1((n + 1) * bs, ia));
1486:     if (n) {
1487:       (*ia)[0] = oshift;
1488:       for (j = 1; j < bs; j++) (*ia)[j] = (tia[1] - tia[0]) * bs + (*ia)[j - 1];
1489:     }

1491:     for (i = 1; i < n; i++) {
1492:       (*ia)[i * bs] = (tia[i] - tia[i - 1]) * bs + (*ia)[i * bs - 1];
1493:       for (j = 1; j < bs; j++) (*ia)[i * bs + j] = (tia[i + 1] - tia[i]) * bs + (*ia)[i * bs + j - 1];
1494:     }
1495:     if (n) (*ia)[n * bs] = (tia[n] - tia[n - 1]) * bs + (*ia)[n * bs - 1];

1497:     if (inja) {
1498:       PetscCall(PetscMalloc1(nz * bs * bs, ja));
1499:       cnt = 0;
1500:       for (i = 0; i < n; i++) {
1501:         for (j = 0; j < bs; j++) {
1502:           for (k = tia[i]; k < tia[i + 1]; k++) {
1503:             for (l = 0; l < bs; l++) (*ja)[cnt++] = bs * tja[k] + l;
1504:           }
1505:         }
1506:       }
1507:     }

1509:     if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
1510:       PetscCall(PetscFree(tia));
1511:       PetscCall(PetscFree(tja));
1512:     }
1513:   } else if (oshift == 1) {
1514:     if (symmetric) {
1515:       nz = tia[A->rmap->n / bs];
1516:       /*  add 1 to i and j indices */
1517:       for (i = 0; i < A->rmap->n / bs + 1; i++) tia[i] = tia[i] + 1;
1518:       *ia = tia;
1519:       if (ja) {
1520:         for (i = 0; i < nz; i++) tja[i] = tja[i] + 1;
1521:         *ja = tja;
1522:       }
1523:     } else {
1524:       nz = a->i[A->rmap->n / bs];
1525:       /* malloc space and  add 1 to i and j indices */
1526:       PetscCall(PetscMalloc1(A->rmap->n / bs + 1, ia));
1527:       for (i = 0; i < A->rmap->n / bs + 1; i++) (*ia)[i] = a->i[i] + 1;
1528:       if (ja) {
1529:         PetscCall(PetscMalloc1(nz, ja));
1530:         for (i = 0; i < nz; i++) (*ja)[i] = a->j[i] + 1;
1531:       }
1532:     }
1533:   } else {
1534:     *ia = tia;
1535:     if (ja) *ja = tja;
1536:   }
1537:   PetscFunctionReturn(PETSC_SUCCESS);
1538: }

1540: static PetscErrorCode MatRestoreRowIJ_SeqBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done)
1541: {
1542:   PetscFunctionBegin;
1543:   if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
1544:   if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
1545:     PetscCall(PetscFree(*ia));
1546:     if (ja) PetscCall(PetscFree(*ja));
1547:   }
1548:   PetscFunctionReturn(PETSC_SUCCESS);
1549: }

1551: PetscErrorCode MatDestroy_SeqBAIJ(Mat A)
1552: {
1553:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;

1555:   PetscFunctionBegin;
1556:   if (A->hash_active) {
1557:     PetscInt bs;
1558:     A->ops[0] = a->cops;
1559:     PetscCall(PetscHMapIJVDestroy(&a->ht));
1560:     PetscCall(MatGetBlockSize(A, &bs));
1561:     if (bs > 1) PetscCall(PetscHSetIJDestroy(&a->bht));
1562:     PetscCall(PetscFree(a->dnz));
1563:     PetscCall(PetscFree(a->bdnz));
1564:     A->hash_active = PETSC_FALSE;
1565:   }
1566:   PetscCall(PetscLogObjectState((PetscObject)A, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, A->rmap->N, A->cmap->n, a->nz));
1567:   PetscCall(MatSeqXAIJFreeAIJ(A, &a->a, &a->j, &a->i));
1568:   PetscCall(ISDestroy(&a->row));
1569:   PetscCall(ISDestroy(&a->col));
1570:   if (a->free_diag) PetscCall(PetscFree(a->diag));
1571:   PetscCall(PetscFree(a->idiag));
1572:   if (a->free_imax_ilen) PetscCall(PetscFree2(a->imax, a->ilen));
1573:   PetscCall(PetscFree(a->solve_work));
1574:   PetscCall(PetscFree(a->mult_work));
1575:   PetscCall(PetscFree(a->sor_workt));
1576:   PetscCall(PetscFree(a->sor_work));
1577:   PetscCall(ISDestroy(&a->icol));
1578:   PetscCall(PetscFree(a->saved_values));
1579:   PetscCall(PetscFree2(a->compressedrow.i, a->compressedrow.rindex));

1581:   PetscCall(MatDestroy(&a->sbaijMat));
1582:   PetscCall(MatDestroy(&a->parent));
1583:   PetscCall(PetscFree(A->data));

1585:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
1586:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJGetArray_C", NULL));
1587:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJRestoreArray_C", NULL));
1588:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatStoreValues_C", NULL));
1589:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatRetrieveValues_C", NULL));
1590:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJSetColumnIndices_C", NULL));
1591:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_seqaij_C", NULL));
1592:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_seqsbaij_C", NULL));
1593:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJSetPreallocation_C", NULL));
1594:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJSetPreallocationCSR_C", NULL));
1595:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_seqbstrm_C", NULL));
1596:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatIsTranspose_C", NULL));
1597: #if defined(PETSC_HAVE_HYPRE)
1598:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_hypre_C", NULL));
1599: #endif
1600:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_is_C", NULL));
1601:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
1602:   PetscFunctionReturn(PETSC_SUCCESS);
1603: }

1605: static PetscErrorCode MatSetOption_SeqBAIJ(Mat A, MatOption op, PetscBool flg)
1606: {
1607:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;

1609:   PetscFunctionBegin;
1610:   switch (op) {
1611:   case MAT_ROW_ORIENTED:
1612:     a->roworiented = flg;
1613:     break;
1614:   case MAT_KEEP_NONZERO_PATTERN:
1615:     a->keepnonzeropattern = flg;
1616:     break;
1617:   case MAT_NEW_NONZERO_LOCATIONS:
1618:     a->nonew = (flg ? 0 : 1);
1619:     break;
1620:   case MAT_NEW_NONZERO_LOCATION_ERR:
1621:     a->nonew = (flg ? -1 : 0);
1622:     break;
1623:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1624:     a->nonew = (flg ? -2 : 0);
1625:     break;
1626:   case MAT_UNUSED_NONZERO_LOCATION_ERR:
1627:     a->nounused = (flg ? -1 : 0);
1628:     break;
1629:   case MAT_FORCE_DIAGONAL_ENTRIES:
1630:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1631:   case MAT_USE_HASH_TABLE:
1632:   case MAT_SORTED_FULL:
1633:     PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op]));
1634:     break;
1635:   case MAT_SPD:
1636:   case MAT_SYMMETRIC:
1637:   case MAT_STRUCTURALLY_SYMMETRIC:
1638:   case MAT_HERMITIAN:
1639:   case MAT_SYMMETRY_ETERNAL:
1640:   case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
1641:   case MAT_SUBMAT_SINGLEIS:
1642:   case MAT_STRUCTURE_ONLY:
1643:   case MAT_SPD_ETERNAL:
1644:     /* if the diagonal matrix is square it inherits some of the properties above */
1645:     break;
1646:   default:
1647:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %d", op);
1648:   }
1649:   PetscFunctionReturn(PETSC_SUCCESS);
1650: }

1652: /* used for both SeqBAIJ and SeqSBAIJ matrices */
1653: PetscErrorCode MatGetRow_SeqBAIJ_private(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v, PetscInt *ai, PetscInt *aj, PetscScalar *aa)
1654: {
1655:   PetscInt     itmp, i, j, k, M, bn, bp, *idx_i, bs, bs2;
1656:   MatScalar   *aa_i;
1657:   PetscScalar *v_i;

1659:   PetscFunctionBegin;
1660:   bs  = A->rmap->bs;
1661:   bs2 = bs * bs;
1662:   PetscCheck(row >= 0 && row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " out of range", row);

1664:   bn  = row / bs; /* Block number */
1665:   bp  = row % bs; /* Block Position */
1666:   M   = ai[bn + 1] - ai[bn];
1667:   *nz = bs * M;

1669:   if (v) {
1670:     *v = NULL;
1671:     if (*nz) {
1672:       PetscCall(PetscMalloc1(*nz, v));
1673:       for (i = 0; i < M; i++) { /* for each block in the block row */
1674:         v_i  = *v + i * bs;
1675:         aa_i = aa + bs2 * (ai[bn] + i);
1676:         for (j = bp, k = 0; j < bs2; j += bs, k++) v_i[k] = aa_i[j];
1677:       }
1678:     }
1679:   }

1681:   if (idx) {
1682:     *idx = NULL;
1683:     if (*nz) {
1684:       PetscCall(PetscMalloc1(*nz, idx));
1685:       for (i = 0; i < M; i++) { /* for each block in the block row */
1686:         idx_i = *idx + i * bs;
1687:         itmp  = bs * aj[ai[bn] + i];
1688:         for (j = 0; j < bs; j++) idx_i[j] = itmp++;
1689:       }
1690:     }
1691:   }
1692:   PetscFunctionReturn(PETSC_SUCCESS);
1693: }

1695: PetscErrorCode MatGetRow_SeqBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1696: {
1697:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;

1699:   PetscFunctionBegin;
1700:   PetscCall(MatGetRow_SeqBAIJ_private(A, row, nz, idx, v, a->i, a->j, a->a));
1701:   PetscFunctionReturn(PETSC_SUCCESS);
1702: }

1704: PetscErrorCode MatRestoreRow_SeqBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1705: {
1706:   PetscFunctionBegin;
1707:   if (idx) PetscCall(PetscFree(*idx));
1708:   if (v) PetscCall(PetscFree(*v));
1709:   PetscFunctionReturn(PETSC_SUCCESS);
1710: }

1712: static PetscErrorCode MatTranspose_SeqBAIJ(Mat A, MatReuse reuse, Mat *B)
1713: {
1714:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *at;
1715:   Mat          C;
1716:   PetscInt     i, j, k, *aj = a->j, *ai = a->i, bs = A->rmap->bs, mbs = a->mbs, nbs = a->nbs, *atfill;
1717:   PetscInt     bs2 = a->bs2, *ati, *atj, anzj, kr;
1718:   MatScalar   *ata, *aa = a->a;

1720:   PetscFunctionBegin;
1721:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1722:   PetscCall(PetscCalloc1(1 + nbs, &atfill));
1723:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1724:     for (i = 0; i < ai[mbs]; i++) atfill[aj[i]] += 1; /* count num of non-zeros in row aj[i] */

1726:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1727:     PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->N, A->cmap->n, A->rmap->N));
1728:     PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
1729:     PetscCall(MatSeqBAIJSetPreallocation(C, bs, 0, atfill));

1731:     at  = (Mat_SeqBAIJ *)C->data;
1732:     ati = at->i;
1733:     for (i = 0; i < nbs; i++) at->ilen[i] = at->imax[i] = ati[i + 1] - ati[i];
1734:   } else {
1735:     C   = *B;
1736:     at  = (Mat_SeqBAIJ *)C->data;
1737:     ati = at->i;
1738:   }

1740:   atj = at->j;
1741:   ata = at->a;

1743:   /* Copy ati into atfill so we have locations of the next free space in atj */
1744:   PetscCall(PetscArraycpy(atfill, ati, nbs));

1746:   /* Walk through A row-wise and mark nonzero entries of A^T. */
1747:   for (i = 0; i < mbs; i++) {
1748:     anzj = ai[i + 1] - ai[i];
1749:     for (j = 0; j < anzj; j++) {
1750:       atj[atfill[*aj]] = i;
1751:       for (kr = 0; kr < bs; kr++) {
1752:         for (k = 0; k < bs; k++) ata[bs2 * atfill[*aj] + k * bs + kr] = *aa++;
1753:       }
1754:       atfill[*aj++] += 1;
1755:     }
1756:   }
1757:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1758:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

1760:   /* Clean up temporary space and complete requests. */
1761:   PetscCall(PetscFree(atfill));

1763:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1764:     PetscCall(MatSetBlockSizes(C, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
1765:     *B = C;
1766:   } else {
1767:     PetscCall(MatHeaderMerge(A, &C));
1768:   }
1769:   PetscFunctionReturn(PETSC_SUCCESS);
1770: }

1772: static PetscErrorCode MatIsTranspose_SeqBAIJ(Mat A, Mat B, PetscReal tol, PetscBool *f)
1773: {
1774:   Mat Btrans;

1776:   PetscFunctionBegin;
1777:   *f = PETSC_FALSE;
1778:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Btrans));
1779:   PetscCall(MatEqual_SeqBAIJ(B, Btrans, f));
1780:   PetscCall(MatDestroy(&Btrans));
1781:   PetscFunctionReturn(PETSC_SUCCESS);
1782: }

1784: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
1785: PetscErrorCode MatView_SeqBAIJ_Binary(Mat mat, PetscViewer viewer)
1786: {
1787:   Mat_SeqBAIJ *A = (Mat_SeqBAIJ *)mat->data;
1788:   PetscInt     header[4], M, N, m, bs, nz, cnt, i, j, k, l;
1789:   PetscInt    *rowlens, *colidxs;
1790:   PetscScalar *matvals;

1792:   PetscFunctionBegin;
1793:   PetscCall(PetscViewerSetUp(viewer));

1795:   M  = mat->rmap->N;
1796:   N  = mat->cmap->N;
1797:   m  = mat->rmap->n;
1798:   bs = mat->rmap->bs;
1799:   nz = bs * bs * A->nz;

1801:   /* write matrix header */
1802:   header[0] = MAT_FILE_CLASSID;
1803:   header[1] = M;
1804:   header[2] = N;
1805:   header[3] = nz;
1806:   PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));

1808:   /* store row lengths */
1809:   PetscCall(PetscMalloc1(m, &rowlens));
1810:   for (cnt = 0, i = 0; i < A->mbs; i++)
1811:     for (j = 0; j < bs; j++) rowlens[cnt++] = bs * (A->i[i + 1] - A->i[i]);
1812:   PetscCall(PetscViewerBinaryWrite(viewer, rowlens, m, PETSC_INT));
1813:   PetscCall(PetscFree(rowlens));

1815:   /* store column indices  */
1816:   PetscCall(PetscMalloc1(nz, &colidxs));
1817:   for (cnt = 0, i = 0; i < A->mbs; i++)
1818:     for (k = 0; k < bs; k++)
1819:       for (j = A->i[i]; j < A->i[i + 1]; j++)
1820:         for (l = 0; l < bs; l++) colidxs[cnt++] = bs * A->j[j] + l;
1821:   PetscCheck(cnt == nz, PETSC_COMM_SELF, PETSC_ERR_LIB, "Internal PETSc error: cnt = %" PetscInt_FMT " nz = %" PetscInt_FMT, cnt, nz);
1822:   PetscCall(PetscViewerBinaryWrite(viewer, colidxs, nz, PETSC_INT));
1823:   PetscCall(PetscFree(colidxs));

1825:   /* store nonzero values */
1826:   PetscCall(PetscMalloc1(nz, &matvals));
1827:   for (cnt = 0, i = 0; i < A->mbs; i++)
1828:     for (k = 0; k < bs; k++)
1829:       for (j = A->i[i]; j < A->i[i + 1]; j++)
1830:         for (l = 0; l < bs; l++) matvals[cnt++] = A->a[bs * (bs * j + l) + k];
1831:   PetscCheck(cnt == nz, PETSC_COMM_SELF, PETSC_ERR_LIB, "Internal PETSc error: cnt = %" PetscInt_FMT " nz = %" PetscInt_FMT, cnt, nz);
1832:   PetscCall(PetscViewerBinaryWrite(viewer, matvals, nz, PETSC_SCALAR));
1833:   PetscCall(PetscFree(matvals));

1835:   /* write block size option to the viewer's .info file */
1836:   PetscCall(MatView_Binary_BlockSizes(mat, viewer));
1837:   PetscFunctionReturn(PETSC_SUCCESS);
1838: }

1840: static PetscErrorCode MatView_SeqBAIJ_ASCII_structonly(Mat A, PetscViewer viewer)
1841: {
1842:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1843:   PetscInt     i, bs = A->rmap->bs, k;

1845:   PetscFunctionBegin;
1846:   PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1847:   for (i = 0; i < a->mbs; i++) {
1848:     PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT "-%" PetscInt_FMT ":", i * bs, i * bs + bs - 1));
1849:     for (k = a->i[i]; k < a->i[i + 1]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT "-%" PetscInt_FMT ") ", bs * a->j[k], bs * a->j[k] + bs - 1));
1850:     PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1851:   }
1852:   PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1853:   PetscFunctionReturn(PETSC_SUCCESS);
1854: }

1856: static PetscErrorCode MatView_SeqBAIJ_ASCII(Mat A, PetscViewer viewer)
1857: {
1858:   Mat_SeqBAIJ      *a = (Mat_SeqBAIJ *)A->data;
1859:   PetscInt          i, j, bs = A->rmap->bs, k, l, bs2 = a->bs2;
1860:   PetscViewerFormat format;

1862:   PetscFunctionBegin;
1863:   if (A->structure_only) {
1864:     PetscCall(MatView_SeqBAIJ_ASCII_structonly(A, viewer));
1865:     PetscFunctionReturn(PETSC_SUCCESS);
1866:   }

1868:   PetscCall(PetscViewerGetFormat(viewer, &format));
1869:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1870:     PetscCall(PetscViewerASCIIPrintf(viewer, "  block size is %" PetscInt_FMT "\n", bs));
1871:   } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
1872:     const char *matname;
1873:     Mat         aij;
1874:     PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &aij));
1875:     PetscCall(PetscObjectGetName((PetscObject)A, &matname));
1876:     PetscCall(PetscObjectSetName((PetscObject)aij, matname));
1877:     PetscCall(MatView(aij, viewer));
1878:     PetscCall(MatDestroy(&aij));
1879:   } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
1880:     PetscFunctionReturn(PETSC_SUCCESS);
1881:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1882:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1883:     for (i = 0; i < a->mbs; i++) {
1884:       for (j = 0; j < bs; j++) {
1885:         PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j));
1886:         for (k = a->i[i]; k < a->i[i + 1]; k++) {
1887:           for (l = 0; l < bs; l++) {
1888: #if defined(PETSC_USE_COMPLEX)
1889:             if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0 && PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
1890:               PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %gi) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), (double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j])));
1891:             } else if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0 && PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
1892:               PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %gi) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), -(double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j])));
1893:             } else if (PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
1894:               PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j])));
1895:             }
1896: #else
1897:             if (a->a[bs2 * k + l * bs + j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)a->a[bs2 * k + l * bs + j]));
1898: #endif
1899:           }
1900:         }
1901:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1902:       }
1903:     }
1904:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1905:   } else {
1906:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1907:     for (i = 0; i < a->mbs; i++) {
1908:       for (j = 0; j < bs; j++) {
1909:         PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j));
1910:         for (k = a->i[i]; k < a->i[i + 1]; k++) {
1911:           for (l = 0; l < bs; l++) {
1912: #if defined(PETSC_USE_COMPLEX)
1913:             if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0) {
1914:               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])));
1915:             } else if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0) {
1916:               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])));
1917:             } else {
1918:               PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j])));
1919:             }
1920: #else
1921:             PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)a->a[bs2 * k + l * bs + j]));
1922: #endif
1923:           }
1924:         }
1925:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1926:       }
1927:     }
1928:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1929:   }
1930:   PetscCall(PetscViewerFlush(viewer));
1931:   PetscFunctionReturn(PETSC_SUCCESS);
1932: }

1934: #include <petscdraw.h>
1935: static PetscErrorCode MatView_SeqBAIJ_Draw_Zoom(PetscDraw draw, void *Aa)
1936: {
1937:   Mat               A = (Mat)Aa;
1938:   Mat_SeqBAIJ      *a = (Mat_SeqBAIJ *)A->data;
1939:   PetscInt          row, i, j, k, l, mbs = a->mbs, color, bs = A->rmap->bs, bs2 = a->bs2;
1940:   PetscReal         xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1941:   MatScalar        *aa;
1942:   PetscViewer       viewer;
1943:   PetscViewerFormat format;

1945:   PetscFunctionBegin;
1946:   PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1947:   PetscCall(PetscViewerGetFormat(viewer, &format));
1948:   PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));

1950:   /* loop over matrix elements drawing boxes */

1952:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1953:     PetscDrawCollectiveBegin(draw);
1954:     /* Blue for negative, Cyan for zero and  Red for positive */
1955:     color = PETSC_DRAW_BLUE;
1956:     for (i = 0, row = 0; i < mbs; i++, row += bs) {
1957:       for (j = a->i[i]; j < a->i[i + 1]; j++) {
1958:         y_l = A->rmap->N - row - 1.0;
1959:         y_r = y_l + 1.0;
1960:         x_l = a->j[j] * bs;
1961:         x_r = x_l + 1.0;
1962:         aa  = a->a + j * bs2;
1963:         for (k = 0; k < bs; k++) {
1964:           for (l = 0; l < bs; l++) {
1965:             if (PetscRealPart(*aa++) >= 0.) continue;
1966:             PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
1967:           }
1968:         }
1969:       }
1970:     }
1971:     color = PETSC_DRAW_CYAN;
1972:     for (i = 0, row = 0; i < mbs; i++, row += bs) {
1973:       for (j = a->i[i]; j < a->i[i + 1]; j++) {
1974:         y_l = A->rmap->N - row - 1.0;
1975:         y_r = y_l + 1.0;
1976:         x_l = a->j[j] * bs;
1977:         x_r = x_l + 1.0;
1978:         aa  = a->a + j * bs2;
1979:         for (k = 0; k < bs; k++) {
1980:           for (l = 0; l < bs; l++) {
1981:             if (PetscRealPart(*aa++) != 0.) continue;
1982:             PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
1983:           }
1984:         }
1985:       }
1986:     }
1987:     color = PETSC_DRAW_RED;
1988:     for (i = 0, row = 0; i < mbs; i++, row += bs) {
1989:       for (j = a->i[i]; j < a->i[i + 1]; j++) {
1990:         y_l = A->rmap->N - row - 1.0;
1991:         y_r = y_l + 1.0;
1992:         x_l = a->j[j] * bs;
1993:         x_r = x_l + 1.0;
1994:         aa  = a->a + j * bs2;
1995:         for (k = 0; k < bs; k++) {
1996:           for (l = 0; l < bs; l++) {
1997:             if (PetscRealPart(*aa++) <= 0.) continue;
1998:             PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
1999:           }
2000:         }
2001:       }
2002:     }
2003:     PetscDrawCollectiveEnd(draw);
2004:   } else {
2005:     /* use contour shading to indicate magnitude of values */
2006:     /* first determine max of all nonzero values */
2007:     PetscReal minv = 0.0, maxv = 0.0;
2008:     PetscDraw popup;

2010:     for (i = 0; i < a->nz * a->bs2; i++) {
2011:       if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]);
2012:     }
2013:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
2014:     PetscCall(PetscDrawGetPopup(draw, &popup));
2015:     PetscCall(PetscDrawScalePopup(popup, 0.0, maxv));

2017:     PetscDrawCollectiveBegin(draw);
2018:     for (i = 0, row = 0; i < mbs; i++, row += bs) {
2019:       for (j = a->i[i]; j < a->i[i + 1]; j++) {
2020:         y_l = A->rmap->N - row - 1.0;
2021:         y_r = y_l + 1.0;
2022:         x_l = a->j[j] * bs;
2023:         x_r = x_l + 1.0;
2024:         aa  = a->a + j * bs2;
2025:         for (k = 0; k < bs; k++) {
2026:           for (l = 0; l < bs; l++) {
2027:             MatScalar v = *aa++;
2028:             color       = PetscDrawRealToColor(PetscAbsScalar(v), minv, maxv);
2029:             PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
2030:           }
2031:         }
2032:       }
2033:     }
2034:     PetscDrawCollectiveEnd(draw);
2035:   }
2036:   PetscFunctionReturn(PETSC_SUCCESS);
2037: }

2039: static PetscErrorCode MatView_SeqBAIJ_Draw(Mat A, PetscViewer viewer)
2040: {
2041:   PetscReal xl, yl, xr, yr, w, h;
2042:   PetscDraw draw;
2043:   PetscBool isnull;

2045:   PetscFunctionBegin;
2046:   PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2047:   PetscCall(PetscDrawIsNull(draw, &isnull));
2048:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

2050:   xr = A->cmap->n;
2051:   yr = A->rmap->N;
2052:   h  = yr / 10.0;
2053:   w  = xr / 10.0;
2054:   xr += w;
2055:   yr += h;
2056:   xl = -w;
2057:   yl = -h;
2058:   PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
2059:   PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
2060:   PetscCall(PetscDrawZoom(draw, MatView_SeqBAIJ_Draw_Zoom, A));
2061:   PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
2062:   PetscCall(PetscDrawSave(draw));
2063:   PetscFunctionReturn(PETSC_SUCCESS);
2064: }

2066: PetscErrorCode MatView_SeqBAIJ(Mat A, PetscViewer viewer)
2067: {
2068:   PetscBool iascii, isbinary, isdraw;

2070:   PetscFunctionBegin;
2071:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2072:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2073:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2074:   if (iascii) {
2075:     PetscCall(MatView_SeqBAIJ_ASCII(A, viewer));
2076:   } else if (isbinary) {
2077:     PetscCall(MatView_SeqBAIJ_Binary(A, viewer));
2078:   } else if (isdraw) {
2079:     PetscCall(MatView_SeqBAIJ_Draw(A, viewer));
2080:   } else {
2081:     Mat B;
2082:     PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &B));
2083:     PetscCall(MatView(B, viewer));
2084:     PetscCall(MatDestroy(&B));
2085:   }
2086:   PetscFunctionReturn(PETSC_SUCCESS);
2087: }

2089: PetscErrorCode MatGetValues_SeqBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], PetscScalar v[])
2090: {
2091:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2092:   PetscInt    *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j;
2093:   PetscInt    *ai = a->i, *ailen = a->ilen;
2094:   PetscInt     brow, bcol, ridx, cidx, bs = A->rmap->bs, bs2 = a->bs2;
2095:   MatScalar   *ap, *aa = a->a;

2097:   PetscFunctionBegin;
2098:   for (k = 0; k < m; k++) { /* loop over rows */
2099:     row  = im[k];
2100:     brow = row / bs;
2101:     if (row < 0) {
2102:       v += n;
2103:       continue;
2104:     } /* negative row */
2105:     PetscCheck(row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " too large", row);
2106:     rp   = PetscSafePointerPlusOffset(aj, ai[brow]);
2107:     ap   = PetscSafePointerPlusOffset(aa, bs2 * ai[brow]);
2108:     nrow = ailen[brow];
2109:     for (l = 0; l < n; l++) { /* loop over columns */
2110:       if (in[l] < 0) {
2111:         v++;
2112:         continue;
2113:       } /* negative column */
2114:       PetscCheck(in[l] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column %" PetscInt_FMT " too large", in[l]);
2115:       col  = in[l];
2116:       bcol = col / bs;
2117:       cidx = col % bs;
2118:       ridx = row % bs;
2119:       high = nrow;
2120:       low  = 0; /* assume unsorted */
2121:       while (high - low > 5) {
2122:         t = (low + high) / 2;
2123:         if (rp[t] > bcol) high = t;
2124:         else low = t;
2125:       }
2126:       for (i = low; i < high; i++) {
2127:         if (rp[i] > bcol) break;
2128:         if (rp[i] == bcol) {
2129:           *v++ = ap[bs2 * i + bs * cidx + ridx];
2130:           goto finished;
2131:         }
2132:       }
2133:       *v++ = 0.0;
2134:     finished:;
2135:     }
2136:   }
2137:   PetscFunctionReturn(PETSC_SUCCESS);
2138: }

2140: PetscErrorCode MatSetValuesBlocked_SeqBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
2141: {
2142:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data;
2143:   PetscInt          *rp, k, low, high, t, ii, jj, row, nrow, i, col, l, rmax, N, lastcol = -1;
2144:   PetscInt          *imax = a->imax, *ai = a->i, *ailen = a->ilen;
2145:   PetscInt          *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs, stepval;
2146:   PetscBool          roworiented = a->roworiented;
2147:   const PetscScalar *value       = v;
2148:   MatScalar         *ap = NULL, *aa = a->a, *bap;

2150:   PetscFunctionBegin;
2151:   if (roworiented) {
2152:     stepval = (n - 1) * bs;
2153:   } else {
2154:     stepval = (m - 1) * bs;
2155:   }
2156:   for (k = 0; k < m; k++) { /* loop over added rows */
2157:     row = im[k];
2158:     if (row < 0) continue;
2159:     PetscCheck(row < a->mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block row index too large %" PetscInt_FMT " max %" PetscInt_FMT, row, a->mbs - 1);
2160:     rp = aj + ai[row];
2161:     if (!A->structure_only) ap = aa + bs2 * ai[row];
2162:     rmax = imax[row];
2163:     nrow = ailen[row];
2164:     low  = 0;
2165:     high = nrow;
2166:     for (l = 0; l < n; l++) { /* loop over added columns */
2167:       if (in[l] < 0) continue;
2168:       PetscCheck(in[l] < a->nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block column index too large %" PetscInt_FMT " max %" PetscInt_FMT, in[l], a->nbs - 1);
2169:       col = in[l];
2170:       if (!A->structure_only) {
2171:         if (roworiented) {
2172:           value = v + (k * (stepval + bs) + l) * bs;
2173:         } else {
2174:           value = v + (l * (stepval + bs) + k) * bs;
2175:         }
2176:       }
2177:       if (col <= lastcol) low = 0;
2178:       else high = nrow;
2179:       lastcol = col;
2180:       while (high - low > 7) {
2181:         t = (low + high) / 2;
2182:         if (rp[t] > col) high = t;
2183:         else low = t;
2184:       }
2185:       for (i = low; i < high; i++) {
2186:         if (rp[i] > col) break;
2187:         if (rp[i] == col) {
2188:           if (A->structure_only) goto noinsert2;
2189:           bap = ap + bs2 * i;
2190:           if (roworiented) {
2191:             if (is == ADD_VALUES) {
2192:               for (ii = 0; ii < bs; ii++, value += stepval) {
2193:                 for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
2194:               }
2195:             } else {
2196:               for (ii = 0; ii < bs; ii++, value += stepval) {
2197:                 for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
2198:               }
2199:             }
2200:           } else {
2201:             if (is == ADD_VALUES) {
2202:               for (ii = 0; ii < bs; ii++, value += bs + stepval) {
2203:                 for (jj = 0; jj < bs; jj++) bap[jj] += value[jj];
2204:                 bap += bs;
2205:               }
2206:             } else {
2207:               for (ii = 0; ii < bs; ii++, value += bs + stepval) {
2208:                 for (jj = 0; jj < bs; jj++) bap[jj] = value[jj];
2209:                 bap += bs;
2210:               }
2211:             }
2212:           }
2213:           goto noinsert2;
2214:         }
2215:       }
2216:       if (nonew == 1) goto noinsert2;
2217:       PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new blocked index new nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col);
2218:       if (A->structure_only) {
2219:         MatSeqXAIJReallocateAIJ_structure_only(A, a->mbs, bs2, nrow, row, col, rmax, ai, aj, rp, imax, nonew, MatScalar);
2220:       } else {
2221:         MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
2222:       }
2223:       N = nrow++ - 1;
2224:       high++;
2225:       /* shift up all the later entries in this row */
2226:       PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
2227:       rp[i] = col;
2228:       if (!A->structure_only) {
2229:         PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
2230:         bap = ap + bs2 * i;
2231:         if (roworiented) {
2232:           for (ii = 0; ii < bs; ii++, value += stepval) {
2233:             for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
2234:           }
2235:         } else {
2236:           for (ii = 0; ii < bs; ii++, value += stepval) {
2237:             for (jj = 0; jj < bs; jj++) *bap++ = *value++;
2238:           }
2239:         }
2240:       }
2241:     noinsert2:;
2242:       low = i;
2243:     }
2244:     ailen[row] = nrow;
2245:   }
2246:   PetscFunctionReturn(PETSC_SUCCESS);
2247: }

2249: PetscErrorCode MatAssemblyEnd_SeqBAIJ(Mat A, MatAssemblyType mode)
2250: {
2251:   Mat_SeqBAIJ *a      = (Mat_SeqBAIJ *)A->data;
2252:   PetscInt     fshift = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax;
2253:   PetscInt     m = A->rmap->N, *ip, N, *ailen = a->ilen;
2254:   PetscInt     mbs = a->mbs, bs2 = a->bs2, rmax = 0;
2255:   MatScalar   *aa    = a->a, *ap;
2256:   PetscReal    ratio = 0.6;

2258:   PetscFunctionBegin;
2259:   if (mode == MAT_FLUSH_ASSEMBLY || (A->was_assembled && A->ass_nonzerostate == A->nonzerostate)) PetscFunctionReturn(PETSC_SUCCESS);

2261:   if (m) rmax = ailen[0];
2262:   for (i = 1; i < mbs; i++) {
2263:     /* move each row back by the amount of empty slots (fshift) before it*/
2264:     fshift += imax[i - 1] - ailen[i - 1];
2265:     rmax = PetscMax(rmax, ailen[i]);
2266:     if (fshift) {
2267:       ip = aj + ai[i];
2268:       ap = aa + bs2 * ai[i];
2269:       N  = ailen[i];
2270:       PetscCall(PetscArraymove(ip - fshift, ip, N));
2271:       if (!A->structure_only) PetscCall(PetscArraymove(ap - bs2 * fshift, ap, bs2 * N));
2272:     }
2273:     ai[i] = ai[i - 1] + ailen[i - 1];
2274:   }
2275:   if (mbs) {
2276:     fshift += imax[mbs - 1] - ailen[mbs - 1];
2277:     ai[mbs] = ai[mbs - 1] + ailen[mbs - 1];
2278:   }

2280:   /* reset ilen and imax for each row */
2281:   a->nonzerorowcnt = 0;
2282:   if (A->structure_only) {
2283:     PetscCall(PetscFree2(a->imax, a->ilen));
2284:   } else { /* !A->structure_only */
2285:     for (i = 0; i < mbs; i++) {
2286:       ailen[i] = imax[i] = ai[i + 1] - ai[i];
2287:       a->nonzerorowcnt += ((ai[i + 1] - ai[i]) > 0);
2288:     }
2289:   }
2290:   a->nz = ai[mbs];

2292:   /* diagonals may have moved, so kill the diagonal pointers */
2293:   a->idiagvalid = PETSC_FALSE;
2294:   if (fshift && a->diag) PetscCall(PetscFree(a->diag));
2295:   if (fshift) PetscCheck(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);
2296:   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->cmap->n, A->rmap->bs, fshift * bs2, a->nz * bs2));
2297:   PetscCall(PetscInfo(A, "Number of mallocs during MatSetValues is %" PetscInt_FMT "\n", a->reallocs));
2298:   PetscCall(PetscInfo(A, "Most nonzeros blocks in any row is %" PetscInt_FMT "\n", rmax));

2300:   A->info.mallocs += a->reallocs;
2301:   a->reallocs         = 0;
2302:   A->info.nz_unneeded = (PetscReal)fshift * bs2;
2303:   a->rmax             = rmax;

2305:   if (!A->structure_only) PetscCall(MatCheckCompressedRow(A, a->nonzerorowcnt, &a->compressedrow, a->i, mbs, ratio));
2306:   PetscFunctionReturn(PETSC_SUCCESS);
2307: }

2309: /*
2310:    This function returns an array of flags which indicate the locations of contiguous
2311:    blocks that should be zeroed. for eg: if bs = 3  and is = [0,1,2,3,5,6,7,8,9]
2312:    then the resulting sizes = [3,1,1,3,1] corresponding to sets [(0,1,2),(3),(5),(6,7,8),(9)]
2313:    Assume: sizes should be long enough to hold all the values.
2314: */
2315: static PetscErrorCode MatZeroRows_SeqBAIJ_Check_Blocks(PetscInt idx[], PetscInt n, PetscInt bs, PetscInt sizes[], PetscInt *bs_max)
2316: {
2317:   PetscInt j = 0;

2319:   PetscFunctionBegin;
2320:   for (PetscInt i = 0; i < n; j++) {
2321:     PetscInt row = idx[i];
2322:     if (row % bs != 0) { /* Not the beginning of a block */
2323:       sizes[j] = 1;
2324:       i++;
2325:     } else if (i + bs > n) { /* complete block doesn't exist (at idx end) */
2326:       sizes[j] = 1;          /* Also makes sure at least 'bs' values exist for next else */
2327:       i++;
2328:     } else { /* Beginning of the block, so check if the complete block exists */
2329:       PetscBool flg = PETSC_TRUE;
2330:       for (PetscInt k = 1; k < bs; k++) {
2331:         if (row + k != idx[i + k]) { /* break in the block */
2332:           flg = PETSC_FALSE;
2333:           break;
2334:         }
2335:       }
2336:       if (flg) { /* No break in the bs */
2337:         sizes[j] = bs;
2338:         i += bs;
2339:       } else {
2340:         sizes[j] = 1;
2341:         i++;
2342:       }
2343:     }
2344:   }
2345:   *bs_max = j;
2346:   PetscFunctionReturn(PETSC_SUCCESS);
2347: }

2349: PetscErrorCode MatZeroRows_SeqBAIJ(Mat A, PetscInt is_n, const PetscInt is_idx[], PetscScalar diag, Vec x, Vec b)
2350: {
2351:   Mat_SeqBAIJ       *baij = (Mat_SeqBAIJ *)A->data;
2352:   PetscInt           i, j, k, count, *rows;
2353:   PetscInt           bs = A->rmap->bs, bs2 = baij->bs2, *sizes, row, bs_max;
2354:   PetscScalar        zero = 0.0;
2355:   MatScalar         *aa;
2356:   const PetscScalar *xx;
2357:   PetscScalar       *bb;

2359:   PetscFunctionBegin;
2360:   /* fix right-hand side if needed */
2361:   if (x && b) {
2362:     PetscCall(VecGetArrayRead(x, &xx));
2363:     PetscCall(VecGetArray(b, &bb));
2364:     for (i = 0; i < is_n; i++) bb[is_idx[i]] = diag * xx[is_idx[i]];
2365:     PetscCall(VecRestoreArrayRead(x, &xx));
2366:     PetscCall(VecRestoreArray(b, &bb));
2367:   }

2369:   /* Make a copy of the IS and  sort it */
2370:   /* allocate memory for rows,sizes */
2371:   PetscCall(PetscMalloc2(is_n, &rows, 2 * is_n, &sizes));

2373:   /* copy IS values to rows, and sort them */
2374:   for (i = 0; i < is_n; i++) rows[i] = is_idx[i];
2375:   PetscCall(PetscSortInt(is_n, rows));

2377:   if (baij->keepnonzeropattern) {
2378:     for (i = 0; i < is_n; i++) sizes[i] = 1;
2379:     bs_max = is_n;
2380:   } else {
2381:     PetscCall(MatZeroRows_SeqBAIJ_Check_Blocks(rows, is_n, bs, sizes, &bs_max));
2382:     A->nonzerostate++;
2383:   }

2385:   for (i = 0, j = 0; i < bs_max; j += sizes[i], i++) {
2386:     row = rows[j];
2387:     PetscCheck(row >= 0 && row <= A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", row);
2388:     count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
2389:     aa    = ((MatScalar *)baij->a) + baij->i[row / bs] * bs2 + (row % bs);
2390:     if (sizes[i] == bs && !baij->keepnonzeropattern) {
2391:       if (diag != (PetscScalar)0.0) {
2392:         if (baij->ilen[row / bs] > 0) {
2393:           baij->ilen[row / bs]       = 1;
2394:           baij->j[baij->i[row / bs]] = row / bs;

2396:           PetscCall(PetscArrayzero(aa, count * bs));
2397:         }
2398:         /* Now insert all the diagonal values for this bs */
2399:         for (k = 0; k < bs; k++) PetscUseTypeMethod(A, setvalues, 1, rows + j + k, 1, rows + j + k, &diag, INSERT_VALUES);
2400:       } else { /* (diag == 0.0) */
2401:         baij->ilen[row / bs] = 0;
2402:       } /* end (diag == 0.0) */
2403:     } else { /* (sizes[i] != bs) */
2404:       PetscAssert(sizes[i] == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Internal Error. Value should be 1");
2405:       for (k = 0; k < count; k++) {
2406:         aa[0] = zero;
2407:         aa += bs;
2408:       }
2409:       if (diag != (PetscScalar)0.0) PetscUseTypeMethod(A, setvalues, 1, rows + j, 1, rows + j, &diag, INSERT_VALUES);
2410:     }
2411:   }

2413:   PetscCall(PetscFree2(rows, sizes));
2414:   PetscCall(MatAssemblyEnd_SeqBAIJ(A, MAT_FINAL_ASSEMBLY));
2415:   PetscFunctionReturn(PETSC_SUCCESS);
2416: }

2418: static PetscErrorCode MatZeroRowsColumns_SeqBAIJ(Mat A, PetscInt is_n, const PetscInt is_idx[], PetscScalar diag, Vec x, Vec b)
2419: {
2420:   Mat_SeqBAIJ       *baij = (Mat_SeqBAIJ *)A->data;
2421:   PetscInt           i, j, k, count;
2422:   PetscInt           bs = A->rmap->bs, bs2 = baij->bs2, row, col;
2423:   PetscScalar        zero = 0.0;
2424:   MatScalar         *aa;
2425:   const PetscScalar *xx;
2426:   PetscScalar       *bb;
2427:   PetscBool         *zeroed, vecs = PETSC_FALSE;

2429:   PetscFunctionBegin;
2430:   /* fix right-hand side if needed */
2431:   if (x && b) {
2432:     PetscCall(VecGetArrayRead(x, &xx));
2433:     PetscCall(VecGetArray(b, &bb));
2434:     vecs = PETSC_TRUE;
2435:   }

2437:   /* zero the columns */
2438:   PetscCall(PetscCalloc1(A->rmap->n, &zeroed));
2439:   for (i = 0; i < is_n; i++) {
2440:     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]);
2441:     zeroed[is_idx[i]] = PETSC_TRUE;
2442:   }
2443:   for (i = 0; i < A->rmap->N; i++) {
2444:     if (!zeroed[i]) {
2445:       row = i / bs;
2446:       for (j = baij->i[row]; j < baij->i[row + 1]; j++) {
2447:         for (k = 0; k < bs; k++) {
2448:           col = bs * baij->j[j] + k;
2449:           if (zeroed[col]) {
2450:             aa = ((MatScalar *)baij->a) + j * bs2 + (i % bs) + bs * k;
2451:             if (vecs) bb[i] -= aa[0] * xx[col];
2452:             aa[0] = 0.0;
2453:           }
2454:         }
2455:       }
2456:     } else if (vecs) bb[i] = diag * xx[i];
2457:   }
2458:   PetscCall(PetscFree(zeroed));
2459:   if (vecs) {
2460:     PetscCall(VecRestoreArrayRead(x, &xx));
2461:     PetscCall(VecRestoreArray(b, &bb));
2462:   }

2464:   /* zero the rows */
2465:   for (i = 0; i < is_n; i++) {
2466:     row   = is_idx[i];
2467:     count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
2468:     aa    = ((MatScalar *)baij->a) + baij->i[row / bs] * bs2 + (row % bs);
2469:     for (k = 0; k < count; k++) {
2470:       aa[0] = zero;
2471:       aa += bs;
2472:     }
2473:     if (diag != (PetscScalar)0.0) PetscUseTypeMethod(A, setvalues, 1, &row, 1, &row, &diag, INSERT_VALUES);
2474:   }
2475:   PetscCall(MatAssemblyEnd_SeqBAIJ(A, MAT_FINAL_ASSEMBLY));
2476:   PetscFunctionReturn(PETSC_SUCCESS);
2477: }

2479: PetscErrorCode MatSetValues_SeqBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
2480: {
2481:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2482:   PetscInt    *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N, lastcol = -1;
2483:   PetscInt    *imax = a->imax, *ai = a->i, *ailen = a->ilen;
2484:   PetscInt    *aj = a->j, nonew = a->nonew, bs = A->rmap->bs, brow, bcol;
2485:   PetscInt     ridx, cidx, bs2                 = a->bs2;
2486:   PetscBool    roworiented = a->roworiented;
2487:   MatScalar   *ap = NULL, value = 0.0, *aa = a->a, *bap;

2489:   PetscFunctionBegin;
2490:   for (k = 0; k < m; k++) { /* loop over added rows */
2491:     row  = im[k];
2492:     brow = row / bs;
2493:     if (row < 0) continue;
2494:     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);
2495:     rp = PetscSafePointerPlusOffset(aj, ai[brow]);
2496:     if (!A->structure_only) ap = PetscSafePointerPlusOffset(aa, bs2 * ai[brow]);
2497:     rmax = imax[brow];
2498:     nrow = ailen[brow];
2499:     low  = 0;
2500:     high = nrow;
2501:     for (l = 0; l < n; l++) { /* loop over added columns */
2502:       if (in[l] < 0) continue;
2503:       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);
2504:       col  = in[l];
2505:       bcol = col / bs;
2506:       ridx = row % bs;
2507:       cidx = col % bs;
2508:       if (!A->structure_only) {
2509:         if (roworiented) {
2510:           value = v[l + k * n];
2511:         } else {
2512:           value = v[k + l * m];
2513:         }
2514:       }
2515:       if (col <= lastcol) low = 0;
2516:       else high = nrow;
2517:       lastcol = col;
2518:       while (high - low > 7) {
2519:         t = (low + high) / 2;
2520:         if (rp[t] > bcol) high = t;
2521:         else low = t;
2522:       }
2523:       for (i = low; i < high; i++) {
2524:         if (rp[i] > bcol) break;
2525:         if (rp[i] == bcol) {
2526:           bap = PetscSafePointerPlusOffset(ap, bs2 * i + bs * cidx + ridx);
2527:           if (!A->structure_only) {
2528:             if (is == ADD_VALUES) *bap += value;
2529:             else *bap = value;
2530:           }
2531:           goto noinsert1;
2532:         }
2533:       }
2534:       if (nonew == 1) goto noinsert1;
2535:       PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col);
2536:       if (A->structure_only) {
2537:         MatSeqXAIJReallocateAIJ_structure_only(A, a->mbs, bs2, nrow, brow, bcol, rmax, ai, aj, rp, imax, nonew, MatScalar);
2538:       } else {
2539:         MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
2540:       }
2541:       N = nrow++ - 1;
2542:       high++;
2543:       /* shift up all the later entries in this row */
2544:       PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
2545:       rp[i] = bcol;
2546:       if (!A->structure_only) {
2547:         PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
2548:         PetscCall(PetscArrayzero(ap + bs2 * i, bs2));
2549:         ap[bs2 * i + bs * cidx + ridx] = value;
2550:       }
2551:       a->nz++;
2552:       A->nonzerostate++;
2553:     noinsert1:;
2554:       low = i;
2555:     }
2556:     ailen[brow] = nrow;
2557:   }
2558:   PetscFunctionReturn(PETSC_SUCCESS);
2559: }

2561: static PetscErrorCode MatILUFactor_SeqBAIJ(Mat inA, IS row, IS col, const MatFactorInfo *info)
2562: {
2563:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)inA->data;
2564:   Mat          outA;
2565:   PetscBool    row_identity, col_identity;

2567:   PetscFunctionBegin;
2568:   PetscCheck(info->levels == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only levels = 0 supported for in-place ILU");
2569:   PetscCall(ISIdentity(row, &row_identity));
2570:   PetscCall(ISIdentity(col, &col_identity));
2571:   PetscCheck(row_identity && col_identity, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Row and column permutations must be identity for in-place ILU");

2573:   outA            = inA;
2574:   inA->factortype = MAT_FACTOR_LU;
2575:   PetscCall(PetscFree(inA->solvertype));
2576:   PetscCall(PetscStrallocpy(MATSOLVERPETSC, &inA->solvertype));

2578:   PetscCall(MatMarkDiagonal_SeqBAIJ(inA));

2580:   PetscCall(PetscObjectReference((PetscObject)row));
2581:   PetscCall(ISDestroy(&a->row));
2582:   a->row = row;
2583:   PetscCall(PetscObjectReference((PetscObject)col));
2584:   PetscCall(ISDestroy(&a->col));
2585:   a->col = col;

2587:   /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */
2588:   PetscCall(ISDestroy(&a->icol));
2589:   PetscCall(ISInvertPermutation(col, PETSC_DECIDE, &a->icol));

2591:   PetscCall(MatSeqBAIJSetNumericFactorization_inplace(inA, (PetscBool)(row_identity && col_identity)));
2592:   if (!a->solve_work) PetscCall(PetscMalloc1(inA->rmap->N + inA->rmap->bs, &a->solve_work));
2593:   PetscCall(MatLUFactorNumeric(outA, inA, info));
2594:   PetscFunctionReturn(PETSC_SUCCESS);
2595: }

2597: static PetscErrorCode MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat, const PetscInt *indices)
2598: {
2599:   Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data;

2601:   PetscFunctionBegin;
2602:   baij->nz = baij->maxnz;
2603:   PetscCall(PetscArraycpy(baij->j, indices, baij->nz));
2604:   PetscCall(PetscArraycpy(baij->ilen, baij->imax, baij->mbs));
2605:   PetscFunctionReturn(PETSC_SUCCESS);
2606: }

2608: /*@
2609:   MatSeqBAIJSetColumnIndices - Set the column indices for all the block rows in the matrix.

2611:   Input Parameters:
2612: + mat     - the `MATSEQBAIJ` matrix
2613: - indices - the block column indices

2615:   Level: advanced

2617:   Notes:
2618:   This can be called if you have precomputed the nonzero structure of the
2619:   matrix and want to provide it to the matrix object to improve the performance
2620:   of the `MatSetValues()` operation.

2622:   You MUST have set the correct numbers of nonzeros per row in the call to
2623:   `MatCreateSeqBAIJ()`, and the columns indices MUST be sorted.

2625:   MUST be called before any calls to `MatSetValues()`

2627: .seealso: [](ch_matrices), `Mat`, `MATSEQBAIJ`, `MatSetValues()`
2628: @*/
2629: PetscErrorCode MatSeqBAIJSetColumnIndices(Mat mat, PetscInt *indices)
2630: {
2631:   PetscFunctionBegin;
2633:   PetscAssertPointer(indices, 2);
2634:   PetscUseMethod(mat, "MatSeqBAIJSetColumnIndices_C", (Mat, const PetscInt *), (mat, (const PetscInt *)indices));
2635:   PetscFunctionReturn(PETSC_SUCCESS);
2636: }

2638: static PetscErrorCode MatGetRowMaxAbs_SeqBAIJ(Mat A, Vec v, PetscInt idx[])
2639: {
2640:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2641:   PetscInt     i, j, n, row, bs, *ai, *aj, mbs;
2642:   PetscReal    atmp;
2643:   PetscScalar *x, zero = 0.0;
2644:   MatScalar   *aa;
2645:   PetscInt     ncols, brow, krow, kcol;

2647:   PetscFunctionBegin;
2648:   /* why is this not a macro???????????????????????????????????????????????????????????????? */
2649:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2650:   bs  = A->rmap->bs;
2651:   aa  = a->a;
2652:   ai  = a->i;
2653:   aj  = a->j;
2654:   mbs = a->mbs;

2656:   PetscCall(VecSet(v, zero));
2657:   PetscCall(VecGetArray(v, &x));
2658:   PetscCall(VecGetLocalSize(v, &n));
2659:   PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2660:   for (i = 0; i < mbs; i++) {
2661:     ncols = ai[1] - ai[0];
2662:     ai++;
2663:     brow = bs * i;
2664:     for (j = 0; j < ncols; j++) {
2665:       for (kcol = 0; kcol < bs; kcol++) {
2666:         for (krow = 0; krow < bs; krow++) {
2667:           atmp = PetscAbsScalar(*aa);
2668:           aa++;
2669:           row = brow + krow; /* row index */
2670:           if (PetscAbsScalar(x[row]) < atmp) {
2671:             x[row] = atmp;
2672:             if (idx) idx[row] = bs * (*aj) + kcol;
2673:           }
2674:         }
2675:       }
2676:       aj++;
2677:     }
2678:   }
2679:   PetscCall(VecRestoreArray(v, &x));
2680:   PetscFunctionReturn(PETSC_SUCCESS);
2681: }

2683: static PetscErrorCode MatGetRowSumAbs_SeqBAIJ(Mat A, Vec v)
2684: {
2685:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2686:   PetscInt     i, j, n, row, bs, *ai, mbs;
2687:   PetscReal    atmp;
2688:   PetscScalar *x, zero = 0.0;
2689:   MatScalar   *aa;
2690:   PetscInt     ncols, brow, krow, kcol;

2692:   PetscFunctionBegin;
2693:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2694:   bs  = A->rmap->bs;
2695:   aa  = a->a;
2696:   ai  = a->i;
2697:   mbs = a->mbs;

2699:   PetscCall(VecSet(v, zero));
2700:   PetscCall(VecGetArrayWrite(v, &x));
2701:   PetscCall(VecGetLocalSize(v, &n));
2702:   PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2703:   for (i = 0; i < mbs; i++) {
2704:     ncols = ai[1] - ai[0];
2705:     ai++;
2706:     brow = bs * i;
2707:     for (j = 0; j < ncols; j++) {
2708:       for (kcol = 0; kcol < bs; kcol++) {
2709:         for (krow = 0; krow < bs; krow++) {
2710:           atmp = PetscAbsScalar(*aa);
2711:           aa++;
2712:           row = brow + krow; /* row index */
2713:           x[row] += atmp;
2714:         }
2715:       }
2716:     }
2717:   }
2718:   PetscCall(VecRestoreArrayWrite(v, &x));
2719:   PetscFunctionReturn(PETSC_SUCCESS);
2720: }

2722: static PetscErrorCode MatCopy_SeqBAIJ(Mat A, Mat B, MatStructure str)
2723: {
2724:   PetscFunctionBegin;
2725:   /* If the two matrices have the same copy implementation, use fast copy. */
2726:   if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
2727:     Mat_SeqBAIJ *a    = (Mat_SeqBAIJ *)A->data;
2728:     Mat_SeqBAIJ *b    = (Mat_SeqBAIJ *)B->data;
2729:     PetscInt     ambs = a->mbs, bmbs = b->mbs, abs = A->rmap->bs, bbs = B->rmap->bs, bs2 = abs * abs;

2731:     PetscCheck(a->i[ambs] == b->i[bmbs], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of nonzero blocks in matrices A %" PetscInt_FMT " and B %" PetscInt_FMT " are different", a->i[ambs], b->i[bmbs]);
2732:     PetscCheck(abs == bbs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Block size A %" PetscInt_FMT " and B %" PetscInt_FMT " are different", abs, bbs);
2733:     PetscCall(PetscArraycpy(b->a, a->a, bs2 * a->i[ambs]));
2734:     PetscCall(PetscObjectStateIncrease((PetscObject)B));
2735:   } else {
2736:     PetscCall(MatCopy_Basic(A, B, str));
2737:   }
2738:   PetscFunctionReturn(PETSC_SUCCESS);
2739: }

2741: static PetscErrorCode MatSeqBAIJGetArray_SeqBAIJ(Mat A, PetscScalar *array[])
2742: {
2743:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;

2745:   PetscFunctionBegin;
2746:   *array = a->a;
2747:   PetscFunctionReturn(PETSC_SUCCESS);
2748: }

2750: static PetscErrorCode MatSeqBAIJRestoreArray_SeqBAIJ(Mat A, PetscScalar *array[])
2751: {
2752:   PetscFunctionBegin;
2753:   *array = NULL;
2754:   PetscFunctionReturn(PETSC_SUCCESS);
2755: }

2757: PetscErrorCode MatAXPYGetPreallocation_SeqBAIJ(Mat Y, Mat X, PetscInt *nnz)
2758: {
2759:   PetscInt     bs = Y->rmap->bs, mbs = Y->rmap->N / bs;
2760:   Mat_SeqBAIJ *x = (Mat_SeqBAIJ *)X->data;
2761:   Mat_SeqBAIJ *y = (Mat_SeqBAIJ *)Y->data;

2763:   PetscFunctionBegin;
2764:   /* Set the number of nonzeros in the new matrix */
2765:   PetscCall(MatAXPYGetPreallocation_SeqX_private(mbs, x->i, x->j, y->i, y->j, nnz));
2766:   PetscFunctionReturn(PETSC_SUCCESS);
2767: }

2769: PetscErrorCode MatAXPY_SeqBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str)
2770: {
2771:   Mat_SeqBAIJ *x = (Mat_SeqBAIJ *)X->data, *y = (Mat_SeqBAIJ *)Y->data;
2772:   PetscInt     bs = Y->rmap->bs, bs2 = bs * bs;
2773:   PetscBLASInt one = 1;

2775:   PetscFunctionBegin;
2776:   if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) {
2777:     PetscBool e = x->nz == y->nz && x->mbs == y->mbs && bs == X->rmap->bs ? PETSC_TRUE : PETSC_FALSE;
2778:     if (e) {
2779:       PetscCall(PetscArraycmp(x->i, y->i, x->mbs + 1, &e));
2780:       if (e) {
2781:         PetscCall(PetscArraycmp(x->j, y->j, x->i[x->mbs], &e));
2782:         if (e) str = SAME_NONZERO_PATTERN;
2783:       }
2784:     }
2785:     if (!e) PetscCheck(str != SAME_NONZERO_PATTERN, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatStructure is not SAME_NONZERO_PATTERN");
2786:   }
2787:   if (str == SAME_NONZERO_PATTERN) {
2788:     PetscScalar  alpha = a;
2789:     PetscBLASInt bnz;
2790:     PetscCall(PetscBLASIntCast(x->nz * bs2, &bnz));
2791:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one));
2792:     PetscCall(PetscObjectStateIncrease((PetscObject)Y));
2793:   } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
2794:     PetscCall(MatAXPY_Basic(Y, a, X, str));
2795:   } else {
2796:     Mat       B;
2797:     PetscInt *nnz;
2798:     PetscCheck(bs == X->rmap->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrices must have same block size");
2799:     PetscCall(PetscMalloc1(Y->rmap->N, &nnz));
2800:     PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B));
2801:     PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name));
2802:     PetscCall(MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N));
2803:     PetscCall(MatSetBlockSizesFromMats(B, Y, Y));
2804:     PetscCall(MatSetType(B, (MatType)((PetscObject)Y)->type_name));
2805:     PetscCall(MatAXPYGetPreallocation_SeqBAIJ(Y, X, nnz));
2806:     PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, nnz));
2807:     PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
2808:     PetscCall(MatHeaderMerge(Y, &B));
2809:     PetscCall(PetscFree(nnz));
2810:   }
2811:   PetscFunctionReturn(PETSC_SUCCESS);
2812: }

2814: PETSC_INTERN PetscErrorCode MatConjugate_SeqBAIJ(Mat A)
2815: {
2816: #if PetscDefined(USE_COMPLEX)
2817:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2818:   PetscInt     i, nz = a->bs2 * a->i[a->mbs];
2819:   MatScalar   *aa = a->a;

2821:   PetscFunctionBegin;
2822:   for (i = 0; i < nz; i++) aa[i] = PetscConj(aa[i]);
2823:   PetscFunctionReturn(PETSC_SUCCESS);
2824: #else
2825:   (void)A;
2826:   return PETSC_SUCCESS;
2827: #endif
2828: }

2830: static PetscErrorCode MatRealPart_SeqBAIJ(Mat A)
2831: {
2832: #if PetscDefined(USE_COMPLEX)
2833:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2834:   PetscInt     i, nz = a->bs2 * a->i[a->mbs];
2835:   MatScalar   *aa = a->a;

2837:   PetscFunctionBegin;
2838:   for (i = 0; i < nz; i++) aa[i] = PetscRealPart(aa[i]);
2839:   PetscFunctionReturn(PETSC_SUCCESS);
2840: #else
2841:   (void)A;
2842:   return PETSC_SUCCESS;
2843: #endif
2844: }

2846: static PetscErrorCode MatImaginaryPart_SeqBAIJ(Mat A)
2847: {
2848: #if PetscDefined(USE_COMPLEX)
2849:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2850:   PetscInt     i, nz = a->bs2 * a->i[a->mbs];
2851:   MatScalar   *aa = a->a;

2853:   PetscFunctionBegin;
2854:   for (i = 0; i < nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2855:   PetscFunctionReturn(PETSC_SUCCESS);
2856: #else
2857:   (void)A;
2858:   return PETSC_SUCCESS;
2859: #endif
2860: }

2862: /*
2863:     Code almost identical to MatGetColumnIJ_SeqAIJ() should share common code
2864: */
2865: static PetscErrorCode MatGetColumnIJ_SeqBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done)
2866: {
2867:   Mat_SeqBAIJ *a  = (Mat_SeqBAIJ *)A->data;
2868:   PetscInt     bs = A->rmap->bs, i, *collengths, *cia, *cja, n = A->cmap->n / bs, m = A->rmap->n / bs;
2869:   PetscInt     nz = a->i[m], row, *jj, mr, col;

2871:   PetscFunctionBegin;
2872:   *nn = n;
2873:   if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
2874:   PetscCheck(!symmetric, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for BAIJ matrices");
2875:   PetscCall(PetscCalloc1(n, &collengths));
2876:   PetscCall(PetscMalloc1(n + 1, &cia));
2877:   PetscCall(PetscMalloc1(nz, &cja));
2878:   jj = a->j;
2879:   for (i = 0; i < nz; i++) collengths[jj[i]]++;
2880:   cia[0] = oshift;
2881:   for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i];
2882:   PetscCall(PetscArrayzero(collengths, n));
2883:   jj = a->j;
2884:   for (row = 0; row < m; row++) {
2885:     mr = a->i[row + 1] - a->i[row];
2886:     for (i = 0; i < mr; i++) {
2887:       col = *jj++;

2889:       cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
2890:     }
2891:   }
2892:   PetscCall(PetscFree(collengths));
2893:   *ia = cia;
2894:   *ja = cja;
2895:   PetscFunctionReturn(PETSC_SUCCESS);
2896: }

2898: static PetscErrorCode MatRestoreColumnIJ_SeqBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done)
2899: {
2900:   PetscFunctionBegin;
2901:   if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
2902:   PetscCall(PetscFree(*ia));
2903:   PetscCall(PetscFree(*ja));
2904:   PetscFunctionReturn(PETSC_SUCCESS);
2905: }

2907: /*
2908:  MatGetColumnIJ_SeqBAIJ_Color() and MatRestoreColumnIJ_SeqBAIJ_Color() are customized from
2909:  MatGetColumnIJ_SeqBAIJ() and MatRestoreColumnIJ_SeqBAIJ() by adding an output
2910:  spidx[], index of a->a, to be used in MatTransposeColoringCreate() and MatFDColoringCreate()
2911:  */
2912: PetscErrorCode MatGetColumnIJ_SeqBAIJ_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done)
2913: {
2914:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2915:   PetscInt     i, *collengths, *cia, *cja, n = a->nbs, m = a->mbs;
2916:   PetscInt     nz = a->i[m], row, *jj, mr, col;
2917:   PetscInt    *cspidx;

2919:   PetscFunctionBegin;
2920:   *nn = n;
2921:   if (!ia) PetscFunctionReturn(PETSC_SUCCESS);

2923:   PetscCall(PetscCalloc1(n, &collengths));
2924:   PetscCall(PetscMalloc1(n + 1, &cia));
2925:   PetscCall(PetscMalloc1(nz, &cja));
2926:   PetscCall(PetscMalloc1(nz, &cspidx));
2927:   jj = a->j;
2928:   for (i = 0; i < nz; i++) collengths[jj[i]]++;
2929:   cia[0] = oshift;
2930:   for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i];
2931:   PetscCall(PetscArrayzero(collengths, n));
2932:   jj = a->j;
2933:   for (row = 0; row < m; row++) {
2934:     mr = a->i[row + 1] - a->i[row];
2935:     for (i = 0; i < mr; i++) {
2936:       col                                         = *jj++;
2937:       cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */
2938:       cja[cia[col] + collengths[col]++ - oshift]  = row + oshift;
2939:     }
2940:   }
2941:   PetscCall(PetscFree(collengths));
2942:   *ia    = cia;
2943:   *ja    = cja;
2944:   *spidx = cspidx;
2945:   PetscFunctionReturn(PETSC_SUCCESS);
2946: }

2948: PetscErrorCode MatRestoreColumnIJ_SeqBAIJ_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done)
2949: {
2950:   PetscFunctionBegin;
2951:   PetscCall(MatRestoreColumnIJ_SeqBAIJ(A, oshift, symmetric, inodecompressed, n, ia, ja, done));
2952:   PetscCall(PetscFree(*spidx));
2953:   PetscFunctionReturn(PETSC_SUCCESS);
2954: }

2956: static PetscErrorCode MatShift_SeqBAIJ(Mat Y, PetscScalar a)
2957: {
2958:   Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)Y->data;

2960:   PetscFunctionBegin;
2961:   if (!Y->preallocated || !aij->nz) PetscCall(MatSeqBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL));
2962:   PetscCall(MatShift_Basic(Y, a));
2963:   PetscFunctionReturn(PETSC_SUCCESS);
2964: }

2966: PetscErrorCode MatEliminateZeros_SeqBAIJ(Mat A, PetscBool keep)
2967: {
2968:   Mat_SeqBAIJ *a      = (Mat_SeqBAIJ *)A->data;
2969:   PetscInt     fshift = 0, fshift_prev = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax, j, k;
2970:   PetscInt     m = A->rmap->N, *ailen = a->ilen;
2971:   PetscInt     mbs = a->mbs, bs2 = a->bs2, rmax = 0;
2972:   MatScalar   *aa = a->a, *ap;
2973:   PetscBool    zero;

2975:   PetscFunctionBegin;
2976:   PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot eliminate zeros for unassembled matrix");
2977:   if (m) rmax = ailen[0];
2978:   for (i = 1; i <= mbs; i++) {
2979:     for (k = ai[i - 1]; k < ai[i]; k++) {
2980:       zero = PETSC_TRUE;
2981:       ap   = aa + bs2 * k;
2982:       for (j = 0; j < bs2 && zero; j++) {
2983:         if (ap[j] != 0.0) zero = PETSC_FALSE;
2984:       }
2985:       if (zero && (aj[k] != i - 1 || !keep)) fshift++;
2986:       else {
2987:         if (zero && aj[k] == i - 1) PetscCall(PetscInfo(A, "Keep the diagonal block at row %" PetscInt_FMT "\n", i - 1));
2988:         aj[k - fshift] = aj[k];
2989:         PetscCall(PetscArraymove(ap - bs2 * fshift, ap, bs2));
2990:       }
2991:     }
2992:     ai[i - 1] -= fshift_prev;
2993:     fshift_prev  = fshift;
2994:     ailen[i - 1] = imax[i - 1] = ai[i] - fshift - ai[i - 1];
2995:     a->nonzerorowcnt += ((ai[i] - fshift - ai[i - 1]) > 0);
2996:     rmax = PetscMax(rmax, ailen[i - 1]);
2997:   }
2998:   if (fshift) {
2999:     if (mbs) {
3000:       ai[mbs] -= fshift;
3001:       a->nz = ai[mbs];
3002:     }
3003:     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));
3004:     A->nonzerostate++;
3005:     A->info.nz_unneeded += (PetscReal)fshift;
3006:     a->rmax = rmax;
3007:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
3008:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
3009:   }
3010:   PetscFunctionReturn(PETSC_SUCCESS);
3011: }

3013: static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ,
3014:                                        MatGetRow_SeqBAIJ,
3015:                                        MatRestoreRow_SeqBAIJ,
3016:                                        MatMult_SeqBAIJ_N,
3017:                                        /* 4*/ MatMultAdd_SeqBAIJ_N,
3018:                                        MatMultTranspose_SeqBAIJ,
3019:                                        MatMultTransposeAdd_SeqBAIJ,
3020:                                        NULL,
3021:                                        NULL,
3022:                                        NULL,
3023:                                        /* 10*/ NULL,
3024:                                        MatLUFactor_SeqBAIJ,
3025:                                        NULL,
3026:                                        NULL,
3027:                                        MatTranspose_SeqBAIJ,
3028:                                        /* 15*/ MatGetInfo_SeqBAIJ,
3029:                                        MatEqual_SeqBAIJ,
3030:                                        MatGetDiagonal_SeqBAIJ,
3031:                                        MatDiagonalScale_SeqBAIJ,
3032:                                        MatNorm_SeqBAIJ,
3033:                                        /* 20*/ NULL,
3034:                                        MatAssemblyEnd_SeqBAIJ,
3035:                                        MatSetOption_SeqBAIJ,
3036:                                        MatZeroEntries_SeqBAIJ,
3037:                                        /* 24*/ MatZeroRows_SeqBAIJ,
3038:                                        NULL,
3039:                                        NULL,
3040:                                        NULL,
3041:                                        NULL,
3042:                                        /* 29*/ MatSetUp_Seq_Hash,
3043:                                        NULL,
3044:                                        NULL,
3045:                                        NULL,
3046:                                        NULL,
3047:                                        /* 34*/ MatDuplicate_SeqBAIJ,
3048:                                        NULL,
3049:                                        NULL,
3050:                                        MatILUFactor_SeqBAIJ,
3051:                                        NULL,
3052:                                        /* 39*/ MatAXPY_SeqBAIJ,
3053:                                        MatCreateSubMatrices_SeqBAIJ,
3054:                                        MatIncreaseOverlap_SeqBAIJ,
3055:                                        MatGetValues_SeqBAIJ,
3056:                                        MatCopy_SeqBAIJ,
3057:                                        /* 44*/ NULL,
3058:                                        MatScale_SeqBAIJ,
3059:                                        MatShift_SeqBAIJ,
3060:                                        NULL,
3061:                                        MatZeroRowsColumns_SeqBAIJ,
3062:                                        /* 49*/ NULL,
3063:                                        MatGetRowIJ_SeqBAIJ,
3064:                                        MatRestoreRowIJ_SeqBAIJ,
3065:                                        MatGetColumnIJ_SeqBAIJ,
3066:                                        MatRestoreColumnIJ_SeqBAIJ,
3067:                                        /* 54*/ MatFDColoringCreate_SeqXAIJ,
3068:                                        NULL,
3069:                                        NULL,
3070:                                        NULL,
3071:                                        MatSetValuesBlocked_SeqBAIJ,
3072:                                        /* 59*/ MatCreateSubMatrix_SeqBAIJ,
3073:                                        MatDestroy_SeqBAIJ,
3074:                                        MatView_SeqBAIJ,
3075:                                        NULL,
3076:                                        NULL,
3077:                                        /* 64*/ NULL,
3078:                                        NULL,
3079:                                        NULL,
3080:                                        NULL,
3081:                                        NULL,
3082:                                        /* 69*/ MatGetRowMaxAbs_SeqBAIJ,
3083:                                        NULL,
3084:                                        MatConvert_Basic,
3085:                                        NULL,
3086:                                        NULL,
3087:                                        /* 74*/ NULL,
3088:                                        MatFDColoringApply_BAIJ,
3089:                                        NULL,
3090:                                        NULL,
3091:                                        NULL,
3092:                                        /* 79*/ NULL,
3093:                                        NULL,
3094:                                        NULL,
3095:                                        NULL,
3096:                                        MatLoad_SeqBAIJ,
3097:                                        /* 84*/ NULL,
3098:                                        NULL,
3099:                                        NULL,
3100:                                        NULL,
3101:                                        NULL,
3102:                                        /* 89*/ NULL,
3103:                                        NULL,
3104:                                        NULL,
3105:                                        NULL,
3106:                                        NULL,
3107:                                        /* 94*/ NULL,
3108:                                        NULL,
3109:                                        NULL,
3110:                                        NULL,
3111:                                        NULL,
3112:                                        /* 99*/ NULL,
3113:                                        NULL,
3114:                                        NULL,
3115:                                        MatConjugate_SeqBAIJ,
3116:                                        NULL,
3117:                                        /*104*/ NULL,
3118:                                        MatRealPart_SeqBAIJ,
3119:                                        MatImaginaryPart_SeqBAIJ,
3120:                                        NULL,
3121:                                        NULL,
3122:                                        /*109*/ NULL,
3123:                                        NULL,
3124:                                        NULL,
3125:                                        NULL,
3126:                                        MatMissingDiagonal_SeqBAIJ,
3127:                                        /*114*/ NULL,
3128:                                        NULL,
3129:                                        NULL,
3130:                                        NULL,
3131:                                        NULL,
3132:                                        /*119*/ NULL,
3133:                                        NULL,
3134:                                        MatMultHermitianTranspose_SeqBAIJ,
3135:                                        MatMultHermitianTransposeAdd_SeqBAIJ,
3136:                                        NULL,
3137:                                        /*124*/ NULL,
3138:                                        MatGetColumnReductions_SeqBAIJ,
3139:                                        MatInvertBlockDiagonal_SeqBAIJ,
3140:                                        NULL,
3141:                                        NULL,
3142:                                        /*129*/ NULL,
3143:                                        NULL,
3144:                                        NULL,
3145:                                        NULL,
3146:                                        NULL,
3147:                                        /*134*/ NULL,
3148:                                        NULL,
3149:                                        NULL,
3150:                                        NULL,
3151:                                        NULL,
3152:                                        /*139*/ MatSetBlockSizes_Default,
3153:                                        NULL,
3154:                                        NULL,
3155:                                        MatFDColoringSetUp_SeqXAIJ,
3156:                                        NULL,
3157:                                        /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqBAIJ,
3158:                                        MatDestroySubMatrices_SeqBAIJ,
3159:                                        NULL,
3160:                                        NULL,
3161:                                        NULL,
3162:                                        NULL,
3163:                                        /*150*/ NULL,
3164:                                        MatEliminateZeros_SeqBAIJ,
3165:                                        MatGetRowSumAbs_SeqBAIJ};

3167: static PetscErrorCode MatStoreValues_SeqBAIJ(Mat mat)
3168: {
3169:   Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data;
3170:   PetscInt     nz  = aij->i[aij->mbs] * aij->bs2;

3172:   PetscFunctionBegin;
3173:   PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");

3175:   /* allocate space for values if not already there */
3176:   if (!aij->saved_values) PetscCall(PetscMalloc1(nz + 1, &aij->saved_values));

3178:   /* copy values over */
3179:   PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz));
3180:   PetscFunctionReturn(PETSC_SUCCESS);
3181: }

3183: static PetscErrorCode MatRetrieveValues_SeqBAIJ(Mat mat)
3184: {
3185:   Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data;
3186:   PetscInt     nz  = aij->i[aij->mbs] * aij->bs2;

3188:   PetscFunctionBegin;
3189:   PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3190:   PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first");

3192:   /* copy values over */
3193:   PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz));
3194:   PetscFunctionReturn(PETSC_SUCCESS);
3195: }

3197: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *);
3198: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat, MatType, MatReuse, Mat *);

3200: PetscErrorCode MatSeqBAIJSetPreallocation_SeqBAIJ(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[])
3201: {
3202:   Mat_SeqBAIJ *b = (Mat_SeqBAIJ *)B->data;
3203:   PetscInt     i, mbs, nbs, bs2;
3204:   PetscBool    flg = PETSC_FALSE, skipallocation = PETSC_FALSE, realalloc = PETSC_FALSE;

3206:   PetscFunctionBegin;
3207:   if (B->hash_active) {
3208:     PetscInt bs;
3209:     B->ops[0] = b->cops;
3210:     PetscCall(PetscHMapIJVDestroy(&b->ht));
3211:     PetscCall(MatGetBlockSize(B, &bs));
3212:     if (bs > 1) PetscCall(PetscHSetIJDestroy(&b->bht));
3213:     PetscCall(PetscFree(b->dnz));
3214:     PetscCall(PetscFree(b->bdnz));
3215:     B->hash_active = PETSC_FALSE;
3216:   }
3217:   if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
3218:   if (nz == MAT_SKIP_ALLOCATION) {
3219:     skipallocation = PETSC_TRUE;
3220:     nz             = 0;
3221:   }

3223:   PetscCall(MatSetBlockSize(B, PetscAbs(bs)));
3224:   PetscCall(PetscLayoutSetUp(B->rmap));
3225:   PetscCall(PetscLayoutSetUp(B->cmap));
3226:   PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));

3228:   B->preallocated = PETSC_TRUE;

3230:   mbs = B->rmap->n / bs;
3231:   nbs = B->cmap->n / bs;
3232:   bs2 = bs * bs;

3234:   PetscCheck(mbs * bs == B->rmap->n && nbs * bs == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows %" PetscInt_FMT ", cols %" PetscInt_FMT " must be divisible by blocksize %" PetscInt_FMT, B->rmap->N, B->cmap->n, bs);

3236:   if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
3237:   PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz);
3238:   if (nnz) {
3239:     for (i = 0; i < mbs; i++) {
3240:       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]);
3241:       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 " rowlength %" PetscInt_FMT, i, nnz[i], nbs);
3242:     }
3243:   }

3245:   PetscOptionsBegin(PetscObjectComm((PetscObject)B), NULL, "Optimize options for SEQBAIJ matrix 2 ", "Mat");
3246:   PetscCall(PetscOptionsBool("-mat_no_unroll", "Do not optimize for block size (slow)", NULL, flg, &flg, NULL));
3247:   PetscOptionsEnd();

3249:   if (!flg) {
3250:     switch (bs) {
3251:     case 1:
3252:       B->ops->mult    = MatMult_SeqBAIJ_1;
3253:       B->ops->multadd = MatMultAdd_SeqBAIJ_1;
3254:       break;
3255:     case 2:
3256:       B->ops->mult    = MatMult_SeqBAIJ_2;
3257:       B->ops->multadd = MatMultAdd_SeqBAIJ_2;
3258:       break;
3259:     case 3:
3260:       B->ops->mult    = MatMult_SeqBAIJ_3;
3261:       B->ops->multadd = MatMultAdd_SeqBAIJ_3;
3262:       break;
3263:     case 4:
3264:       B->ops->mult    = MatMult_SeqBAIJ_4;
3265:       B->ops->multadd = MatMultAdd_SeqBAIJ_4;
3266:       break;
3267:     case 5:
3268:       B->ops->mult    = MatMult_SeqBAIJ_5;
3269:       B->ops->multadd = MatMultAdd_SeqBAIJ_5;
3270:       break;
3271:     case 6:
3272:       B->ops->mult    = MatMult_SeqBAIJ_6;
3273:       B->ops->multadd = MatMultAdd_SeqBAIJ_6;
3274:       break;
3275:     case 7:
3276:       B->ops->mult    = MatMult_SeqBAIJ_7;
3277:       B->ops->multadd = MatMultAdd_SeqBAIJ_7;
3278:       break;
3279:     case 9: {
3280:       PetscInt version = 1;
3281:       PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)B)->prefix, "-mat_baij_mult_version", &version, NULL));
3282:       switch (version) {
3283: #if defined(PETSC_HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && defined(PETSC_USE_REAL_DOUBLE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_64BIT_INDICES)
3284:       case 1:
3285:         B->ops->mult    = MatMult_SeqBAIJ_9_AVX2;
3286:         B->ops->multadd = MatMultAdd_SeqBAIJ_9_AVX2;
3287:         PetscCall(PetscInfo((PetscObject)B, "Using AVX2 for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3288:         break;
3289: #endif
3290:       default:
3291:         B->ops->mult    = MatMult_SeqBAIJ_N;
3292:         B->ops->multadd = MatMultAdd_SeqBAIJ_N;
3293:         PetscCall(PetscInfo((PetscObject)B, "Using BLAS for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3294:         break;
3295:       }
3296:       break;
3297:     }
3298:     case 11:
3299:       B->ops->mult    = MatMult_SeqBAIJ_11;
3300:       B->ops->multadd = MatMultAdd_SeqBAIJ_11;
3301:       break;
3302:     case 12: {
3303:       PetscInt version = 1;
3304:       PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)B)->prefix, "-mat_baij_mult_version", &version, NULL));
3305:       switch (version) {
3306:       case 1:
3307:         B->ops->mult    = MatMult_SeqBAIJ_12_ver1;
3308:         B->ops->multadd = MatMultAdd_SeqBAIJ_12_ver1;
3309:         PetscCall(PetscInfo((PetscObject)B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3310:         break;
3311:       case 2:
3312:         B->ops->mult    = MatMult_SeqBAIJ_12_ver2;
3313:         B->ops->multadd = MatMultAdd_SeqBAIJ_12_ver2;
3314:         PetscCall(PetscInfo((PetscObject)B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3315:         break;
3316: #if defined(PETSC_HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && defined(PETSC_USE_REAL_DOUBLE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_64BIT_INDICES)
3317:       case 3:
3318:         B->ops->mult    = MatMult_SeqBAIJ_12_AVX2;
3319:         B->ops->multadd = MatMultAdd_SeqBAIJ_12_ver1;
3320:         PetscCall(PetscInfo((PetscObject)B, "Using AVX2 for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3321:         break;
3322: #endif
3323:       default:
3324:         B->ops->mult    = MatMult_SeqBAIJ_N;
3325:         B->ops->multadd = MatMultAdd_SeqBAIJ_N;
3326:         PetscCall(PetscInfo((PetscObject)B, "Using BLAS for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3327:         break;
3328:       }
3329:       break;
3330:     }
3331:     case 15: {
3332:       PetscInt version = 1;
3333:       PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)B)->prefix, "-mat_baij_mult_version", &version, NULL));
3334:       switch (version) {
3335:       case 1:
3336:         B->ops->mult = MatMult_SeqBAIJ_15_ver1;
3337:         PetscCall(PetscInfo((PetscObject)B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3338:         break;
3339:       case 2:
3340:         B->ops->mult = MatMult_SeqBAIJ_15_ver2;
3341:         PetscCall(PetscInfo((PetscObject)B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3342:         break;
3343:       case 3:
3344:         B->ops->mult = MatMult_SeqBAIJ_15_ver3;
3345:         PetscCall(PetscInfo((PetscObject)B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3346:         break;
3347:       case 4:
3348:         B->ops->mult = MatMult_SeqBAIJ_15_ver4;
3349:         PetscCall(PetscInfo((PetscObject)B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3350:         break;
3351:       default:
3352:         B->ops->mult = MatMult_SeqBAIJ_N;
3353:         PetscCall(PetscInfo((PetscObject)B, "Using BLAS for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3354:         break;
3355:       }
3356:       B->ops->multadd = MatMultAdd_SeqBAIJ_N;
3357:       break;
3358:     }
3359:     default:
3360:       B->ops->mult    = MatMult_SeqBAIJ_N;
3361:       B->ops->multadd = MatMultAdd_SeqBAIJ_N;
3362:       PetscCall(PetscInfo((PetscObject)B, "Using BLAS for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3363:       break;
3364:     }
3365:   }
3366:   B->ops->sor = MatSOR_SeqBAIJ;
3367:   b->mbs      = mbs;
3368:   b->nbs      = nbs;
3369:   if (!skipallocation) {
3370:     if (!b->imax) {
3371:       PetscCall(PetscMalloc2(mbs, &b->imax, mbs, &b->ilen));

3373:       b->free_imax_ilen = PETSC_TRUE;
3374:     }
3375:     /* b->ilen will count nonzeros in each block row so far. */
3376:     for (i = 0; i < mbs; i++) b->ilen[i] = 0;
3377:     if (!nnz) {
3378:       if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
3379:       else if (nz < 0) nz = 1;
3380:       nz = PetscMin(nz, nbs);
3381:       for (i = 0; i < mbs; i++) b->imax[i] = nz;
3382:       PetscCall(PetscIntMultError(nz, mbs, &nz));
3383:     } else {
3384:       PetscInt64 nz64 = 0;
3385:       for (i = 0; i < mbs; i++) {
3386:         b->imax[i] = nnz[i];
3387:         nz64 += nnz[i];
3388:       }
3389:       PetscCall(PetscIntCast(nz64, &nz));
3390:     }

3392:     /* allocate the matrix space */
3393:     PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i));
3394:     if (B->structure_only) {
3395:       PetscCall(PetscMalloc1(nz, &b->j));
3396:       PetscCall(PetscMalloc1(B->rmap->N + 1, &b->i));
3397:     } else {
3398:       PetscInt nzbs2 = 0;
3399:       PetscCall(PetscIntMultError(nz, bs2, &nzbs2));
3400:       PetscCall(PetscMalloc3(nzbs2, &b->a, nz, &b->j, B->rmap->N + 1, &b->i));
3401:       PetscCall(PetscArrayzero(b->a, nz * bs2));
3402:     }
3403:     PetscCall(PetscArrayzero(b->j, nz));

3405:     if (B->structure_only) {
3406:       b->singlemalloc = PETSC_FALSE;
3407:       b->free_a       = PETSC_FALSE;
3408:     } else {
3409:       b->singlemalloc = PETSC_TRUE;
3410:       b->free_a       = PETSC_TRUE;
3411:     }
3412:     b->free_ij = PETSC_TRUE;

3414:     b->i[0] = 0;
3415:     for (i = 1; i < mbs + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1];

3417:   } else {
3418:     b->free_a  = PETSC_FALSE;
3419:     b->free_ij = PETSC_FALSE;
3420:   }

3422:   b->bs2              = bs2;
3423:   b->mbs              = mbs;
3424:   b->nz               = 0;
3425:   b->maxnz            = nz;
3426:   B->info.nz_unneeded = (PetscReal)b->maxnz * bs2;
3427:   B->was_assembled    = PETSC_FALSE;
3428:   B->assembled        = PETSC_FALSE;
3429:   if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
3430:   PetscFunctionReturn(PETSC_SUCCESS);
3431: }

3433: static PetscErrorCode MatSeqBAIJSetPreallocationCSR_SeqBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[])
3434: {
3435:   PetscInt     i, m, nz, nz_max = 0, *nnz;
3436:   PetscScalar *values      = NULL;
3437:   PetscBool    roworiented = ((Mat_SeqBAIJ *)B->data)->roworiented;

3439:   PetscFunctionBegin;
3440:   PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid block size specified, must be positive but it is %" PetscInt_FMT, bs);
3441:   PetscCall(PetscLayoutSetBlockSize(B->rmap, bs));
3442:   PetscCall(PetscLayoutSetBlockSize(B->cmap, bs));
3443:   PetscCall(PetscLayoutSetUp(B->rmap));
3444:   PetscCall(PetscLayoutSetUp(B->cmap));
3445:   PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
3446:   m = B->rmap->n / bs;

3448:   PetscCheck(ii[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]);
3449:   PetscCall(PetscMalloc1(m + 1, &nnz));
3450:   for (i = 0; i < m; i++) {
3451:     nz = ii[i + 1] - ii[i];
3452:     PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz);
3453:     nz_max = PetscMax(nz_max, nz);
3454:     nnz[i] = nz;
3455:   }
3456:   PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, nnz));
3457:   PetscCall(PetscFree(nnz));

3459:   values = (PetscScalar *)V;
3460:   if (!values) PetscCall(PetscCalloc1(bs * bs * (nz_max + 1), &values));
3461:   for (i = 0; i < m; i++) {
3462:     PetscInt        ncols = ii[i + 1] - ii[i];
3463:     const PetscInt *icols = jj + ii[i];
3464:     if (bs == 1 || !roworiented) {
3465:       const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0);
3466:       PetscCall(MatSetValuesBlocked_SeqBAIJ(B, 1, &i, ncols, icols, svals, INSERT_VALUES));
3467:     } else {
3468:       PetscInt j;
3469:       for (j = 0; j < ncols; j++) {
3470:         const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0);
3471:         PetscCall(MatSetValuesBlocked_SeqBAIJ(B, 1, &i, 1, &icols[j], svals, INSERT_VALUES));
3472:       }
3473:     }
3474:   }
3475:   if (!V) PetscCall(PetscFree(values));
3476:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
3477:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
3478:   PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
3479:   PetscFunctionReturn(PETSC_SUCCESS);
3480: }

3482: /*@C
3483:   MatSeqBAIJGetArray - gives read/write access to the array where the data for a `MATSEQBAIJ` matrix is stored

3485:   Not Collective

3487:   Input Parameter:
3488: . A - a `MATSEQBAIJ` matrix

3490:   Output Parameter:
3491: . array - pointer to the data

3493:   Level: intermediate

3495: .seealso: [](ch_matrices), `Mat`, `MATSEQBAIJ`, `MatSeqBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()`
3496: @*/
3497: PetscErrorCode MatSeqBAIJGetArray(Mat A, PetscScalar **array)
3498: {
3499:   PetscFunctionBegin;
3500:   PetscUseMethod(A, "MatSeqBAIJGetArray_C", (Mat, PetscScalar **), (A, array));
3501:   PetscFunctionReturn(PETSC_SUCCESS);
3502: }

3504: /*@C
3505:   MatSeqBAIJRestoreArray - returns access to the array where the data for a `MATSEQBAIJ` matrix is stored obtained by `MatSeqBAIJGetArray()`

3507:   Not Collective

3509:   Input Parameters:
3510: + A     - a `MATSEQBAIJ` matrix
3511: - array - pointer to the data

3513:   Level: intermediate

3515: .seealso: [](ch_matrices), `Mat`, `MatSeqBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()`
3516: @*/
3517: PetscErrorCode MatSeqBAIJRestoreArray(Mat A, PetscScalar **array)
3518: {
3519:   PetscFunctionBegin;
3520:   PetscUseMethod(A, "MatSeqBAIJRestoreArray_C", (Mat, PetscScalar **), (A, array));
3521:   PetscFunctionReturn(PETSC_SUCCESS);
3522: }

3524: /*MC
3525:    MATSEQBAIJ - MATSEQBAIJ = "seqbaij" - A matrix type to be used for sequential block sparse matrices, based on
3526:    block sparse compressed row format.

3528:    Options Database Keys:
3529: + -mat_type seqbaij - sets the matrix type to `MATSEQBAIJ` during a call to `MatSetFromOptions()`
3530: - -mat_baij_mult_version version - indicate the version of the matrix-vector product to use (0 often indicates using BLAS)

3532:    Level: beginner

3534:    Notes:
3535:    `MatSetOptions`(,`MAT_STRUCTURE_ONLY`,`PETSC_TRUE`) may be called for this matrix type. In this no
3536:    space is allocated for the nonzero entries and any entries passed with `MatSetValues()` are ignored

3538:    Run with `-info` to see what version of the matrix-vector product is being used

3540: .seealso: [](ch_matrices), `Mat`, `MatCreateSeqBAIJ()`
3541: M*/

3543: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqBSTRM(Mat, MatType, MatReuse, Mat *);

3545: PETSC_EXTERN PetscErrorCode MatCreate_SeqBAIJ(Mat B)
3546: {
3547:   PetscMPIInt  size;
3548:   Mat_SeqBAIJ *b;

3550:   PetscFunctionBegin;
3551:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
3552:   PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");

3554:   PetscCall(PetscNew(&b));
3555:   B->data   = (void *)b;
3556:   B->ops[0] = MatOps_Values;

3558:   b->row          = NULL;
3559:   b->col          = NULL;
3560:   b->icol         = NULL;
3561:   b->reallocs     = 0;
3562:   b->saved_values = NULL;

3564:   b->roworiented        = PETSC_TRUE;
3565:   b->nonew              = 0;
3566:   b->diag               = NULL;
3567:   B->spptr              = NULL;
3568:   B->info.nz_unneeded   = (PetscReal)b->maxnz * b->bs2;
3569:   b->keepnonzeropattern = PETSC_FALSE;

3571:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJGetArray_C", MatSeqBAIJGetArray_SeqBAIJ));
3572:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJRestoreArray_C", MatSeqBAIJRestoreArray_SeqBAIJ));
3573:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqBAIJ));
3574:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqBAIJ));
3575:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJSetColumnIndices_C", MatSeqBAIJSetColumnIndices_SeqBAIJ));
3576:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_seqaij_C", MatConvert_SeqBAIJ_SeqAIJ));
3577:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_seqsbaij_C", MatConvert_SeqBAIJ_SeqSBAIJ));
3578:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJSetPreallocation_C", MatSeqBAIJSetPreallocation_SeqBAIJ));
3579:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJSetPreallocationCSR_C", MatSeqBAIJSetPreallocationCSR_SeqBAIJ));
3580:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_SeqBAIJ));
3581: #if defined(PETSC_HAVE_HYPRE)
3582:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_hypre_C", MatConvert_AIJ_HYPRE));
3583: #endif
3584:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_is_C", MatConvert_XAIJ_IS));
3585:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQBAIJ));
3586:   PetscFunctionReturn(PETSC_SUCCESS);
3587: }

3589: PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat C, Mat A, MatDuplicateOption cpvalues, PetscBool mallocmatspace)
3590: {
3591:   Mat_SeqBAIJ *c = (Mat_SeqBAIJ *)C->data, *a = (Mat_SeqBAIJ *)A->data;
3592:   PetscInt     i, mbs = a->mbs, nz = a->nz, bs2 = a->bs2;

3594:   PetscFunctionBegin;
3595:   PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Cannot duplicate unassembled matrix");
3596:   PetscCheck(a->i[mbs] == nz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt matrix");

3598:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3599:     c->imax           = a->imax;
3600:     c->ilen           = a->ilen;
3601:     c->free_imax_ilen = PETSC_FALSE;
3602:   } else {
3603:     PetscCall(PetscMalloc2(mbs, &c->imax, mbs, &c->ilen));
3604:     for (i = 0; i < mbs; i++) {
3605:       c->imax[i] = a->imax[i];
3606:       c->ilen[i] = a->ilen[i];
3607:     }
3608:     c->free_imax_ilen = PETSC_TRUE;
3609:   }

3611:   /* allocate the matrix space */
3612:   if (mallocmatspace) {
3613:     if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3614:       PetscCall(PetscCalloc1(bs2 * nz, &c->a));

3616:       c->i            = a->i;
3617:       c->j            = a->j;
3618:       c->singlemalloc = PETSC_FALSE;
3619:       c->free_a       = PETSC_TRUE;
3620:       c->free_ij      = PETSC_FALSE;
3621:       c->parent       = A;
3622:       C->preallocated = PETSC_TRUE;
3623:       C->assembled    = PETSC_TRUE;

3625:       PetscCall(PetscObjectReference((PetscObject)A));
3626:       PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
3627:       PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
3628:     } else {
3629:       PetscCall(PetscMalloc3(bs2 * nz, &c->a, nz, &c->j, mbs + 1, &c->i));

3631:       c->singlemalloc = PETSC_TRUE;
3632:       c->free_a       = PETSC_TRUE;
3633:       c->free_ij      = PETSC_TRUE;

3635:       PetscCall(PetscArraycpy(c->i, a->i, mbs + 1));
3636:       if (mbs > 0) {
3637:         PetscCall(PetscArraycpy(c->j, a->j, nz));
3638:         if (cpvalues == MAT_COPY_VALUES) {
3639:           PetscCall(PetscArraycpy(c->a, a->a, bs2 * nz));
3640:         } else {
3641:           PetscCall(PetscArrayzero(c->a, bs2 * nz));
3642:         }
3643:       }
3644:       C->preallocated = PETSC_TRUE;
3645:       C->assembled    = PETSC_TRUE;
3646:     }
3647:   }

3649:   c->roworiented = a->roworiented;
3650:   c->nonew       = a->nonew;

3652:   PetscCall(PetscLayoutReference(A->rmap, &C->rmap));
3653:   PetscCall(PetscLayoutReference(A->cmap, &C->cmap));

3655:   c->bs2 = a->bs2;
3656:   c->mbs = a->mbs;
3657:   c->nbs = a->nbs;

3659:   if (a->diag) {
3660:     if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3661:       c->diag      = a->diag;
3662:       c->free_diag = PETSC_FALSE;
3663:     } else {
3664:       PetscCall(PetscMalloc1(mbs + 1, &c->diag));
3665:       for (i = 0; i < mbs; i++) c->diag[i] = a->diag[i];
3666:       c->free_diag = PETSC_TRUE;
3667:     }
3668:   } else c->diag = NULL;

3670:   c->nz         = a->nz;
3671:   c->maxnz      = a->nz; /* Since we allocate exactly the right amount */
3672:   c->solve_work = NULL;
3673:   c->mult_work  = NULL;
3674:   c->sor_workt  = NULL;
3675:   c->sor_work   = NULL;

3677:   c->compressedrow.use   = a->compressedrow.use;
3678:   c->compressedrow.nrows = a->compressedrow.nrows;
3679:   if (a->compressedrow.use) {
3680:     i = a->compressedrow.nrows;
3681:     PetscCall(PetscMalloc2(i + 1, &c->compressedrow.i, i + 1, &c->compressedrow.rindex));
3682:     PetscCall(PetscArraycpy(c->compressedrow.i, a->compressedrow.i, i + 1));
3683:     PetscCall(PetscArraycpy(c->compressedrow.rindex, a->compressedrow.rindex, i));
3684:   } else {
3685:     c->compressedrow.use    = PETSC_FALSE;
3686:     c->compressedrow.i      = NULL;
3687:     c->compressedrow.rindex = NULL;
3688:   }
3689:   c->nonzerorowcnt = a->nonzerorowcnt;
3690:   C->nonzerostate  = A->nonzerostate;

3692:   PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist));
3693:   PetscFunctionReturn(PETSC_SUCCESS);
3694: }

3696: PetscErrorCode MatDuplicate_SeqBAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B)
3697: {
3698:   PetscFunctionBegin;
3699:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
3700:   PetscCall(MatSetSizes(*B, A->rmap->N, A->cmap->n, A->rmap->N, A->cmap->n));
3701:   PetscCall(MatSetType(*B, MATSEQBAIJ));
3702:   PetscCall(MatDuplicateNoCreate_SeqBAIJ(*B, A, cpvalues, PETSC_TRUE));
3703:   PetscFunctionReturn(PETSC_SUCCESS);
3704: }

3706: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
3707: PetscErrorCode MatLoad_SeqBAIJ_Binary(Mat mat, PetscViewer viewer)
3708: {
3709:   PetscInt     header[4], M, N, nz, bs, m, n, mbs, nbs, rows, cols, sum, i, j, k;
3710:   PetscInt    *rowidxs, *colidxs;
3711:   PetscScalar *matvals;

3713:   PetscFunctionBegin;
3714:   PetscCall(PetscViewerSetUp(viewer));

3716:   /* read matrix header */
3717:   PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
3718:   PetscCheck(header[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
3719:   M  = header[1];
3720:   N  = header[2];
3721:   nz = header[3];
3722:   PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
3723:   PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
3724:   PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix stored in special format on disk, cannot load as SeqBAIJ");

3726:   /* set block sizes from the viewer's .info file */
3727:   PetscCall(MatLoad_Binary_BlockSizes(mat, viewer));
3728:   /* set local and global sizes if not set already */
3729:   if (mat->rmap->n < 0) mat->rmap->n = M;
3730:   if (mat->cmap->n < 0) mat->cmap->n = N;
3731:   if (mat->rmap->N < 0) mat->rmap->N = M;
3732:   if (mat->cmap->N < 0) mat->cmap->N = N;
3733:   PetscCall(PetscLayoutSetUp(mat->rmap));
3734:   PetscCall(PetscLayoutSetUp(mat->cmap));

3736:   /* check if the matrix sizes are correct */
3737:   PetscCall(MatGetSize(mat, &rows, &cols));
3738:   PetscCheck(M == rows && N == cols, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols);
3739:   PetscCall(MatGetBlockSize(mat, &bs));
3740:   PetscCall(MatGetLocalSize(mat, &m, &n));
3741:   mbs = m / bs;
3742:   nbs = n / bs;

3744:   /* read in row lengths, column indices and nonzero values */
3745:   PetscCall(PetscMalloc1(m + 1, &rowidxs));
3746:   PetscCall(PetscViewerBinaryRead(viewer, rowidxs + 1, m, NULL, PETSC_INT));
3747:   rowidxs[0] = 0;
3748:   for (i = 0; i < m; i++) rowidxs[i + 1] += rowidxs[i];
3749:   sum = rowidxs[m];
3750:   PetscCheck(sum == nz, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT, nz, sum);

3752:   /* read in column indices and nonzero values */
3753:   PetscCall(PetscMalloc2(rowidxs[m], &colidxs, nz, &matvals));
3754:   PetscCall(PetscViewerBinaryRead(viewer, colidxs, rowidxs[m], NULL, PETSC_INT));
3755:   PetscCall(PetscViewerBinaryRead(viewer, matvals, rowidxs[m], NULL, PETSC_SCALAR));

3757:   {               /* preallocate matrix storage */
3758:     PetscBT   bt; /* helper bit set to count nonzeros */
3759:     PetscInt *nnz;
3760:     PetscBool sbaij;

3762:     PetscCall(PetscBTCreate(nbs, &bt));
3763:     PetscCall(PetscCalloc1(mbs, &nnz));
3764:     PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQSBAIJ, &sbaij));
3765:     for (i = 0; i < mbs; i++) {
3766:       PetscCall(PetscBTMemzero(nbs, bt));
3767:       for (k = 0; k < bs; k++) {
3768:         PetscInt row = bs * i + k;
3769:         for (j = rowidxs[row]; j < rowidxs[row + 1]; j++) {
3770:           PetscInt col = colidxs[j];
3771:           if (!sbaij || col >= row)
3772:             if (!PetscBTLookupSet(bt, col / bs)) nnz[i]++;
3773:         }
3774:       }
3775:     }
3776:     PetscCall(PetscBTDestroy(&bt));
3777:     PetscCall(MatSeqBAIJSetPreallocation(mat, bs, 0, nnz));
3778:     PetscCall(MatSeqSBAIJSetPreallocation(mat, bs, 0, nnz));
3779:     PetscCall(PetscFree(nnz));
3780:   }

3782:   /* store matrix values */
3783:   for (i = 0; i < m; i++) {
3784:     PetscInt row = i, s = rowidxs[i], e = rowidxs[i + 1];
3785:     PetscUseTypeMethod(mat, setvalues, 1, &row, e - s, colidxs + s, matvals + s, INSERT_VALUES);
3786:   }

3788:   PetscCall(PetscFree(rowidxs));
3789:   PetscCall(PetscFree2(colidxs, matvals));
3790:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3791:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3792:   PetscFunctionReturn(PETSC_SUCCESS);
3793: }

3795: PetscErrorCode MatLoad_SeqBAIJ(Mat mat, PetscViewer viewer)
3796: {
3797:   PetscBool isbinary;

3799:   PetscFunctionBegin;
3800:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
3801:   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);
3802:   PetscCall(MatLoad_SeqBAIJ_Binary(mat, viewer));
3803:   PetscFunctionReturn(PETSC_SUCCESS);
3804: }

3806: /*@C
3807:   MatCreateSeqBAIJ - Creates a sparse matrix in `MATSEQAIJ` (block
3808:   compressed row) format.  For good matrix assembly performance the
3809:   user should preallocate the matrix storage by setting the parameter `nz`
3810:   (or the array `nnz`).

3812:   Collective

3814:   Input Parameters:
3815: + comm - MPI communicator, set to `PETSC_COMM_SELF`
3816: . bs   - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
3817:          blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
3818: . m    - number of rows
3819: . n    - number of columns
3820: . nz   - number of nonzero blocks  per block row (same for all rows)
3821: - nnz  - array containing the number of nonzero blocks in the various block rows
3822:          (possibly different for each block row) or `NULL`

3824:   Output Parameter:
3825: . A - the matrix

3827:   Options Database Keys:
3828: + -mat_no_unroll  - uses code that does not unroll the loops in the block calculations (much slower)
3829: - -mat_block_size - size of the blocks to use

3831:   Level: intermediate

3833:   Notes:
3834:   It is recommended that one use `MatCreateFromOptions()` or the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
3835:   MatXXXXSetPreallocation() paradigm instead of this routine directly.
3836:   [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`]

3838:   The number of rows and columns must be divisible by blocksize.

3840:   If the `nnz` parameter is given then the `nz` parameter is ignored

3842:   A nonzero block is any block that as 1 or more nonzeros in it

3844:   The `MATSEQBAIJ` format is fully compatible with standard Fortran
3845:   storage.  That is, the stored row and column indices can begin at
3846:   either one (as in Fortran) or zero.

3848:   Specify the preallocated storage with either `nz` or `nnz` (not both).
3849:   Set `nz` = `PETSC_DEFAULT` and `nnz` = `NULL` for PETSc to control dynamic memory
3850:   allocation.  See [Sparse Matrices](sec_matsparse) for details.
3851:   matrices.

3853: .seealso: [](ch_matrices), `Mat`, [Sparse Matrices](sec_matsparse), `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`
3854: @*/
3855: PetscErrorCode MatCreateSeqBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A)
3856: {
3857:   PetscFunctionBegin;
3858:   PetscCall(MatCreate(comm, A));
3859:   PetscCall(MatSetSizes(*A, m, n, m, n));
3860:   PetscCall(MatSetType(*A, MATSEQBAIJ));
3861:   PetscCall(MatSeqBAIJSetPreallocation(*A, bs, nz, (PetscInt *)nnz));
3862:   PetscFunctionReturn(PETSC_SUCCESS);
3863: }

3865: /*@C
3866:   MatSeqBAIJSetPreallocation - Sets the block size and expected nonzeros
3867:   per row in the matrix. For good matrix assembly performance the
3868:   user should preallocate the matrix storage by setting the parameter `nz`
3869:   (or the array `nnz`).

3871:   Collective

3873:   Input Parameters:
3874: + B   - the matrix
3875: . bs  - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
3876:         blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
3877: . nz  - number of block nonzeros per block row (same for all rows)
3878: - nnz - array containing the number of block nonzeros in the various block rows
3879:         (possibly different for each block row) or `NULL`

3881:   Options Database Keys:
3882: + -mat_no_unroll  - uses code that does not unroll the loops in the block calculations (much slower)
3883: - -mat_block_size - size of the blocks to use

3885:   Level: intermediate

3887:   Notes:
3888:   If the `nnz` parameter is given then the `nz` parameter is ignored

3890:   You can call `MatGetInfo()` to get information on how effective the preallocation was;
3891:   for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3892:   You can also run with the option `-info` and look for messages with the string
3893:   malloc in them to see if additional memory allocation was needed.

3895:   The `MATSEQBAIJ` format is fully compatible with standard Fortran
3896:   storage.  That is, the stored row and column indices can begin at
3897:   either one (as in Fortran) or zero.

3899:   Specify the preallocated storage with either `nz` or `nnz` (not both).
3900:   Set `nz` = `PETSC_DEFAULT` and `nnz` = `NULL` for PETSc to control dynamic memory
3901:   allocation.  See [Sparse Matrices](sec_matsparse) for details.

3903: .seealso: [](ch_matrices), `Mat`, [Sparse Matrices](sec_matsparse), `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`, `MatGetInfo()`
3904: @*/
3905: PetscErrorCode MatSeqBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[])
3906: {
3907:   PetscFunctionBegin;
3911:   PetscTryMethod(B, "MatSeqBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[]), (B, bs, nz, nnz));
3912:   PetscFunctionReturn(PETSC_SUCCESS);
3913: }

3915: /*@C
3916:   MatSeqBAIJSetPreallocationCSR - Creates a sparse sequential matrix in `MATSEQBAIJ` format using the given nonzero structure and (optional) numerical values

3918:   Collective

3920:   Input Parameters:
3921: + B  - the matrix
3922: . bs - the blocksize
3923: . i  - the indices into `j` for the start of each local row (indices start with zero)
3924: . j  - the column indices for each local row (indices start with zero) these must be sorted for each row
3925: - v  - optional values in the matrix, use `NULL` if not provided

3927:   Level: advanced

3929:   Notes:
3930:   The `i`,`j`,`v` values are COPIED with this routine; to avoid the copy use `MatCreateSeqBAIJWithArrays()`

3932:   The order of the entries in values is specified by the `MatOption` `MAT_ROW_ORIENTED`.  For example, C programs
3933:   may want to use the default `MAT_ROW_ORIENTED` of `PETSC_TRUE` and use an array v[nnz][bs][bs] where the second index is
3934:   over rows within a block and the last index is over columns within a block row.  Fortran programs will likely set
3935:   `MAT_ROW_ORIENTED` of `PETSC_FALSE` and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
3936:   block column and the second index is over columns within a block.

3938:   Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries and usually the numerical values as well

3940: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateSeqBAIJ()`, `MatSetValues()`, `MatSeqBAIJSetPreallocation()`, `MATSEQBAIJ`
3941: @*/
3942: PetscErrorCode MatSeqBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[])
3943: {
3944:   PetscFunctionBegin;
3948:   PetscTryMethod(B, "MatSeqBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v));
3949:   PetscFunctionReturn(PETSC_SUCCESS);
3950: }

3952: /*@
3953:   MatCreateSeqBAIJWithArrays - Creates a `MATSEQBAIJ` matrix using matrix elements provided by the user.

3955:   Collective

3957:   Input Parameters:
3958: + comm - must be an MPI communicator of size 1
3959: . bs   - size of block
3960: . m    - number of rows
3961: . n    - number of columns
3962: . i    - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row block row of the matrix
3963: . j    - column indices
3964: - a    - matrix values

3966:   Output Parameter:
3967: . mat - the matrix

3969:   Level: advanced

3971:   Notes:
3972:   The `i`, `j`, and `a` arrays are not copied by this routine, the user must free these arrays
3973:   once the matrix is destroyed

3975:   You cannot set new nonzero locations into this matrix, that will generate an error.

3977:   The `i` and `j` indices are 0 based

3979:   When block size is greater than 1 the matrix values must be stored using the `MATSEQBAIJ` storage format

3981:   The order of the entries in values is the same as the block compressed sparse row storage format; that is, it is
3982:   the same as a three dimensional array in Fortran values(bs,bs,nnz) that contains the first column of the first
3983:   block, followed by the second column of the first block etc etc.  That is, the blocks are contiguous in memory
3984:   with column-major ordering within blocks.

3986: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateBAIJ()`, `MatCreateSeqBAIJ()`
3987: @*/
3988: PetscErrorCode MatCreateSeqBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat)
3989: {
3990:   Mat_SeqBAIJ *baij;

3992:   PetscFunctionBegin;
3993:   PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "block size %" PetscInt_FMT " > 1 is not supported yet", bs);
3994:   if (m > 0) PetscCheck(i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");

3996:   PetscCall(MatCreate(comm, mat));
3997:   PetscCall(MatSetSizes(*mat, m, n, m, n));
3998:   PetscCall(MatSetType(*mat, MATSEQBAIJ));
3999:   PetscCall(MatSeqBAIJSetPreallocation(*mat, bs, MAT_SKIP_ALLOCATION, NULL));
4000:   baij = (Mat_SeqBAIJ *)(*mat)->data;
4001:   PetscCall(PetscMalloc2(m, &baij->imax, m, &baij->ilen));

4003:   baij->i = i;
4004:   baij->j = j;
4005:   baij->a = a;

4007:   baij->singlemalloc   = PETSC_FALSE;
4008:   baij->nonew          = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
4009:   baij->free_a         = PETSC_FALSE;
4010:   baij->free_ij        = PETSC_FALSE;
4011:   baij->free_imax_ilen = PETSC_TRUE;

4013:   for (PetscInt ii = 0; ii < m; ii++) {
4014:     const PetscInt row_len = i[ii + 1] - i[ii];

4016:     baij->ilen[ii] = baij->imax[ii] = row_len;
4017:     PetscCheck(row_len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT, ii, row_len);
4018:   }
4019:   if (PetscDefined(USE_DEBUG)) {
4020:     for (PetscInt ii = 0; ii < baij->i[m]; ii++) {
4021:       PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]);
4022:       PetscCheck(j[ii] <= n - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]);
4023:     }
4024:   }

4026:   PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
4027:   PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
4028:   PetscFunctionReturn(PETSC_SUCCESS);
4029: }

4031: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
4032: {
4033:   PetscFunctionBegin;
4034:   PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIBAIJ(comm, inmat, n, scall, outmat));
4035:   PetscFunctionReturn(PETSC_SUCCESS);
4036: }