1: static char help[] = "Regression test for the MatSOR_SeqAIJ_Inode() block-diagonal cache.\n\n\
  2:    This test catches the typo in MatInvertDiagonalForSOR_SeqAIJ_Inode() introduced\n\
  3:    in MR !8797 where the cache validity check at inode.c:2432 reads a->idiagState\n\
  4:    (the point-SOR cache token) instead of a->inode.ibdiagState (the inode block-SOR\n\
  5:    cache token).  With the bug, the inverted block diagonal is recomputed on every\n\
  6:    MatSOR() call, producing a runtime regression but not a correctness failure --\n\
  7:    so we observe the cache invariant directly by poisoning the cached buffer and\n\
  8:    checking whether it survives a second MatSOR() call on an unchanged matrix.\n\n";

 10: #include <petscmat.h>

 12: /* The bug lives in private inode-cache fields.  Including the private header is
 13:    the cleanest way to test the cache invariant; this is the established practice
 14:    for tests that target a specific implementation detail rather than public API. */
 15: #include <../src/mat/impls/aij/seq/aij.h>

 17: int main(int argc, char **args)
 18: {
 19:   Mat         A;
 20:   Vec         b, x_pristine, x_poisoned;
 21:   Mat_SeqAIJ *a;
 22:   PetscInt    nblock = 20, blocksize = 3, n;
 23:   PetscBool   same;

 25:   PetscFunctionBeginUser;
 26:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 27:   n = nblock * blocksize;

 29:   PetscCall(MatCreate(PETSC_COMM_SELF, &A));
 30:   PetscCall(MatSetSizes(A, n, n, n, n));
 31:   PetscCall(MatSetType(A, MATSEQAIJ));
 32:   PetscCall(MatSeqAIJSetPreallocation(A, 3 * blocksize, NULL));

 34:   /* Block-tridiagonal matrix: each row has one diagonal entry of value 4.0
 35:      and -1.0 entries at every other column in its own block and the two
 36:      neighbouring blocks.  All `blocksize` rows of a block share the same column
 37:      pattern, so MatSeqAIJCheckInode() forms `nblock` inodes of size `blocksize`. */
 38:   for (PetscInt Ib = 0; Ib < nblock; Ib++) {
 39:     PetscInt    Jlo = PetscMax(Ib - 1, 0), Jhi = PetscMin(Ib + 1, nblock - 1);
 40:     PetscInt    ncols = (Jhi - Jlo + 1) * blocksize, cols[9];
 41:     PetscScalar vals[9];

 43:     for (PetscInt Jb = Jlo, k = 0; Jb <= Jhi; Jb++)
 44:       for (PetscInt jj = 0; jj < blocksize; jj++, k++) cols[k] = Jb * blocksize + jj;

 46:     for (PetscInt ii = 0; ii < blocksize; ii++) {
 47:       PetscInt row = Ib * blocksize + ii;

 49:       for (PetscInt k = 0; k < ncols; k++) vals[k] = (cols[k] == row) ? 4.0 : -1.0;
 50:       PetscCall(MatSetValues(A, 1, &row, ncols, cols, vals, INSERT_VALUES));
 51:     }
 52:   }
 53:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 54:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 56:   PetscCall(MatCreateVecs(A, &b, &x_pristine));
 57:   PetscCall(VecDuplicate(x_pristine, &x_poisoned));
 58:   PetscCall(VecSet(b, 1.0));
 59:   PetscCall(VecSet(x_pristine, 0.0));
 60:   PetscCall(VecSet(x_poisoned, 0.0));

 62:   /* First MatSOR() call: triggers MatInvertDiagonalForSOR_SeqAIJ_Inode() and
 63:      populates a->inode.ibdiag with correctly inverted blocks. */
 64:   PetscCall(MatSOR(A, b, 1.0, SOR_FORWARD_SWEEP, 0.0, 1, 1, x_pristine));

 66:   /* Confirm we actually reached the inode SOR path -- otherwise the test is
 67:      vacuous and would silently pass against a buggy build. */
 68:   a = (Mat_SeqAIJ *)A->data;
 69:   PetscCheck(a->inode.use && a->inode.size_csr && a->inode.node_count > 0 && a->inode.ibdiag != NULL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inode SOR path was not exercised (use=%d node_count=%" PetscInt_FMT " ibdiag=%p); test cannot detect the bug.",
 70:              (int)a->inode.use, a->inode.node_count, (void *)a->inode.ibdiag);

 72:   /* Poison the cached inverted block diagonal.  A is unchanged, so its
 73:      PetscObject state does not advance and the cache *should* be honoured by
 74:      the next MatSOR() call.  A buggy implementation that rebuilds ibdiag on
 75:      every call will overwrite our zeros with the correct inverse and recover
 76:      x_pristine.  A correct implementation will use the zeros and produce a
 77:      materially different x_poisoned. */
 78:   for (PetscInt i = 0; i < a->inode.bdiagsize; i++) a->inode.ibdiag[i] = 0.0;

 80:   PetscCall(MatSOR(A, b, 1.0, SOR_FORWARD_SWEEP, 0.0, 1, 1, x_poisoned));

 82:   PetscCall(VecEqual(x_pristine, x_poisoned, &same));
 83:   PetscCheck(!same, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatSOR_SeqAIJ_Inode() rebuilt its cached block diagonal despite an unchanged matrix state. The cache validity check in MatInvertDiagonalForSOR_SeqAIJ_Inode() is comparing the wrong PetscObjectState field (regression of MR !8797).");

 85:   PetscCall(MatDestroy(&A));
 86:   PetscCall(VecDestroy(&b));
 87:   PetscCall(VecDestroy(&x_pristine));
 88:   PetscCall(VecDestroy(&x_poisoned));
 89:   PetscCall(PetscFinalize());
 90:   return 0;
 91: }

 93: /*TEST

 95:    test:
 96:      suffix: 0
 97:      args: -mat_no_inode false
 98:      output_file: output/empty.out

100: TEST*/