Actual source code: sbaij.c


  2: /*
  3:     Defines the basic matrix operations for the SBAIJ (compressed row)
  4:   matrix storage format.
  5: */
  6: #include <../src/mat/impls/baij/seq/baij.h>
  7: #include <../src/mat/impls/sbaij/seq/sbaij.h>
  8: #include <petscblaslapack.h>

 10: #include <../src/mat/impls/sbaij/seq/relax.h>
 11: #define USESHORT
 12: #include <../src/mat/impls/sbaij/seq/relax.h>

 14: #if defined(PETSC_HAVE_ELEMENTAL)
 15: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat, MatType, MatReuse, Mat *);
 16: #endif
 17: #if defined(PETSC_HAVE_SCALAPACK)
 18: PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *);
 19: #endif
 20: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat, MatType, MatReuse, Mat *);

 22: /*
 23:      Checks for missing diagonals
 24: */
 25: PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A, PetscBool *missing, PetscInt *dd)
 26: {
 27:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
 28:   PetscInt     *diag, *ii = a->i, i;

 30:   MatMarkDiagonal_SeqSBAIJ(A);
 31:   *missing = PETSC_FALSE;
 32:   if (A->rmap->n > 0 && !ii) {
 33:     *missing = PETSC_TRUE;
 34:     if (dd) *dd = 0;
 35:     PetscInfo(A, "Matrix has no entries therefore is missing diagonal\n");
 36:   } else {
 37:     diag = a->diag;
 38:     for (i = 0; i < a->mbs; i++) {
 39:       if (diag[i] >= ii[i + 1]) {
 40:         *missing = PETSC_TRUE;
 41:         if (dd) *dd = i;
 42:         break;
 43:       }
 44:     }
 45:   }
 46:   return 0;
 47: }

 49: PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A)
 50: {
 51:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
 52:   PetscInt      i, j;

 54:   if (!a->diag) {
 55:     PetscMalloc1(a->mbs, &a->diag);
 56:     a->free_diag = PETSC_TRUE;
 57:   }
 58:   for (i = 0; i < a->mbs; i++) {
 59:     a->diag[i] = a->i[i + 1];
 60:     for (j = a->i[i]; j < a->i[i + 1]; j++) {
 61:       if (a->j[j] == i) {
 62:         a->diag[i] = j;
 63:         break;
 64:       }
 65:     }
 66:   }
 67:   return 0;
 68: }

 70: static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done)
 71: {
 72:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
 73:   PetscInt      i, j, n = a->mbs, nz = a->i[n], *tia, *tja, bs = A->rmap->bs, k, l, cnt;
 74:   PetscInt    **ia = (PetscInt **)inia, **ja = (PetscInt **)inja;

 76:   *nn = n;
 77:   if (!ia) return 0;
 78:   if (symmetric) {
 79:     MatToSymmetricIJ_SeqAIJ(n, a->i, a->j, PETSC_FALSE, 0, 0, &tia, &tja);
 80:     nz = tia[n];
 81:   } else {
 82:     tia = a->i;
 83:     tja = a->j;
 84:   }

 86:   if (!blockcompressed && bs > 1) {
 87:     (*nn) *= bs;
 88:     /* malloc & create the natural set of indices */
 89:     PetscMalloc1((n + 1) * bs, ia);
 90:     if (n) {
 91:       (*ia)[0] = oshift;
 92:       for (j = 1; j < bs; j++) (*ia)[j] = (tia[1] - tia[0]) * bs + (*ia)[j - 1];
 93:     }

 95:     for (i = 1; i < n; i++) {
 96:       (*ia)[i * bs] = (tia[i] - tia[i - 1]) * bs + (*ia)[i * bs - 1];
 97:       for (j = 1; j < bs; j++) (*ia)[i * bs + j] = (tia[i + 1] - tia[i]) * bs + (*ia)[i * bs + j - 1];
 98:     }
 99:     if (n) (*ia)[n * bs] = (tia[n] - tia[n - 1]) * bs + (*ia)[n * bs - 1];

101:     if (inja) {
102:       PetscMalloc1(nz * bs * bs, ja);
103:       cnt = 0;
104:       for (i = 0; i < n; i++) {
105:         for (j = 0; j < bs; j++) {
106:           for (k = tia[i]; k < tia[i + 1]; k++) {
107:             for (l = 0; l < bs; l++) (*ja)[cnt++] = bs * tja[k] + l;
108:           }
109:         }
110:       }
111:     }

113:     if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
114:       PetscFree(tia);
115:       PetscFree(tja);
116:     }
117:   } else if (oshift == 1) {
118:     if (symmetric) {
119:       nz = tia[A->rmap->n / bs];
120:       /*  add 1 to i and j indices */
121:       for (i = 0; i < A->rmap->n / bs + 1; i++) tia[i] = tia[i] + 1;
122:       *ia = tia;
123:       if (ja) {
124:         for (i = 0; i < nz; i++) tja[i] = tja[i] + 1;
125:         *ja = tja;
126:       }
127:     } else {
128:       nz = a->i[A->rmap->n / bs];
129:       /* malloc space and  add 1 to i and j indices */
130:       PetscMalloc1(A->rmap->n / bs + 1, ia);
131:       for (i = 0; i < A->rmap->n / bs + 1; i++) (*ia)[i] = a->i[i] + 1;
132:       if (ja) {
133:         PetscMalloc1(nz, ja);
134:         for (i = 0; i < nz; i++) (*ja)[i] = a->j[i] + 1;
135:       }
136:     }
137:   } else {
138:     *ia = tia;
139:     if (ja) *ja = tja;
140:   }
141:   return 0;
142: }

144: static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done)
145: {
146:   if (!ia) return 0;
147:   if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
148:     PetscFree(*ia);
149:     if (ja) PetscFree(*ja);
150:   }
151:   return 0;
152: }

154: PetscErrorCode MatDestroy_SeqSBAIJ(Mat A)
155: {
156:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;

158: #if defined(PETSC_USE_LOG)
159:   PetscLogObjectState((PetscObject)A, "Rows=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, A->rmap->N, a->nz);
160: #endif
161:   MatSeqXAIJFreeAIJ(A, &a->a, &a->j, &a->i);
162:   if (a->free_diag) PetscFree(a->diag);
163:   ISDestroy(&a->row);
164:   ISDestroy(&a->col);
165:   ISDestroy(&a->icol);
166:   PetscFree(a->idiag);
167:   PetscFree(a->inode.size);
168:   if (a->free_imax_ilen) PetscFree2(a->imax, a->ilen);
169:   PetscFree(a->solve_work);
170:   PetscFree(a->sor_work);
171:   PetscFree(a->solves_work);
172:   PetscFree(a->mult_work);
173:   PetscFree(a->saved_values);
174:   if (a->free_jshort) PetscFree(a->jshort);
175:   PetscFree(a->inew);
176:   MatDestroy(&a->parent);
177:   PetscFree(A->data);

179:   PetscObjectChangeTypeName((PetscObject)A, NULL);
180:   PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJGetArray_C", NULL);
181:   PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJRestoreArray_C", NULL);
182:   PetscObjectComposeFunction((PetscObject)A, "MatStoreValues_C", NULL);
183:   PetscObjectComposeFunction((PetscObject)A, "MatRetrieveValues_C", NULL);
184:   PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetColumnIndices_C", NULL);
185:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_seqaij_C", NULL);
186:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_seqbaij_C", NULL);
187:   PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetPreallocation_C", NULL);
188:   PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetPreallocationCSR_C", NULL);
189: #if defined(PETSC_HAVE_ELEMENTAL)
190:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_elemental_C", NULL);
191: #endif
192: #if defined(PETSC_HAVE_SCALAPACK)
193:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_scalapack_C", NULL);
194: #endif
195:   PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL);
196:   return 0;
197: }

199: PetscErrorCode MatSetOption_SeqSBAIJ(Mat A, MatOption op, PetscBool flg)
200: {
201:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
202: #if defined(PETSC_USE_COMPLEX)
203:   PetscInt bs;
204: #endif

206: #if defined(PETSC_USE_COMPLEX)
207:   MatGetBlockSize(A, &bs);
208: #endif
209:   switch (op) {
210:   case MAT_ROW_ORIENTED:
211:     a->roworiented = flg;
212:     break;
213:   case MAT_KEEP_NONZERO_PATTERN:
214:     a->keepnonzeropattern = flg;
215:     break;
216:   case MAT_NEW_NONZERO_LOCATIONS:
217:     a->nonew = (flg ? 0 : 1);
218:     break;
219:   case MAT_NEW_NONZERO_LOCATION_ERR:
220:     a->nonew = (flg ? -1 : 0);
221:     break;
222:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
223:     a->nonew = (flg ? -2 : 0);
224:     break;
225:   case MAT_UNUSED_NONZERO_LOCATION_ERR:
226:     a->nounused = (flg ? -1 : 0);
227:     break;
228:   case MAT_FORCE_DIAGONAL_ENTRIES:
229:   case MAT_IGNORE_OFF_PROC_ENTRIES:
230:   case MAT_USE_HASH_TABLE:
231:   case MAT_SORTED_FULL:
232:     PetscInfo(A, "Option %s ignored\n", MatOptions[op]);
233:     break;
234:   case MAT_HERMITIAN:
235: #if defined(PETSC_USE_COMPLEX)
236:     if (flg) { /* disable transpose ops */
238:       A->ops->multtranspose    = NULL;
239:       A->ops->multtransposeadd = NULL;
240:       A->symmetric             = PETSC_BOOL3_FALSE;
241:     }
242: #endif
243:     break;
244:   case MAT_SYMMETRIC:
245:   case MAT_SPD:
246: #if defined(PETSC_USE_COMPLEX)
247:     if (flg) { /* An hermitian and symmetric matrix has zero imaginary part (restore back transpose ops) */
248:       A->ops->multtranspose    = A->ops->mult;
249:       A->ops->multtransposeadd = A->ops->multadd;
250:     }
251: #endif
252:     break;
253:     /* These options are handled directly by MatSetOption() */
254:   case MAT_STRUCTURALLY_SYMMETRIC:
255:   case MAT_SYMMETRY_ETERNAL:
256:   case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
257:   case MAT_STRUCTURE_ONLY:
258:   case MAT_SPD_ETERNAL:
259:     /* These options are handled directly by MatSetOption() */
260:     break;
261:   case MAT_IGNORE_LOWER_TRIANGULAR:
262:     a->ignore_ltriangular = flg;
263:     break;
264:   case MAT_ERROR_LOWER_TRIANGULAR:
265:     a->ignore_ltriangular = flg;
266:     break;
267:   case MAT_GETROW_UPPERTRIANGULAR:
268:     a->getrow_utriangular = flg;
269:     break;
270:   case MAT_SUBMAT_SINGLEIS:
271:     break;
272:   default:
273:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %d", op);
274:   }
275:   return 0;
276: }

278: PetscErrorCode MatGetRow_SeqSBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
279: {
280:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;


284:   /* Get the upper triangular part of the row */
285:   MatGetRow_SeqBAIJ_private(A, row, nz, idx, v, a->i, a->j, a->a);
286:   return 0;
287: }

289: PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
290: {
291:   if (nz) *nz = 0;
292:   if (idx) PetscFree(*idx);
293:   if (v) PetscFree(*v);
294:   return 0;
295: }

297: PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A)
298: {
299:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;

301:   a->getrow_utriangular = PETSC_TRUE;
302:   return 0;
303: }

305: PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A)
306: {
307:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;

309:   a->getrow_utriangular = PETSC_FALSE;
310:   return 0;
311: }

313: PetscErrorCode MatTranspose_SeqSBAIJ(Mat A, MatReuse reuse, Mat *B)
314: {
315:   if (reuse == MAT_REUSE_MATRIX) MatTransposeCheckNonzeroState_Private(A, *B);
316:   if (reuse == MAT_INITIAL_MATRIX) {
317:     MatDuplicate(A, MAT_COPY_VALUES, B);
318:   } else if (reuse == MAT_REUSE_MATRIX) {
319:     MatCopy(A, *B, SAME_NONZERO_PATTERN);
320:   }
321:   return 0;
322: }

324: PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A, PetscViewer viewer)
325: {
326:   Mat_SeqSBAIJ     *a = (Mat_SeqSBAIJ *)A->data;
327:   PetscInt          i, j, bs = A->rmap->bs, k, l, bs2 = a->bs2;
328:   PetscViewerFormat format;
329:   PetscInt         *diag;

331:   PetscViewerGetFormat(viewer, &format);
332:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
333:     PetscViewerASCIIPrintf(viewer, "  block size is %" PetscInt_FMT "\n", bs);
334:   } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
335:     Mat         aij;
336:     const char *matname;

338:     if (A->factortype && bs > 1) {
339:       PetscPrintf(PETSC_COMM_SELF, "Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");
340:       return 0;
341:     }
342:     MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &aij);
343:     PetscObjectGetName((PetscObject)A, &matname);
344:     PetscObjectSetName((PetscObject)aij, matname);
345:     MatView(aij, viewer);
346:     MatDestroy(&aij);
347:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
348:     PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
349:     for (i = 0; i < a->mbs; i++) {
350:       for (j = 0; j < bs; j++) {
351:         PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j);
352:         for (k = a->i[i]; k < a->i[i + 1]; k++) {
353:           for (l = 0; l < bs; l++) {
354: #if defined(PETSC_USE_COMPLEX)
355:             if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0 && PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
356:               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]));
357:             } else if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0 && PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
358:               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]));
359:             } else if (PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
360:               PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]));
361:             }
362: #else
363:             if (a->a[bs2 * k + l * bs + j] != 0.0) PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)a->a[bs2 * k + l * bs + j]);
364: #endif
365:           }
366:         }
367:         PetscViewerASCIIPrintf(viewer, "\n");
368:       }
369:     }
370:     PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
371:   } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
372:     return 0;
373:   } else {
374:     PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
375:     if (A->factortype) { /* for factored matrix */

378:       diag = a->diag;
379:       for (i = 0; i < a->mbs; i++) { /* for row block i */
380:         PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i);
381:         /* diagonal entry */
382: #if defined(PETSC_USE_COMPLEX)
383:         if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) {
384:           PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]]), (double)PetscImaginaryPart(1.0 / a->a[diag[i]]));
385:         } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) {
386:           PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]]), -(double)PetscImaginaryPart(1.0 / a->a[diag[i]]));
387:         } else {
388:           PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]]));
389:         }
390: #else
391:         PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[diag[i]], (double)(1.0 / a->a[diag[i]]));
392: #endif
393:         /* off-diagonal entries */
394:         for (k = a->i[i]; k < a->i[i + 1] - 1; k++) {
395: #if defined(PETSC_USE_COMPLEX)
396:           if (PetscImaginaryPart(a->a[k]) > 0.0) {
397:             PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", bs * a->j[k], (double)PetscRealPart(a->a[k]), (double)PetscImaginaryPart(a->a[k]));
398:           } else if (PetscImaginaryPart(a->a[k]) < 0.0) {
399:             PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", bs * a->j[k], (double)PetscRealPart(a->a[k]), -(double)PetscImaginaryPart(a->a[k]));
400:           } else {
401:             PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k], (double)PetscRealPart(a->a[k]));
402:           }
403: #else
404:           PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[k], (double)a->a[k]);
405: #endif
406:         }
407:         PetscViewerASCIIPrintf(viewer, "\n");
408:       }

410:     } else {                         /* for non-factored matrix */
411:       for (i = 0; i < a->mbs; i++) { /* for row block i */
412:         for (j = 0; j < bs; j++) {   /* for row bs*i + j */
413:           PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j);
414:           for (k = a->i[i]; k < a->i[i + 1]; k++) { /* for column block */
415:             for (l = 0; l < bs; l++) {              /* for column */
416: #if defined(PETSC_USE_COMPLEX)
417:               if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0) {
418:                 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]));
419:               } else if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0) {
420:                 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]));
421:               } else {
422:                 PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]));
423:               }
424: #else
425:               PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)a->a[bs2 * k + l * bs + j]);
426: #endif
427:             }
428:           }
429:           PetscViewerASCIIPrintf(viewer, "\n");
430:         }
431:       }
432:     }
433:     PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
434:   }
435:   PetscViewerFlush(viewer);
436:   return 0;
437: }

439: #include <petscdraw.h>
440: static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw, void *Aa)
441: {
442:   Mat           A = (Mat)Aa;
443:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
444:   PetscInt      row, i, j, k, l, mbs = a->mbs, color, bs = A->rmap->bs, bs2 = a->bs2;
445:   PetscReal     xl, yl, xr, yr, x_l, x_r, y_l, y_r;
446:   MatScalar    *aa;
447:   PetscViewer   viewer;

449:   PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer);
450:   PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr);

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

454:   PetscDrawCollectiveBegin(draw);
455:   PetscDrawString(draw, .3 * (xl + xr), .3 * (yl + yr), PETSC_DRAW_BLACK, "symmetric");
456:   /* Blue for negative, Cyan for zero and  Red for positive */
457:   color = PETSC_DRAW_BLUE;
458:   for (i = 0, row = 0; i < mbs; i++, row += bs) {
459:     for (j = a->i[i]; j < a->i[i + 1]; j++) {
460:       y_l = A->rmap->N - row - 1.0;
461:       y_r = y_l + 1.0;
462:       x_l = a->j[j] * bs;
463:       x_r = x_l + 1.0;
464:       aa  = a->a + j * bs2;
465:       for (k = 0; k < bs; k++) {
466:         for (l = 0; l < bs; l++) {
467:           if (PetscRealPart(*aa++) >= 0.) continue;
468:           PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color);
469:         }
470:       }
471:     }
472:   }
473:   color = PETSC_DRAW_CYAN;
474:   for (i = 0, row = 0; i < mbs; i++, row += bs) {
475:     for (j = a->i[i]; j < a->i[i + 1]; j++) {
476:       y_l = A->rmap->N - row - 1.0;
477:       y_r = y_l + 1.0;
478:       x_l = a->j[j] * bs;
479:       x_r = x_l + 1.0;
480:       aa  = a->a + j * bs2;
481:       for (k = 0; k < bs; k++) {
482:         for (l = 0; l < bs; l++) {
483:           if (PetscRealPart(*aa++) != 0.) continue;
484:           PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color);
485:         }
486:       }
487:     }
488:   }
489:   color = PETSC_DRAW_RED;
490:   for (i = 0, row = 0; i < mbs; i++, row += bs) {
491:     for (j = a->i[i]; j < a->i[i + 1]; j++) {
492:       y_l = A->rmap->N - row - 1.0;
493:       y_r = y_l + 1.0;
494:       x_l = a->j[j] * bs;
495:       x_r = x_l + 1.0;
496:       aa  = a->a + j * bs2;
497:       for (k = 0; k < bs; k++) {
498:         for (l = 0; l < bs; l++) {
499:           if (PetscRealPart(*aa++) <= 0.) continue;
500:           PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color);
501:         }
502:       }
503:     }
504:   }
505:   PetscDrawCollectiveEnd(draw);
506:   return 0;
507: }

509: static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A, PetscViewer viewer)
510: {
511:   PetscReal xl, yl, xr, yr, w, h;
512:   PetscDraw draw;
513:   PetscBool isnull;

515:   PetscViewerDrawGetDraw(viewer, 0, &draw);
516:   PetscDrawIsNull(draw, &isnull);
517:   if (isnull) return 0;

519:   xr = A->rmap->N;
520:   yr = A->rmap->N;
521:   h  = yr / 10.0;
522:   w  = xr / 10.0;
523:   xr += w;
524:   yr += h;
525:   xl = -w;
526:   yl = -h;
527:   PetscDrawSetCoordinates(draw, xl, yl, xr, yr);
528:   PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer);
529:   PetscDrawZoom(draw, MatView_SeqSBAIJ_Draw_Zoom, A);
530:   PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL);
531:   PetscDrawSave(draw);
532:   return 0;
533: }

535: /* Used for both MPIBAIJ and MPISBAIJ matrices */
536: #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary

538: PetscErrorCode MatView_SeqSBAIJ(Mat A, PetscViewer viewer)
539: {
540:   PetscBool iascii, isbinary, isdraw;

542:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
543:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
544:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
545:   if (iascii) {
546:     MatView_SeqSBAIJ_ASCII(A, viewer);
547:   } else if (isbinary) {
548:     MatView_SeqSBAIJ_Binary(A, viewer);
549:   } else if (isdraw) {
550:     MatView_SeqSBAIJ_Draw(A, viewer);
551:   } else {
552:     Mat         B;
553:     const char *matname;
554:     MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &B);
555:     PetscObjectGetName((PetscObject)A, &matname);
556:     PetscObjectSetName((PetscObject)B, matname);
557:     MatView(B, viewer);
558:     MatDestroy(&B);
559:   }
560:   return 0;
561: }

563: PetscErrorCode MatGetValues_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], PetscScalar v[])
564: {
565:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
566:   PetscInt     *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j;
567:   PetscInt     *ai = a->i, *ailen = a->ilen;
568:   PetscInt      brow, bcol, ridx, cidx, bs = A->rmap->bs, bs2 = a->bs2;
569:   MatScalar    *ap, *aa = a->a;

571:   for (k = 0; k < m; k++) { /* loop over rows */
572:     row  = im[k];
573:     brow = row / bs;
574:     if (row < 0) {
575:       v += n;
576:       continue;
577:     } /* negative row */
579:     rp   = aj + ai[brow];
580:     ap   = aa + bs2 * ai[brow];
581:     nrow = ailen[brow];
582:     for (l = 0; l < n; l++) { /* loop over columns */
583:       if (in[l] < 0) {
584:         v++;
585:         continue;
586:       } /* negative column */
588:       col  = in[l];
589:       bcol = col / bs;
590:       cidx = col % bs;
591:       ridx = row % bs;
592:       high = nrow;
593:       low  = 0; /* assume unsorted */
594:       while (high - low > 5) {
595:         t = (low + high) / 2;
596:         if (rp[t] > bcol) high = t;
597:         else low = t;
598:       }
599:       for (i = low; i < high; i++) {
600:         if (rp[i] > bcol) break;
601:         if (rp[i] == bcol) {
602:           *v++ = ap[bs2 * i + bs * cidx + ridx];
603:           goto finished;
604:         }
605:       }
606:       *v++ = 0.0;
607:     finished:;
608:     }
609:   }
610:   return 0;
611: }

613: PetscErrorCode MatPermute_SeqSBAIJ(Mat A, IS rowp, IS colp, Mat *B)
614: {
615:   Mat C;

617:   MatConvert(A, MATSEQBAIJ, MAT_INITIAL_MATRIX, &C);
618:   MatPermute(C, rowp, colp, B);
619:   MatDestroy(&C);
620:   if (rowp == colp) MatConvert(*B, MATSEQSBAIJ, MAT_INPLACE_MATRIX, B);
621:   return 0;
622: }

624: PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
625: {
626:   Mat_SeqSBAIJ      *a = (Mat_SeqSBAIJ *)A->data;
627:   PetscInt          *rp, k, low, high, t, ii, jj, row, nrow, i, col, l, rmax, N, lastcol = -1;
628:   PetscInt          *imax = a->imax, *ai = a->i, *ailen = a->ilen;
629:   PetscInt          *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs, stepval;
630:   PetscBool          roworiented = a->roworiented;
631:   const PetscScalar *value       = v;
632:   MatScalar         *ap, *aa = a->a, *bap;

634:   if (roworiented) stepval = (n - 1) * bs;
635:   else stepval = (m - 1) * bs;

637:   for (k = 0; k < m; k++) { /* loop over added rows */
638:     row = im[k];
639:     if (row < 0) continue;
641:     rp   = aj + ai[row];
642:     ap   = aa + bs2 * ai[row];
643:     rmax = imax[row];
644:     nrow = ailen[row];
645:     low  = 0;
646:     high = nrow;
647:     for (l = 0; l < n; l++) { /* loop over added columns */
648:       if (in[l] < 0) continue;
649:       col = in[l];
651:       if (col < row) {
652:         if (a->ignore_ltriangular) continue; /* ignore lower triangular block */
653:         else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
654:       }
655:       if (roworiented) value = v + k * (stepval + bs) * bs + l * bs;
656:       else value = v + l * (stepval + bs) * bs + k * bs;

658:       if (col <= lastcol) low = 0;
659:       else high = nrow;

661:       lastcol = col;
662:       while (high - low > 7) {
663:         t = (low + high) / 2;
664:         if (rp[t] > col) high = t;
665:         else low = t;
666:       }
667:       for (i = low; i < high; i++) {
668:         if (rp[i] > col) break;
669:         if (rp[i] == col) {
670:           bap = ap + bs2 * i;
671:           if (roworiented) {
672:             if (is == ADD_VALUES) {
673:               for (ii = 0; ii < bs; ii++, value += stepval) {
674:                 for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
675:               }
676:             } else {
677:               for (ii = 0; ii < bs; ii++, value += stepval) {
678:                 for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
679:               }
680:             }
681:           } else {
682:             if (is == ADD_VALUES) {
683:               for (ii = 0; ii < bs; ii++, value += stepval) {
684:                 for (jj = 0; jj < bs; jj++) *bap++ += *value++;
685:               }
686:             } else {
687:               for (ii = 0; ii < bs; ii++, value += stepval) {
688:                 for (jj = 0; jj < bs; jj++) *bap++ = *value++;
689:               }
690:             }
691:           }
692:           goto noinsert2;
693:         }
694:       }
695:       if (nonew == 1) goto noinsert2;
697:       MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
698:       N = nrow++ - 1;
699:       high++;
700:       /* shift up all the later entries in this row */
701:       PetscArraymove(rp + i + 1, rp + i, N - i + 1);
702:       PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1));
703:       PetscArrayzero(ap + bs2 * i, bs2);
704:       rp[i] = col;
705:       bap   = ap + bs2 * i;
706:       if (roworiented) {
707:         for (ii = 0; ii < bs; ii++, value += stepval) {
708:           for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
709:         }
710:       } else {
711:         for (ii = 0; ii < bs; ii++, value += stepval) {
712:           for (jj = 0; jj < bs; jj++) *bap++ = *value++;
713:         }
714:       }
715:     noinsert2:;
716:       low = i;
717:     }
718:     ailen[row] = nrow;
719:   }
720:   return 0;
721: }

723: /*
724:     This is not yet used
725: */
726: PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A)
727: {
728:   Mat_SeqSBAIJ   *a  = (Mat_SeqSBAIJ *)A->data;
729:   const PetscInt *ai = a->i, *aj = a->j, *cols;
730:   PetscInt        i = 0, j, blk_size, m = A->rmap->n, node_count = 0, nzx, nzy, *ns, row, nz, cnt, cnt2, *counts;
731:   PetscBool       flag;

733:   PetscMalloc1(m, &ns);
734:   while (i < m) {
735:     nzx = ai[i + 1] - ai[i]; /* Number of nonzeros */
736:     /* Limits the number of elements in a node to 'a->inode.limit' */
737:     for (j = i + 1, blk_size = 1; j < m && blk_size < a->inode.limit; ++j, ++blk_size) {
738:       nzy = ai[j + 1] - ai[j];
739:       if (nzy != (nzx - j + i)) break;
740:       PetscArraycmp(aj + ai[i] + j - i, aj + ai[j], nzy, &flag);
741:       if (!flag) break;
742:     }
743:     ns[node_count++] = blk_size;

745:     i = j;
746:   }
747:   if (!a->inode.size && m && node_count > .9 * m) {
748:     PetscFree(ns);
749:     PetscInfo(A, "Found %" PetscInt_FMT " nodes out of %" PetscInt_FMT " rows. Not using Inode routines\n", node_count, m);
750:   } else {
751:     a->inode.node_count = node_count;

753:     PetscMalloc1(node_count, &a->inode.size);
754:     PetscArraycpy(a->inode.size, ns, node_count);
755:     PetscFree(ns);
756:     PetscInfo(A, "Found %" PetscInt_FMT " nodes of %" PetscInt_FMT ". Limit used: %" PetscInt_FMT ". Using Inode routines\n", node_count, m, a->inode.limit);

758:     /* count collections of adjacent columns in each inode */
759:     row = 0;
760:     cnt = 0;
761:     for (i = 0; i < node_count; i++) {
762:       cols = aj + ai[row] + a->inode.size[i];
763:       nz   = ai[row + 1] - ai[row] - a->inode.size[i];
764:       for (j = 1; j < nz; j++) {
765:         if (cols[j] != cols[j - 1] + 1) cnt++;
766:       }
767:       cnt++;
768:       row += a->inode.size[i];
769:     }
770:     PetscMalloc1(2 * cnt, &counts);
771:     cnt = 0;
772:     row = 0;
773:     for (i = 0; i < node_count; i++) {
774:       cols            = aj + ai[row] + a->inode.size[i];
775:       counts[2 * cnt] = cols[0];
776:       nz              = ai[row + 1] - ai[row] - a->inode.size[i];
777:       cnt2            = 1;
778:       for (j = 1; j < nz; j++) {
779:         if (cols[j] != cols[j - 1] + 1) {
780:           counts[2 * (cnt++) + 1] = cnt2;
781:           counts[2 * cnt]         = cols[j];
782:           cnt2                    = 1;
783:         } else cnt2++;
784:       }
785:       counts[2 * (cnt++) + 1] = cnt2;
786:       row += a->inode.size[i];
787:     }
788:     PetscIntView(2 * cnt, counts, NULL);
789:   }
790:   return 0;
791: }

793: PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A, MatAssemblyType mode)
794: {
795:   Mat_SeqSBAIJ *a      = (Mat_SeqSBAIJ *)A->data;
796:   PetscInt      fshift = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax;
797:   PetscInt      m = A->rmap->N, *ip, N, *ailen = a->ilen;
798:   PetscInt      mbs = a->mbs, bs2 = a->bs2, rmax = 0;
799:   MatScalar    *aa = a->a, *ap;

801:   if (mode == MAT_FLUSH_ASSEMBLY) return 0;

803:   if (m) rmax = ailen[0];
804:   for (i = 1; i < mbs; i++) {
805:     /* move each row back by the amount of empty slots (fshift) before it*/
806:     fshift += imax[i - 1] - ailen[i - 1];
807:     rmax = PetscMax(rmax, ailen[i]);
808:     if (fshift) {
809:       ip = aj + ai[i];
810:       ap = aa + bs2 * ai[i];
811:       N  = ailen[i];
812:       PetscArraymove(ip - fshift, ip, N);
813:       PetscArraymove(ap - bs2 * fshift, ap, bs2 * N);
814:     }
815:     ai[i] = ai[i - 1] + ailen[i - 1];
816:   }
817:   if (mbs) {
818:     fshift += imax[mbs - 1] - ailen[mbs - 1];
819:     ai[mbs] = ai[mbs - 1] + ailen[mbs - 1];
820:   }
821:   /* reset ilen and imax for each row */
822:   for (i = 0; i < mbs; i++) ailen[i] = imax[i] = ai[i + 1] - ai[i];
823:   a->nz = ai[mbs];

825:   /* diagonals may have moved, reset it */
826:   if (a->diag) PetscArraycpy(a->diag, ai, mbs);

829:   PetscInfo(A, "Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT ", block size %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded, %" PetscInt_FMT " used\n", m, A->rmap->N, A->rmap->bs, fshift * bs2, a->nz * bs2);
830:   PetscInfo(A, "Number of mallocs during MatSetValues is %" PetscInt_FMT "\n", a->reallocs);
831:   PetscInfo(A, "Most nonzeros blocks in any row is %" PetscInt_FMT "\n", rmax);

833:   A->info.mallocs += a->reallocs;
834:   a->reallocs         = 0;
835:   A->info.nz_unneeded = (PetscReal)fshift * bs2;
836:   a->idiagvalid       = PETSC_FALSE;
837:   a->rmax             = rmax;

839:   if (A->cmap->n < 65536 && A->cmap->bs == 1) {
840:     if (a->jshort && a->free_jshort) {
841:       /* when matrix data structure is changed, previous jshort must be replaced */
842:       PetscFree(a->jshort);
843:     }
844:     PetscMalloc1(a->i[A->rmap->n], &a->jshort);
845:     for (i = 0; i < a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i];
846:     A->ops->mult   = MatMult_SeqSBAIJ_1_ushort;
847:     A->ops->sor    = MatSOR_SeqSBAIJ_ushort;
848:     a->free_jshort = PETSC_TRUE;
849:   }
850:   return 0;
851: }

853: /*
854:    This function returns an array of flags which indicate the locations of contiguous
855:    blocks that should be zeroed. for eg: if bs = 3  and is = [0,1,2,3,5,6,7,8,9]
856:    then the resulting sizes = [3,1,1,3,1] corresponding to sets [(0,1,2),(3),(5),(6,7,8),(9)]
857:    Assume: sizes should be long enough to hold all the values.
858: */
859: PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[], PetscInt n, PetscInt bs, PetscInt sizes[], PetscInt *bs_max)
860: {
861:   PetscInt  i, j, k, row;
862:   PetscBool flg;

864:   for (i = 0, j = 0; i < n; j++) {
865:     row = idx[i];
866:     if (row % bs != 0) { /* Not the beginning of a block */
867:       sizes[j] = 1;
868:       i++;
869:     } else if (i + bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
870:       sizes[j] = 1;          /* Also makes sure at least 'bs' values exist for next else */
871:       i++;
872:     } else { /* Beginning of the block, so check if the complete block exists */
873:       flg = PETSC_TRUE;
874:       for (k = 1; k < bs; k++) {
875:         if (row + k != idx[i + k]) { /* break in the block */
876:           flg = PETSC_FALSE;
877:           break;
878:         }
879:       }
880:       if (flg) { /* No break in the bs */
881:         sizes[j] = bs;
882:         i += bs;
883:       } else {
884:         sizes[j] = 1;
885:         i++;
886:       }
887:     }
888:   }
889:   *bs_max = j;
890:   return 0;
891: }

893: /* Only add/insert a(i,j) with i<=j (blocks).
894:    Any a(i,j) with i>j input by user is ingored.
895: */

897: PetscErrorCode MatSetValues_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
898: {
899:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
900:   PetscInt     *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N, lastcol = -1;
901:   PetscInt     *imax = a->imax, *ai = a->i, *ailen = a->ilen, roworiented = a->roworiented;
902:   PetscInt     *aj = a->j, nonew = a->nonew, bs = A->rmap->bs, brow, bcol;
903:   PetscInt      ridx, cidx, bs2                 = a->bs2;
904:   MatScalar    *ap, value, *aa                  = a->a, *bap;

906:   for (k = 0; k < m; k++) { /* loop over added rows */
907:     row  = im[k];           /* row number */
908:     brow = row / bs;        /* block row number */
909:     if (row < 0) continue;
911:     rp   = aj + ai[brow];       /*ptr to beginning of column value of the row block*/
912:     ap   = aa + bs2 * ai[brow]; /*ptr to beginning of element value of the row block*/
913:     rmax = imax[brow];          /* maximum space allocated for this row */
914:     nrow = ailen[brow];         /* actual length of this row */
915:     low  = 0;
916:     high = nrow;
917:     for (l = 0; l < n; l++) { /* loop over added columns */
918:       if (in[l] < 0) continue;
920:       col  = in[l];
921:       bcol = col / bs; /* block col number */

923:       if (brow > bcol) {
924:         if (a->ignore_ltriangular) continue; /* ignore lower triangular values */
925:         else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
926:       }

928:       ridx = row % bs;
929:       cidx = col % bs; /*row and col index inside the block */
930:       if ((brow == bcol && ridx <= cidx) || (brow < bcol)) {
931:         /* element value a(k,l) */
932:         if (roworiented) value = v[l + k * n];
933:         else value = v[k + l * m];

935:         /* move pointer bap to a(k,l) quickly and add/insert value */
936:         if (col <= lastcol) low = 0;
937:         else high = nrow;

939:         lastcol = col;
940:         while (high - low > 7) {
941:           t = (low + high) / 2;
942:           if (rp[t] > bcol) high = t;
943:           else low = t;
944:         }
945:         for (i = low; i < high; i++) {
946:           if (rp[i] > bcol) break;
947:           if (rp[i] == bcol) {
948:             bap = ap + bs2 * i + bs * cidx + ridx;
949:             if (is == ADD_VALUES) *bap += value;
950:             else *bap = value;
951:             /* for diag block, add/insert its symmetric element a(cidx,ridx) */
952:             if (brow == bcol && ridx < cidx) {
953:               bap = ap + bs2 * i + bs * ridx + cidx;
954:               if (is == ADD_VALUES) *bap += value;
955:               else *bap = value;
956:             }
957:             goto noinsert1;
958:           }
959:         }

961:         if (nonew == 1) goto noinsert1;
963:         MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);

965:         N = nrow++ - 1;
966:         high++;
967:         /* shift up all the later entries in this row */
968:         PetscArraymove(rp + i + 1, rp + i, N - i + 1);
969:         PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1));
970:         PetscArrayzero(ap + bs2 * i, bs2);
971:         rp[i]                          = bcol;
972:         ap[bs2 * i + bs * cidx + ridx] = value;
973:         /* for diag block, add/insert its symmetric element a(cidx,ridx) */
974:         if (brow == bcol && ridx < cidx) ap[bs2 * i + bs * ridx + cidx] = value;
975:         A->nonzerostate++;
976:       noinsert1:;
977:         low = i;
978:       }
979:     } /* end of loop over added columns */
980:     ailen[brow] = nrow;
981:   } /* end of loop over added rows */
982:   return 0;
983: }

985: PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA, IS row, const MatFactorInfo *info)
986: {
987:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)inA->data;
988:   Mat           outA;
989:   PetscBool     row_identity;

992:   ISIdentity(row, &row_identity);

996:   outA            = inA;
997:   inA->factortype = MAT_FACTOR_ICC;
998:   PetscFree(inA->solvertype);
999:   PetscStrallocpy(MATSOLVERPETSC, &inA->solvertype);

1001:   MatMarkDiagonal_SeqSBAIJ(inA);
1002:   MatSeqSBAIJSetNumericFactorization_inplace(inA, row_identity);

1004:   PetscObjectReference((PetscObject)row);
1005:   ISDestroy(&a->row);
1006:   a->row = row;
1007:   PetscObjectReference((PetscObject)row);
1008:   ISDestroy(&a->col);
1009:   a->col = row;

1011:   /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
1012:   if (a->icol) ISInvertPermutation(row, PETSC_DECIDE, &a->icol);

1014:   if (!a->solve_work) { PetscMalloc1(inA->rmap->N + inA->rmap->bs, &a->solve_work); }

1016:   MatCholeskyFactorNumeric(outA, inA, info);
1017:   return 0;
1018: }

1020: PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat, PetscInt *indices)
1021: {
1022:   Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data;
1023:   PetscInt      i, nz, n;

1025:   nz = baij->maxnz;
1026:   n  = mat->cmap->n;
1027:   for (i = 0; i < nz; i++) baij->j[i] = indices[i];

1029:   baij->nz = nz;
1030:   for (i = 0; i < n; i++) baij->ilen[i] = baij->imax[i];

1032:   MatSetOption(mat, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
1033:   return 0;
1034: }

1036: /*@
1037:   MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1038:   in a `MATSEQSBAIJ` matrix.

1040:   Input Parameters:
1041:   +  mat     - the `MATSEQSBAIJ` matrix
1042:   -  indices - the column indices

1044:   Level: advanced

1046:   Notes:
1047:   This can be called if you have precomputed the nonzero structure of the
1048:   matrix and want to provide it to the matrix object to improve the performance
1049:   of the `MatSetValues()` operation.

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

1054:   MUST be called before any calls to MatSetValues()

1056:  .seealso: `MATSEQSBAIJ`, `MatCreateSeqSBAIJ`
1057: @*/
1058: PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat, PetscInt *indices)
1059: {
1062:   PetscUseMethod(mat, "MatSeqSBAIJSetColumnIndices_C", (Mat, PetscInt *), (mat, indices));
1063:   return 0;
1064: }

1066: PetscErrorCode MatCopy_SeqSBAIJ(Mat A, Mat B, MatStructure str)
1067: {
1068:   PetscBool isbaij;

1070:   PetscObjectTypeCompareAny((PetscObject)B, &isbaij, MATSEQSBAIJ, MATMPISBAIJ, "");
1072:   /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */
1073:   if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1074:     Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1075:     Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data;

1080:     PetscArraycpy(b->a, a->a, a->bs2 * a->i[a->mbs]);
1081:     PetscObjectStateIncrease((PetscObject)B);
1082:   } else {
1083:     MatGetRowUpperTriangular(A);
1084:     MatCopy_Basic(A, B, str);
1085:     MatRestoreRowUpperTriangular(A);
1086:   }
1087:   return 0;
1088: }

1090: PetscErrorCode MatSetUp_SeqSBAIJ(Mat A)
1091: {
1092:   MatSeqSBAIJSetPreallocation(A, A->rmap->bs, PETSC_DEFAULT, NULL);
1093:   return 0;
1094: }

1096: static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A, PetscScalar *array[])
1097: {
1098:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;

1100:   *array = a->a;
1101:   return 0;
1102: }

1104: static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A, PetscScalar *array[])
1105: {
1106:   *array = NULL;
1107:   return 0;
1108: }

1110: PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y, Mat X, PetscInt *nnz)
1111: {
1112:   PetscInt      bs = Y->rmap->bs, mbs = Y->rmap->N / bs;
1113:   Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ *)X->data;
1114:   Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ *)Y->data;

1116:   /* Set the number of nonzeros in the new matrix */
1117:   MatAXPYGetPreallocation_SeqX_private(mbs, x->i, x->j, y->i, y->j, nnz);
1118:   return 0;
1119: }

1121: PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str)
1122: {
1123:   Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ *)X->data, *y = (Mat_SeqSBAIJ *)Y->data;
1124:   PetscInt      bs = Y->rmap->bs, bs2 = bs * bs;
1125:   PetscBLASInt  one = 1;

1127:   if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) {
1128:     PetscBool e = x->nz == y->nz && x->mbs == y->mbs ? PETSC_TRUE : PETSC_FALSE;
1129:     if (e) {
1130:       PetscArraycmp(x->i, y->i, x->mbs + 1, &e);
1131:       if (e) {
1132:         PetscArraycmp(x->j, y->j, x->i[x->mbs], &e);
1133:         if (e) str = SAME_NONZERO_PATTERN;
1134:       }
1135:     }
1137:   }
1138:   if (str == SAME_NONZERO_PATTERN) {
1139:     PetscScalar  alpha = a;
1140:     PetscBLASInt bnz;
1141:     PetscBLASIntCast(x->nz * bs2, &bnz);
1142:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one));
1143:     PetscObjectStateIncrease((PetscObject)Y);
1144:   } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1145:     MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE);
1146:     MatAXPY_Basic(Y, a, X, str);
1147:     MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE);
1148:   } else {
1149:     Mat       B;
1150:     PetscInt *nnz;
1152:     MatGetRowUpperTriangular(X);
1153:     MatGetRowUpperTriangular(Y);
1154:     PetscMalloc1(Y->rmap->N, &nnz);
1155:     MatCreate(PetscObjectComm((PetscObject)Y), &B);
1156:     PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name);
1157:     MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N);
1158:     MatSetBlockSizesFromMats(B, Y, Y);
1159:     MatSetType(B, ((PetscObject)Y)->type_name);
1160:     MatAXPYGetPreallocation_SeqSBAIJ(Y, X, nnz);
1161:     MatSeqSBAIJSetPreallocation(B, bs, 0, nnz);

1163:     MatAXPY_BasicWithPreallocation(B, Y, a, X, str);

1165:     MatHeaderMerge(Y, &B);
1166:     PetscFree(nnz);
1167:     MatRestoreRowUpperTriangular(X);
1168:     MatRestoreRowUpperTriangular(Y);
1169:   }
1170:   return 0;
1171: }

1173: PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A, PetscReal tol, PetscBool *flg)
1174: {
1175:   *flg = PETSC_TRUE;
1176:   return 0;
1177: }

1179: PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A, PetscBool *flg)
1180: {
1181:   *flg = PETSC_TRUE;
1182:   return 0;
1183: }

1185: PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A, PetscReal tol, PetscBool *flg)
1186: {
1187:   *flg = PETSC_FALSE;
1188:   return 0;
1189: }

1191: PetscErrorCode MatConjugate_SeqSBAIJ(Mat A)
1192: {
1193: #if defined(PETSC_USE_COMPLEX)
1194:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1195:   PetscInt      i, nz = a->bs2 * a->i[a->mbs];
1196:   MatScalar    *aa = a->a;

1198:   for (i = 0; i < nz; i++) aa[i] = PetscConj(aa[i]);
1199: #else
1200: #endif
1201:   return 0;
1202: }

1204: PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1205: {
1206:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1207:   PetscInt      i, nz = a->bs2 * a->i[a->mbs];
1208:   MatScalar    *aa = a->a;

1210:   for (i = 0; i < nz; i++) aa[i] = PetscRealPart(aa[i]);
1211:   return 0;
1212: }

1214: PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1215: {
1216:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1217:   PetscInt      i, nz = a->bs2 * a->i[a->mbs];
1218:   MatScalar    *aa = a->a;

1220:   for (i = 0; i < nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1221:   return 0;
1222: }

1224: PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A, PetscInt is_n, const PetscInt is_idx[], PetscScalar diag, Vec x, Vec b)
1225: {
1226:   Mat_SeqSBAIJ      *baij = (Mat_SeqSBAIJ *)A->data;
1227:   PetscInt           i, j, k, count;
1228:   PetscInt           bs = A->rmap->bs, bs2 = baij->bs2, row, col;
1229:   PetscScalar        zero = 0.0;
1230:   MatScalar         *aa;
1231:   const PetscScalar *xx;
1232:   PetscScalar       *bb;
1233:   PetscBool         *zeroed, vecs = PETSC_FALSE;

1235:   /* fix right hand side if needed */
1236:   if (x && b) {
1237:     VecGetArrayRead(x, &xx);
1238:     VecGetArray(b, &bb);
1239:     vecs = PETSC_TRUE;
1240:   }

1242:   /* zero the columns */
1243:   PetscCalloc1(A->rmap->n, &zeroed);
1244:   for (i = 0; i < is_n; i++) {
1246:     zeroed[is_idx[i]] = PETSC_TRUE;
1247:   }
1248:   if (vecs) {
1249:     for (i = 0; i < A->rmap->N; i++) {
1250:       row = i / bs;
1251:       for (j = baij->i[row]; j < baij->i[row + 1]; j++) {
1252:         for (k = 0; k < bs; k++) {
1253:           col = bs * baij->j[j] + k;
1254:           if (col <= i) continue;
1255:           aa = ((MatScalar *)(baij->a)) + j * bs2 + (i % bs) + bs * k;
1256:           if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0] * xx[col];
1257:           if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0] * xx[i];
1258:         }
1259:       }
1260:     }
1261:     for (i = 0; i < is_n; i++) bb[is_idx[i]] = diag * xx[is_idx[i]];
1262:   }

1264:   for (i = 0; i < A->rmap->N; i++) {
1265:     if (!zeroed[i]) {
1266:       row = i / bs;
1267:       for (j = baij->i[row]; j < baij->i[row + 1]; j++) {
1268:         for (k = 0; k < bs; k++) {
1269:           col = bs * baij->j[j] + k;
1270:           if (zeroed[col]) {
1271:             aa    = ((MatScalar *)(baij->a)) + j * bs2 + (i % bs) + bs * k;
1272:             aa[0] = 0.0;
1273:           }
1274:         }
1275:       }
1276:     }
1277:   }
1278:   PetscFree(zeroed);
1279:   if (vecs) {
1280:     VecRestoreArrayRead(x, &xx);
1281:     VecRestoreArray(b, &bb);
1282:   }

1284:   /* zero the rows */
1285:   for (i = 0; i < is_n; i++) {
1286:     row   = is_idx[i];
1287:     count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
1288:     aa    = ((MatScalar *)(baij->a)) + baij->i[row / bs] * bs2 + (row % bs);
1289:     for (k = 0; k < count; k++) {
1290:       aa[0] = zero;
1291:       aa += bs;
1292:     }
1293:     if (diag != 0.0) PetscUseTypeMethod(A, setvalues, 1, &row, 1, &row, &diag, INSERT_VALUES);
1294:   }
1295:   MatAssemblyEnd_SeqSBAIJ(A, MAT_FINAL_ASSEMBLY);
1296:   return 0;
1297: }

1299: PetscErrorCode MatShift_SeqSBAIJ(Mat Y, PetscScalar a)
1300: {
1301:   Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)Y->data;

1303:   if (!Y->preallocated || !aij->nz) MatSeqSBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL);
1304:   MatShift_Basic(Y, a);
1305:   return 0;
1306: }

1308: /* -------------------------------------------------------------------*/
1309: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1310:                                        MatGetRow_SeqSBAIJ,
1311:                                        MatRestoreRow_SeqSBAIJ,
1312:                                        MatMult_SeqSBAIJ_N,
1313:                                        /*  4*/ MatMultAdd_SeqSBAIJ_N,
1314:                                        MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1315:                                        MatMultAdd_SeqSBAIJ_N,
1316:                                        NULL,
1317:                                        NULL,
1318:                                        NULL,
1319:                                        /* 10*/ NULL,
1320:                                        NULL,
1321:                                        MatCholeskyFactor_SeqSBAIJ,
1322:                                        MatSOR_SeqSBAIJ,
1323:                                        MatTranspose_SeqSBAIJ,
1324:                                        /* 15*/ MatGetInfo_SeqSBAIJ,
1325:                                        MatEqual_SeqSBAIJ,
1326:                                        MatGetDiagonal_SeqSBAIJ,
1327:                                        MatDiagonalScale_SeqSBAIJ,
1328:                                        MatNorm_SeqSBAIJ,
1329:                                        /* 20*/ NULL,
1330:                                        MatAssemblyEnd_SeqSBAIJ,
1331:                                        MatSetOption_SeqSBAIJ,
1332:                                        MatZeroEntries_SeqSBAIJ,
1333:                                        /* 24*/ NULL,
1334:                                        NULL,
1335:                                        NULL,
1336:                                        NULL,
1337:                                        NULL,
1338:                                        /* 29*/ MatSetUp_SeqSBAIJ,
1339:                                        NULL,
1340:                                        NULL,
1341:                                        NULL,
1342:                                        NULL,
1343:                                        /* 34*/ MatDuplicate_SeqSBAIJ,
1344:                                        NULL,
1345:                                        NULL,
1346:                                        NULL,
1347:                                        MatICCFactor_SeqSBAIJ,
1348:                                        /* 39*/ MatAXPY_SeqSBAIJ,
1349:                                        MatCreateSubMatrices_SeqSBAIJ,
1350:                                        MatIncreaseOverlap_SeqSBAIJ,
1351:                                        MatGetValues_SeqSBAIJ,
1352:                                        MatCopy_SeqSBAIJ,
1353:                                        /* 44*/ NULL,
1354:                                        MatScale_SeqSBAIJ,
1355:                                        MatShift_SeqSBAIJ,
1356:                                        NULL,
1357:                                        MatZeroRowsColumns_SeqSBAIJ,
1358:                                        /* 49*/ NULL,
1359:                                        MatGetRowIJ_SeqSBAIJ,
1360:                                        MatRestoreRowIJ_SeqSBAIJ,
1361:                                        NULL,
1362:                                        NULL,
1363:                                        /* 54*/ NULL,
1364:                                        NULL,
1365:                                        NULL,
1366:                                        MatPermute_SeqSBAIJ,
1367:                                        MatSetValuesBlocked_SeqSBAIJ,
1368:                                        /* 59*/ MatCreateSubMatrix_SeqSBAIJ,
1369:                                        NULL,
1370:                                        NULL,
1371:                                        NULL,
1372:                                        NULL,
1373:                                        /* 64*/ NULL,
1374:                                        NULL,
1375:                                        NULL,
1376:                                        NULL,
1377:                                        NULL,
1378:                                        /* 69*/ MatGetRowMaxAbs_SeqSBAIJ,
1379:                                        NULL,
1380:                                        MatConvert_MPISBAIJ_Basic,
1381:                                        NULL,
1382:                                        NULL,
1383:                                        /* 74*/ NULL,
1384:                                        NULL,
1385:                                        NULL,
1386:                                        NULL,
1387:                                        NULL,
1388:                                        /* 79*/ NULL,
1389:                                        NULL,
1390:                                        NULL,
1391:                                        MatGetInertia_SeqSBAIJ,
1392:                                        MatLoad_SeqSBAIJ,
1393:                                        /* 84*/ MatIsSymmetric_SeqSBAIJ,
1394:                                        MatIsHermitian_SeqSBAIJ,
1395:                                        MatIsStructurallySymmetric_SeqSBAIJ,
1396:                                        NULL,
1397:                                        NULL,
1398:                                        /* 89*/ NULL,
1399:                                        NULL,
1400:                                        NULL,
1401:                                        NULL,
1402:                                        NULL,
1403:                                        /* 94*/ NULL,
1404:                                        NULL,
1405:                                        NULL,
1406:                                        NULL,
1407:                                        NULL,
1408:                                        /* 99*/ NULL,
1409:                                        NULL,
1410:                                        NULL,
1411:                                        MatConjugate_SeqSBAIJ,
1412:                                        NULL,
1413:                                        /*104*/ NULL,
1414:                                        MatRealPart_SeqSBAIJ,
1415:                                        MatImaginaryPart_SeqSBAIJ,
1416:                                        MatGetRowUpperTriangular_SeqSBAIJ,
1417:                                        MatRestoreRowUpperTriangular_SeqSBAIJ,
1418:                                        /*109*/ NULL,
1419:                                        NULL,
1420:                                        NULL,
1421:                                        NULL,
1422:                                        MatMissingDiagonal_SeqSBAIJ,
1423:                                        /*114*/ NULL,
1424:                                        NULL,
1425:                                        NULL,
1426:                                        NULL,
1427:                                        NULL,
1428:                                        /*119*/ NULL,
1429:                                        NULL,
1430:                                        NULL,
1431:                                        NULL,
1432:                                        NULL,
1433:                                        /*124*/ NULL,
1434:                                        NULL,
1435:                                        NULL,
1436:                                        NULL,
1437:                                        NULL,
1438:                                        /*129*/ NULL,
1439:                                        NULL,
1440:                                        NULL,
1441:                                        NULL,
1442:                                        NULL,
1443:                                        /*134*/ NULL,
1444:                                        NULL,
1445:                                        NULL,
1446:                                        NULL,
1447:                                        NULL,
1448:                                        /*139*/ MatSetBlockSizes_Default,
1449:                                        NULL,
1450:                                        NULL,
1451:                                        NULL,
1452:                                        NULL,
1453:                                        /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ,
1454:                                        NULL,
1455:                                        NULL,
1456:                                        NULL,
1457:                                        NULL,
1458:                                        NULL,
1459:                                        /*150*/ NULL};

1461: PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1462: {
1463:   Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data;
1464:   PetscInt      nz  = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2;


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

1471:   /* copy values over */
1472:   PetscArraycpy(aij->saved_values, aij->a, nz);
1473:   return 0;
1474: }

1476: PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1477: {
1478:   Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data;
1479:   PetscInt      nz  = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2;


1484:   /* copy values over */
1485:   PetscArraycpy(aij->a, aij->saved_values, nz);
1486:   return 0;
1487: }

1489: static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B, PetscInt bs, PetscInt nz, PetscInt *nnz)
1490: {
1491:   Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data;
1492:   PetscInt      i, mbs, nbs, bs2;
1493:   PetscBool     skipallocation = PETSC_FALSE, flg = PETSC_FALSE, realalloc = PETSC_FALSE;

1495:   if (nz >= 0 || nnz) realalloc = PETSC_TRUE;

1497:   MatSetBlockSize(B, PetscAbs(bs));
1498:   PetscLayoutSetUp(B->rmap);
1499:   PetscLayoutSetUp(B->cmap);
1501:   PetscLayoutGetBlockSize(B->rmap, &bs);

1503:   B->preallocated = PETSC_TRUE;

1505:   mbs = B->rmap->N / bs;
1506:   nbs = B->cmap->n / bs;
1507:   bs2 = bs * bs;


1511:   if (nz == MAT_SKIP_ALLOCATION) {
1512:     skipallocation = PETSC_TRUE;
1513:     nz             = 0;
1514:   }

1516:   if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1518:   if (nnz) {
1519:     for (i = 0; i < mbs; i++) {
1522:     }
1523:   }

1525:   B->ops->mult             = MatMult_SeqSBAIJ_N;
1526:   B->ops->multadd          = MatMultAdd_SeqSBAIJ_N;
1527:   B->ops->multtranspose    = MatMult_SeqSBAIJ_N;
1528:   B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N;

1530:   PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL);
1531:   if (!flg) {
1532:     switch (bs) {
1533:     case 1:
1534:       B->ops->mult             = MatMult_SeqSBAIJ_1;
1535:       B->ops->multadd          = MatMultAdd_SeqSBAIJ_1;
1536:       B->ops->multtranspose    = MatMult_SeqSBAIJ_1;
1537:       B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1538:       break;
1539:     case 2:
1540:       B->ops->mult             = MatMult_SeqSBAIJ_2;
1541:       B->ops->multadd          = MatMultAdd_SeqSBAIJ_2;
1542:       B->ops->multtranspose    = MatMult_SeqSBAIJ_2;
1543:       B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1544:       break;
1545:     case 3:
1546:       B->ops->mult             = MatMult_SeqSBAIJ_3;
1547:       B->ops->multadd          = MatMultAdd_SeqSBAIJ_3;
1548:       B->ops->multtranspose    = MatMult_SeqSBAIJ_3;
1549:       B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1550:       break;
1551:     case 4:
1552:       B->ops->mult             = MatMult_SeqSBAIJ_4;
1553:       B->ops->multadd          = MatMultAdd_SeqSBAIJ_4;
1554:       B->ops->multtranspose    = MatMult_SeqSBAIJ_4;
1555:       B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1556:       break;
1557:     case 5:
1558:       B->ops->mult             = MatMult_SeqSBAIJ_5;
1559:       B->ops->multadd          = MatMultAdd_SeqSBAIJ_5;
1560:       B->ops->multtranspose    = MatMult_SeqSBAIJ_5;
1561:       B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1562:       break;
1563:     case 6:
1564:       B->ops->mult             = MatMult_SeqSBAIJ_6;
1565:       B->ops->multadd          = MatMultAdd_SeqSBAIJ_6;
1566:       B->ops->multtranspose    = MatMult_SeqSBAIJ_6;
1567:       B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1568:       break;
1569:     case 7:
1570:       B->ops->mult             = MatMult_SeqSBAIJ_7;
1571:       B->ops->multadd          = MatMultAdd_SeqSBAIJ_7;
1572:       B->ops->multtranspose    = MatMult_SeqSBAIJ_7;
1573:       B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1574:       break;
1575:     }
1576:   }

1578:   b->mbs = mbs;
1579:   b->nbs = nbs;
1580:   if (!skipallocation) {
1581:     if (!b->imax) {
1582:       PetscMalloc2(mbs, &b->imax, mbs, &b->ilen);

1584:       b->free_imax_ilen = PETSC_TRUE;
1585:     }
1586:     if (!nnz) {
1587:       if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1588:       else if (nz <= 0) nz = 1;
1589:       nz = PetscMin(nbs, nz);
1590:       for (i = 0; i < mbs; i++) b->imax[i] = nz;
1591:       PetscIntMultError(nz, mbs, &nz);
1592:     } else {
1593:       PetscInt64 nz64 = 0;
1594:       for (i = 0; i < mbs; i++) {
1595:         b->imax[i] = nnz[i];
1596:         nz64 += nnz[i];
1597:       }
1598:       PetscIntCast(nz64, &nz);
1599:     }
1600:     /* b->ilen will count nonzeros in each block row so far. */
1601:     for (i = 0; i < mbs; i++) b->ilen[i] = 0;
1602:     /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */

1604:     /* allocate the matrix space */
1605:     MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i);
1606:     PetscMalloc3(bs2 * nz, &b->a, nz, &b->j, B->rmap->N + 1, &b->i);
1607:     PetscArrayzero(b->a, nz * bs2);
1608:     PetscArrayzero(b->j, nz);

1610:     b->singlemalloc = PETSC_TRUE;

1612:     /* pointer to beginning of each row */
1613:     b->i[0] = 0;
1614:     for (i = 1; i < mbs + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1];

1616:     b->free_a  = PETSC_TRUE;
1617:     b->free_ij = PETSC_TRUE;
1618:   } else {
1619:     b->free_a  = PETSC_FALSE;
1620:     b->free_ij = PETSC_FALSE;
1621:   }

1623:   b->bs2     = bs2;
1624:   b->nz      = 0;
1625:   b->maxnz   = nz;
1626:   b->inew    = NULL;
1627:   b->jnew    = NULL;
1628:   b->anew    = NULL;
1629:   b->a2anew  = NULL;
1630:   b->permute = PETSC_FALSE;

1632:   B->was_assembled = PETSC_FALSE;
1633:   B->assembled     = PETSC_FALSE;
1634:   if (realalloc) MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
1635:   return 0;
1636: }

1638: PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[])
1639: {
1640:   PetscInt     i, j, m, nz, anz, nz_max = 0, *nnz;
1641:   PetscScalar *values      = NULL;
1642:   PetscBool    roworiented = ((Mat_SeqSBAIJ *)B->data)->roworiented;

1645:   PetscLayoutSetBlockSize(B->rmap, bs);
1646:   PetscLayoutSetBlockSize(B->cmap, bs);
1647:   PetscLayoutSetUp(B->rmap);
1648:   PetscLayoutSetUp(B->cmap);
1649:   PetscLayoutGetBlockSize(B->rmap, &bs);
1650:   m = B->rmap->n / bs;

1653:   PetscMalloc1(m + 1, &nnz);
1654:   for (i = 0; i < m; i++) {
1655:     nz = ii[i + 1] - ii[i];
1657:     anz = 0;
1658:     for (j = 0; j < nz; j++) {
1659:       /* count only values on the diagonal or above */
1660:       if (jj[ii[i] + j] >= i) {
1661:         anz = nz - j;
1662:         break;
1663:       }
1664:     }
1665:     nz_max = PetscMax(nz_max, anz);
1666:     nnz[i] = anz;
1667:   }
1668:   MatSeqSBAIJSetPreallocation(B, bs, 0, nnz);
1669:   PetscFree(nnz);

1671:   values = (PetscScalar *)V;
1672:   if (!values) PetscCalloc1(bs * bs * nz_max, &values);
1673:   for (i = 0; i < m; i++) {
1674:     PetscInt        ncols = ii[i + 1] - ii[i];
1675:     const PetscInt *icols = jj + ii[i];
1676:     if (!roworiented || bs == 1) {
1677:       const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0);
1678:       MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, ncols, icols, svals, INSERT_VALUES);
1679:     } else {
1680:       for (j = 0; j < ncols; j++) {
1681:         const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0);
1682:         MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, 1, &icols[j], svals, INSERT_VALUES);
1683:       }
1684:     }
1685:   }
1686:   if (!V) PetscFree(values);
1687:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
1688:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
1689:   MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
1690:   return 0;
1691: }

1693: /*
1694:    This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization
1695: */
1696: PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B, PetscBool natural)
1697: {
1698:   PetscBool flg = PETSC_FALSE;
1699:   PetscInt  bs  = B->rmap->bs;

1701:   PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL);
1702:   if (flg) bs = 8;

1704:   if (!natural) {
1705:     switch (bs) {
1706:     case 1:
1707:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace;
1708:       break;
1709:     case 2:
1710:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1711:       break;
1712:     case 3:
1713:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1714:       break;
1715:     case 4:
1716:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1717:       break;
1718:     case 5:
1719:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1720:       break;
1721:     case 6:
1722:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1723:       break;
1724:     case 7:
1725:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1726:       break;
1727:     default:
1728:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N;
1729:       break;
1730:     }
1731:   } else {
1732:     switch (bs) {
1733:     case 1:
1734:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace;
1735:       break;
1736:     case 2:
1737:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1738:       break;
1739:     case 3:
1740:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1741:       break;
1742:     case 4:
1743:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1744:       break;
1745:     case 5:
1746:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1747:       break;
1748:     case 6:
1749:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1750:       break;
1751:     case 7:
1752:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1753:       break;
1754:     default:
1755:       B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
1756:       break;
1757:     }
1758:   }
1759:   return 0;
1760: }

1762: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *);
1763: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType, MatReuse, Mat *);
1764: static PetscErrorCode       MatFactorGetSolverType_petsc(Mat A, MatSolverType *type)
1765: {
1766:   *type = MATSOLVERPETSC;
1767:   return 0;
1768: }

1770: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A, MatFactorType ftype, Mat *B)
1771: {
1772:   PetscInt n = A->rmap->n;

1774: #if defined(PETSC_USE_COMPLEX)
1776: #endif

1778:   MatCreate(PetscObjectComm((PetscObject)A), B);
1779:   MatSetSizes(*B, n, n, n, n);
1780:   if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1781:     MatSetType(*B, MATSEQSBAIJ);
1782:     MatSeqSBAIJSetPreallocation(*B, A->rmap->bs, MAT_SKIP_ALLOCATION, NULL);

1784:     (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ;
1785:     (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqSBAIJ;
1786:     PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);
1787:     PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC]);
1788:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");

1790:   (*B)->factortype     = ftype;
1791:   (*B)->canuseordering = PETSC_TRUE;
1792:   PetscFree((*B)->solvertype);
1793:   PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype);
1794:   PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_petsc);
1795:   return 0;
1796: }

1798: /*@C
1799:    MatSeqSBAIJGetArray - gives access to the array where the data for a `MATSEQSBAIJ` matrix is stored

1801:    Not Collective

1803:    Input Parameter:
1804: .  mat - a `MATSEQSBAIJ` matrix

1806:    Output Parameter:
1807: .   array - pointer to the data

1809:    Level: intermediate

1811: .seealso: `MATSEQSBAIJ`, `MatSeqSBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()`
1812: @*/
1813: PetscErrorCode MatSeqSBAIJGetArray(Mat A, PetscScalar **array)
1814: {
1815:   PetscUseMethod(A, "MatSeqSBAIJGetArray_C", (Mat, PetscScalar **), (A, array));
1816:   return 0;
1817: }

1819: /*@C
1820:    MatSeqSBAIJRestoreArray - returns access to the array where the data for a `MATSEQSBAIJ` matrix is stored obtained by `MatSeqSBAIJGetArray()`

1822:    Not Collective

1824:    Input Parameters:
1825: +  mat - a MATSEQSBAIJ matrix
1826: -  array - pointer to the data

1828:    Level: intermediate

1830: .seealso: `MATSEQSBAIJ`, `MatSeqSBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()`
1831: @*/
1832: PetscErrorCode MatSeqSBAIJRestoreArray(Mat A, PetscScalar **array)
1833: {
1834:   PetscUseMethod(A, "MatSeqSBAIJRestoreArray_C", (Mat, PetscScalar **), (A, array));
1835:   return 0;
1836: }

1838: /*MC
1839:   MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1840:   based on block compressed sparse row format.  Only the upper triangular portion of the matrix is stored.

1842:   For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
1843:   can call `MatSetOption`(`Mat`, `MAT_HERMITIAN`).

1845:   Options Database Keys:
1846:   . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to `MatSetFromOptions()`

1848:   Notes:
1849:     By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not
1850:      stored and it is assumed they symmetric to the upper triangular). If you call `MatSetOption`(`Mat`,`MAT_IGNORE_LOWER_TRIANGULAR`,`PETSC_FALSE`) or use
1851:      the options database -mat_ignore_lower_triangular false it will generate an error if you try to set a value in the lower triangular portion.

1853:     The number of rows in the matrix must be less than or equal to the number of columns

1855:   Level: beginner

1857:   .seealso: `MATSEQSBAIJ`, `MatCreateSeqSBAIJ()`, `MatType`, `MATMPISBAIJ`
1858: M*/
1859: PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1860: {
1861:   Mat_SeqSBAIJ *b;
1862:   PetscMPIInt   size;
1863:   PetscBool     no_unroll = PETSC_FALSE, no_inode = PETSC_FALSE;

1865:   MPI_Comm_size(PetscObjectComm((PetscObject)B), &size);

1868:   PetscNew(&b);
1869:   B->data = (void *)b;
1870:   PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps));

1872:   B->ops->destroy    = MatDestroy_SeqSBAIJ;
1873:   B->ops->view       = MatView_SeqSBAIJ;
1874:   b->row             = NULL;
1875:   b->icol            = NULL;
1876:   b->reallocs        = 0;
1877:   b->saved_values    = NULL;
1878:   b->inode.limit     = 5;
1879:   b->inode.max_limit = 5;

1881:   b->roworiented        = PETSC_TRUE;
1882:   b->nonew              = 0;
1883:   b->diag               = NULL;
1884:   b->solve_work         = NULL;
1885:   b->mult_work          = NULL;
1886:   B->spptr              = NULL;
1887:   B->info.nz_unneeded   = (PetscReal)b->maxnz * b->bs2;
1888:   b->keepnonzeropattern = PETSC_FALSE;

1890:   b->inew    = NULL;
1891:   b->jnew    = NULL;
1892:   b->anew    = NULL;
1893:   b->a2anew  = NULL;
1894:   b->permute = PETSC_FALSE;

1896:   b->ignore_ltriangular = PETSC_TRUE;

1898:   PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_ignore_lower_triangular", &b->ignore_ltriangular, NULL);

1900:   b->getrow_utriangular = PETSC_FALSE;

1902:   PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_getrow_uppertriangular", &b->getrow_utriangular, NULL);

1904:   PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJGetArray_C", MatSeqSBAIJGetArray_SeqSBAIJ);
1905:   PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJRestoreArray_C", MatSeqSBAIJRestoreArray_SeqSBAIJ);
1906:   PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqSBAIJ);
1907:   PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqSBAIJ);
1908:   PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetColumnIndices_C", MatSeqSBAIJSetColumnIndices_SeqSBAIJ);
1909:   PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqaij_C", MatConvert_SeqSBAIJ_SeqAIJ);
1910:   PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqbaij_C", MatConvert_SeqSBAIJ_SeqBAIJ);
1911:   PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocation_C", MatSeqSBAIJSetPreallocation_SeqSBAIJ);
1912:   PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocationCSR_C", MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);
1913: #if defined(PETSC_HAVE_ELEMENTAL)
1914:   PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_elemental_C", MatConvert_SeqSBAIJ_Elemental);
1915: #endif
1916: #if defined(PETSC_HAVE_SCALAPACK)
1917:   PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_scalapack_C", MatConvert_SBAIJ_ScaLAPACK);
1918: #endif

1920:   B->symmetry_eternal            = PETSC_TRUE;
1921:   B->structural_symmetry_eternal = PETSC_TRUE;
1922:   B->symmetric                   = PETSC_BOOL3_TRUE;
1923:   B->structurally_symmetric      = PETSC_BOOL3_TRUE;
1924: #if defined(PETSC_USE_COMPLEX)
1925:   B->hermitian = PETSC_BOOL3_FALSE;
1926: #else
1927:   B->hermitian = PETSC_BOOL3_TRUE;
1928: #endif

1930:   PetscObjectChangeTypeName((PetscObject)B, MATSEQSBAIJ);

1932:   PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "Options for SEQSBAIJ matrix", "Mat");
1933:   PetscOptionsBool("-mat_no_unroll", "Do not optimize for inodes (slower)", NULL, no_unroll, &no_unroll, NULL);
1934:   if (no_unroll) PetscInfo(B, "Not using Inode routines due to -mat_no_unroll\n");
1935:   PetscOptionsBool("-mat_no_inode", "Do not optimize for inodes (slower)", NULL, no_inode, &no_inode, NULL);
1936:   if (no_inode) PetscInfo(B, "Not using Inode routines due to -mat_no_inode\n");
1937:   PetscOptionsInt("-mat_inode_limit", "Do not use inodes larger then this value", NULL, b->inode.limit, &b->inode.limit, NULL);
1938:   PetscOptionsEnd();
1939:   b->inode.use = (PetscBool)(!(no_unroll || no_inode));
1940:   if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit;
1941:   return 0;
1942: }

1944: /*@C
1945:    MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
1946:    compressed row) `MATSEQSBAIJ` format.  For good matrix assembly performance the
1947:    user should preallocate the matrix storage by setting the parameter nz
1948:    (or the array nnz).  By setting these parameters accurately, performance
1949:    during matrix assembly can be increased by more than a factor of 50.

1951:    Collective on Mat

1953:    Input Parameters:
1954: +  B - the symmetric matrix
1955: .  bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
1956:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
1957: .  nz - number of block nonzeros per block row (same for all rows)
1958: -  nnz - array containing the number of block nonzeros in the upper triangular plus
1959:          diagonal portion of each block (possibly different for each block row) or NULL

1961:    Options Database Keys:
1962: +   -mat_no_unroll - uses code that does not unroll the loops in the
1963:                      block calculations (much slower)
1964: -   -mat_block_size - size of the blocks to use (only works if a negative bs is passed in

1966:    Level: intermediate

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

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

1978:    If the nnz parameter is given then the nz parameter is ignored

1980: .seealso: [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()`
1981: @*/
1982: PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[])
1983: {
1987:   PetscTryMethod(B, "MatSeqSBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[]), (B, bs, nz, nnz));
1988:   return 0;
1989: }

1991: /*@C
1992:    MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATSEQSBAIJ` format using the given nonzero structure and (optional) numerical values

1994:    Input Parameters:
1995: +  B - the matrix
1996: .  bs - size of block, the blocks are ALWAYS square.
1997: .  i - the indices into j for the start of each local row (starts with zero)
1998: .  j - the column indices for each local row (starts with zero) these must be sorted for each row
1999: -  v - optional values in the matrix

2001:    Level: advanced

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

2010:    Any entries below the diagonal are ignored

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

2015: .seealso: `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValuesBlocked()`, `MatSeqSBAIJSetPreallocation()`, `MATSEQSBAIJ`
2016: @*/
2017: PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[])
2018: {
2022:   PetscTryMethod(B, "MatSeqSBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v));
2023:   return 0;
2024: }

2026: /*@C
2027:    MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
2028:    compressed row) `MATSEQSBAIJ` format.  For good matrix assembly performance the
2029:    user should preallocate the matrix storage by setting the parameter nz
2030:    (or the array nnz).  By setting these parameters accurately, performance
2031:    during matrix assembly can be increased by more than a factor of 50.

2033:    Collective

2035:    Input Parameters:
2036: +  comm - MPI communicator, set to `PETSC_COMM_SELF`
2037: .  bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
2038:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2039: .  m - number of rows, or number of columns
2040: .  nz - number of block nonzeros per block row (same for all rows)
2041: -  nnz - array containing the number of block nonzeros in the upper triangular plus
2042:          diagonal portion of each block (possibly different for each block row) or NULL

2044:    Output Parameter:
2045: .  A - the symmetric matrix

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

2052:    Level: intermediate

2054:    It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
2055:    MatXXXXSetPreallocation() paradigm instead of this routine directly.
2056:    [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`]

2058:    Notes:
2059:    The number of rows and columns must be divisible by blocksize.
2060:    This matrix type does not support complex Hermitian operation.

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

2066:    If the nnz parameter is given then the nz parameter is ignored

2068: .seealso: [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()`
2069: @*/
2070: PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A)
2071: {
2072:   MatCreate(comm, A);
2073:   MatSetSizes(*A, m, n, m, n);
2074:   MatSetType(*A, MATSEQSBAIJ);
2075:   MatSeqSBAIJSetPreallocation(*A, bs, nz, (PetscInt *)nnz);
2076:   return 0;
2077: }

2079: PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B)
2080: {
2081:   Mat           C;
2082:   Mat_SeqSBAIJ *c, *a  = (Mat_SeqSBAIJ *)A->data;
2083:   PetscInt      i, mbs = a->mbs, nz = a->nz, bs2 = a->bs2;


2087:   *B = NULL;
2088:   MatCreate(PetscObjectComm((PetscObject)A), &C);
2089:   MatSetSizes(C, A->rmap->N, A->cmap->n, A->rmap->N, A->cmap->n);
2090:   MatSetBlockSizesFromMats(C, A, A);
2091:   MatSetType(C, MATSEQSBAIJ);
2092:   c = (Mat_SeqSBAIJ *)C->data;

2094:   C->preallocated       = PETSC_TRUE;
2095:   C->factortype         = A->factortype;
2096:   c->row                = NULL;
2097:   c->icol               = NULL;
2098:   c->saved_values       = NULL;
2099:   c->keepnonzeropattern = a->keepnonzeropattern;
2100:   C->assembled          = PETSC_TRUE;

2102:   PetscLayoutReference(A->rmap, &C->rmap);
2103:   PetscLayoutReference(A->cmap, &C->cmap);
2104:   c->bs2 = a->bs2;
2105:   c->mbs = a->mbs;
2106:   c->nbs = a->nbs;

2108:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2109:     c->imax           = a->imax;
2110:     c->ilen           = a->ilen;
2111:     c->free_imax_ilen = PETSC_FALSE;
2112:   } else {
2113:     PetscMalloc2((mbs + 1), &c->imax, (mbs + 1), &c->ilen);
2114:     for (i = 0; i < mbs; i++) {
2115:       c->imax[i] = a->imax[i];
2116:       c->ilen[i] = a->ilen[i];
2117:     }
2118:     c->free_imax_ilen = PETSC_TRUE;
2119:   }

2121:   /* allocate the matrix space */
2122:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2123:     PetscMalloc1(bs2 * nz, &c->a);
2124:     c->i            = a->i;
2125:     c->j            = a->j;
2126:     c->singlemalloc = PETSC_FALSE;
2127:     c->free_a       = PETSC_TRUE;
2128:     c->free_ij      = PETSC_FALSE;
2129:     c->parent       = A;
2130:     PetscObjectReference((PetscObject)A);
2131:     MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
2132:     MatSetOption(C, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
2133:   } else {
2134:     PetscMalloc3(bs2 * nz, &c->a, nz, &c->j, mbs + 1, &c->i);
2135:     PetscArraycpy(c->i, a->i, mbs + 1);
2136:     c->singlemalloc = PETSC_TRUE;
2137:     c->free_a       = PETSC_TRUE;
2138:     c->free_ij      = PETSC_TRUE;
2139:   }
2140:   if (mbs > 0) {
2141:     if (cpvalues != MAT_SHARE_NONZERO_PATTERN) PetscArraycpy(c->j, a->j, nz);
2142:     if (cpvalues == MAT_COPY_VALUES) {
2143:       PetscArraycpy(c->a, a->a, bs2 * nz);
2144:     } else {
2145:       PetscArrayzero(c->a, bs2 * nz);
2146:     }
2147:     if (a->jshort) {
2148:       /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */
2149:       /* if the parent matrix is reassembled, this child matrix will never notice */
2150:       PetscMalloc1(nz, &c->jshort);
2151:       PetscArraycpy(c->jshort, a->jshort, nz);

2153:       c->free_jshort = PETSC_TRUE;
2154:     }
2155:   }

2157:   c->roworiented = a->roworiented;
2158:   c->nonew       = a->nonew;

2160:   if (a->diag) {
2161:     if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2162:       c->diag      = a->diag;
2163:       c->free_diag = PETSC_FALSE;
2164:     } else {
2165:       PetscMalloc1(mbs, &c->diag);
2166:       for (i = 0; i < mbs; i++) c->diag[i] = a->diag[i];
2167:       c->free_diag = PETSC_TRUE;
2168:     }
2169:   }
2170:   c->nz         = a->nz;
2171:   c->maxnz      = a->nz; /* Since we allocate exactly the right amount */
2172:   c->solve_work = NULL;
2173:   c->mult_work  = NULL;

2175:   *B = C;
2176:   PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist);
2177:   return 0;
2178: }

2180: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
2181: #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary

2183: PetscErrorCode MatLoad_SeqSBAIJ(Mat mat, PetscViewer viewer)
2184: {
2185:   PetscBool isbinary;

2187:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
2189:   MatLoad_SeqSBAIJ_Binary(mat, viewer);
2190:   return 0;
2191: }

2193: /*@
2194:      MatCreateSeqSBAIJWithArrays - Creates an sequential `MATSEQSBAIJ` matrix using matrix elements
2195:               (upper triangular entries in CSR format) provided by the user.

2197:      Collective

2199:    Input Parameters:
2200: +  comm - must be an MPI communicator of size 1
2201: .  bs - size of block
2202: .  m - number of rows
2203: .  n - number of columns
2204: .  i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that row block row of the matrix
2205: .  j - column indices
2206: -  a - matrix values

2208:    Output Parameter:
2209: .  mat - the matrix

2211:    Level: advanced

2213:    Notes:
2214:        The i, j, and a arrays are not copied by this routine, the user must free these arrays
2215:     once the matrix is destroyed

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

2219:        The i and j indices are 0 based

2221:        When block size is greater than 1 the matrix values must be stored using the SBAIJ storage format (see the SBAIJ source code to determine this). For block size of 1
2222:        it is the regular CSR format excluding the lower triangular elements.

2224: .seealso: `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSBAIJ()`, `MatCreateSeqSBAIJ()`
2225: @*/
2226: PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat)
2227: {
2228:   PetscInt      ii;
2229:   Mat_SeqSBAIJ *sbaij;


2234:   MatCreate(comm, mat);
2235:   MatSetSizes(*mat, m, n, m, n);
2236:   MatSetType(*mat, MATSEQSBAIJ);
2237:   MatSeqSBAIJSetPreallocation(*mat, bs, MAT_SKIP_ALLOCATION, NULL);
2238:   sbaij = (Mat_SeqSBAIJ *)(*mat)->data;
2239:   PetscMalloc2(m, &sbaij->imax, m, &sbaij->ilen);

2241:   sbaij->i = i;
2242:   sbaij->j = j;
2243:   sbaij->a = a;

2245:   sbaij->singlemalloc   = PETSC_FALSE;
2246:   sbaij->nonew          = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2247:   sbaij->free_a         = PETSC_FALSE;
2248:   sbaij->free_ij        = PETSC_FALSE;
2249:   sbaij->free_imax_ilen = PETSC_TRUE;

2251:   for (ii = 0; ii < m; ii++) {
2252:     sbaij->ilen[ii] = sbaij->imax[ii] = i[ii + 1] - i[ii];
2254:   }
2255:   if (PetscDefined(USE_DEBUG)) {
2256:     for (ii = 0; ii < sbaij->i[m]; ii++) {
2259:     }
2260:   }

2262:   MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
2263:   MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
2264:   return 0;
2265: }

2267: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
2268: {
2269:   MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm, inmat, n, scall, outmat);
2270:   return 0;
2271: }