1: #include <petscsys.h>
  2: #include <petscmat.h>
  3: #include <petscoptions.h>

  5: static const char help[] = "Test MatGetValue/Row for hypre matrix on device\n";

  7: int main(int argc, char **argv)
  8: {
  9:   PetscInt  n = 10;   // elements
 10:   PetscReal L = 10.0; // domain length

 12:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 14:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 15:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-L", &L, NULL));

 17:   const PetscInt  N = n + 1; // nodes
 18:   const PetscReal h = L / (PetscReal)n;

 20:   Mat M;
 21:   PetscCall(MatCreate(PETSC_COMM_WORLD, &M));
 22:   PetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, N, N));
 23:   PetscCall(MatSetType(M, MATAIJ));
 24:   PetscCall(MatSetFromOptions(M));
 25:   // Tridiagonal pattern: up to 3 nnz/row (ends have 2). Use same for on/off diag for simplicity.
 26:   PetscCall(MatSetUp(M));
 27:   PetscCall(MatMPIAIJSetPreallocation(M, 3, NULL, 3, NULL));
 28:   PetscCall(MatSeqAIJSetPreallocation(M, 3, NULL));

 30:   // Ownership range for rows (nodes)
 31:   PetscInt rstart, rend;
 32:   PetscCall(MatGetOwnershipRange(M, &rstart, &rend));

 34:   // Element matrix (scaled)
 35:   const PetscScalar s     = (PetscScalar)(h / 6.0);
 36:   const PetscScalar ae[4] = {2 * s, 1 * s, 1 * s, 2 * s};

 38:   // Assemble by looping over elements that start at locally-owned row i
 39:   PetscInt idx[2];
 40:   for (PetscInt i = PetscMax(rstart, 0); i < PetscMin(rend, n); ++i) {
 41:     idx[0] = i;
 42:     idx[1] = i + 1;
 43:     PetscCall(MatSetValues(M, 2, idx, 2, idx, ae, ADD_VALUES));
 44:   }

 46:   PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
 47:   PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));

 49:   // --- Verification of MatGetRow ------------------------------------------
 50:   {
 51:     const PetscReal tol = 100 * PETSC_MACHINE_EPSILON; // tight but safe
 52:     PetscBool       ok  = PETSC_TRUE;

 54:     for (PetscInt i = rstart; i < rend; ++i) {
 55:       const PetscInt    *cols;
 56:       const PetscScalar *getrowvals;
 57:       PetscScalar       *getvals;
 58:       PetscInt           ncols = 0, expected_ncols = 0;
 59:       PetscInt           expected_cols[3];
 60:       PetscScalar        expected_vals[3];

 62:       // Build expected data
 63:       if (i > 0) {
 64:         expected_cols[expected_ncols] = i - 1;
 65:         expected_vals[expected_ncols] = s;
 66:         ++expected_ncols;
 67:       }
 68:       expected_cols[expected_ncols] = i;
 69:       expected_vals[expected_ncols] = (i == 0 || i == N - 1) ? (2 * s) : (4 * s); // diag: h/3 at ends, 2h/3 interior
 70:       ++expected_ncols;
 71:       if (i + 1 < N) {
 72:         expected_cols[expected_ncols] = i + 1;
 73:         expected_vals[expected_ncols] = s;
 74:         ++expected_ncols;
 75:       }

 77:       PetscCall(PetscMalloc1(expected_ncols, &getvals));
 78:       PetscCall(MatGetRow(M, i, &ncols, &cols, &getrowvals));
 79:       PetscCall(MatGetValues(M, 1, &i, expected_ncols, (PetscInt *)expected_cols, getvals));

 81:       // Compare counts
 82:       if (ncols != expected_ncols) {
 83:         ok = PETSC_FALSE;
 84:         goto rowdone;
 85:       }

 87:       // Compare values (match by column)
 88:       for (PetscInt k = 0; k < ncols; ++k) {
 89:         PetscInt expected_k = -1;
 90:         /* Matrix is small. Just do a linear search */
 91:         for (PetscInt l = 0; l < expected_ncols; ++l) {
 92:           if (expected_cols[l] == cols[k]) {
 93:             expected_k = l;
 94:             break;
 95:           }
 96:         }
 97:         if (expected_k < 0) {
 98:           ok = PETSC_FALSE;
 99:           goto rowdone;
100:         }
101:         if ((PetscAbsScalar(getrowvals[k] - expected_vals[expected_k]) > tol) || (PetscAbsScalar(getvals[expected_k] - expected_vals[expected_k]) > tol)) {
102:           ok = PETSC_FALSE;
103:           goto rowdone;
104:         }
105:       }

107:     rowdone:
108:       PetscCall(MatRestoreRow(M, i, &ncols, &cols, &getrowvals));
109:       PetscCall(PetscFree(getvals));
110:       if (!ok) break;
111:     }

113:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ok, 1, MPI_C_BOOL, MPI_LAND, PETSC_COMM_WORLD));
114:     if (ok) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mass matrix check: OK\n"));
115:     else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mass matrix check: FAILED\n"));
116:   }
117:   // --------------------------------------------------------------------------

119:   PetscCall(MatDestroy(&M));
120:   PetscCall(PetscFinalize());
121:   return 0;
122: }

124: /*TEST

126:    test:
127:       requires: hypre
128:       suffix: 1
129:       args: -mat_type hypre

131:    test:
132:       requires: hypre
133:       suffix: 2
134:       args: -mat_type hypre
135:       nsize: 4

137: TEST*/