1: static char help[] = "Tests MatView()/MatLoad() with binary viewers for SBAIJ matrices.\n\n";

  3: #include <petscmat.h>

  5: #include <petsc/private/hashtable.h>
  6: static PetscReal MakeValue(PetscInt i, PetscInt j, PetscInt M)
  7: {
  8:   PetscHash_t h = PetscHashCombine(PetscHashInt(i), PetscHashInt(j));
  9:   return (PetscReal)((h % 5 == 0) ? (1 + i + j * M) : 0);
 10: }

 12: static PetscErrorCode CheckValuesAIJ(Mat A)
 13: {
 14:   PetscInt    M, N, rstart, rend, i, j;
 15:   PetscReal   v, w;
 16:   PetscScalar val;
 17:   PetscBool   seqsbaij, mpisbaij, sbaij;

 19:   PetscFunctionBegin;
 20:   PetscCall(MatGetSize(A, &M, &N));
 21:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 22:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &seqsbaij));
 23:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISBAIJ, &mpisbaij));
 24:   sbaij = (seqsbaij || mpisbaij) ? PETSC_TRUE : PETSC_FALSE;
 25:   for (i = rstart; i < rend; i++) {
 26:     for (j = (sbaij ? i : 0); j < N; j++) {
 27:       PetscCall(MatGetValue(A, i, j, &val));
 28:       v = MakeValue(i, j, M);
 29:       w = PetscRealPart(val);
 30:       PetscCheck(PetscAbsReal(v - w) <= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Matrix entry (%" PetscInt_FMT ",%" PetscInt_FMT ") should be %g, got %g", i, j, (double)v, (double)w);
 31:     }
 32:   }
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 34: }

 36: int main(int argc, char **args)
 37: {
 38:   Mat         A;
 39:   PetscInt    M = 24, N = 24, bs = 3;
 40:   PetscInt    rstart, rend, i, j;
 41:   PetscViewer view;

 43:   PetscFunctionBeginUser;
 44:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 45:   /*
 46:       Create a parallel SBAIJ matrix shared by all processors
 47:   */
 48:   PetscCall(MatCreateSBAIJ(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &A));

 50:   /*
 51:       Set values into the matrix
 52:   */
 53:   PetscCall(MatGetSize(A, &M, &N));
 54:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 55:   for (i = rstart; i < rend; i++) {
 56:     for (j = 0; j < N; j++) {
 57:       PetscReal v = MakeValue(i, j, M);
 58:       if (PetscAbsReal(v) > 0) PetscCall(MatSetValue(A, i, j, v, INSERT_VALUES));
 59:     }
 60:   }
 61:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 62:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 63:   PetscCall(MatViewFromOptions(A, NULL, "-mat_base_view"));

 65:   /*
 66:       Store the binary matrix to a file
 67:   */
 68:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view));
 69:   for (i = 0; i < 3; i++) PetscCall(MatView(A, view));
 70:   PetscCall(PetscViewerDestroy(&view));
 71:   PetscCall(MatDestroy(&A));

 73:   /*
 74:       Now reload the matrix and check its values
 75:   */
 76:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
 77:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 78:   PetscCall(MatSetType(A, MATSBAIJ));
 79:   for (i = 0; i < 3; i++) {
 80:     if (i > 0) PetscCall(MatZeroEntries(A));
 81:     PetscCall(MatLoad(A, view));
 82:     PetscCall(CheckValuesAIJ(A));
 83:   }
 84:   PetscCall(PetscViewerDestroy(&view));
 85:   PetscCall(MatViewFromOptions(A, NULL, "-mat_load_view"));
 86:   PetscCall(MatDestroy(&A));

 88:   /*
 89:       Reload in SEQSBAIJ matrix and check its values
 90:   */
 91:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "matrix.dat", FILE_MODE_READ, &view));
 92:   PetscCall(MatCreate(PETSC_COMM_SELF, &A));
 93:   PetscCall(MatSetType(A, MATSEQSBAIJ));
 94:   for (i = 0; i < 3; i++) {
 95:     if (i > 0) PetscCall(MatZeroEntries(A));
 96:     PetscCall(MatLoad(A, view));
 97:     PetscCall(CheckValuesAIJ(A));
 98:   }
 99:   PetscCall(PetscViewerDestroy(&view));
100:   PetscCall(MatDestroy(&A));

102:   /*
103:      Reload in MPISBAIJ matrix and check its values
104:   */
105:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
106:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
107:   PetscCall(MatSetType(A, MATMPISBAIJ));
108:   for (i = 0; i < 3; i++) {
109:     if (i > 0) PetscCall(MatZeroEntries(A));
110:     PetscCall(MatLoad(A, view));
111:     PetscCall(CheckValuesAIJ(A));
112:   }
113:   PetscCall(PetscViewerDestroy(&view));
114:   PetscCall(MatDestroy(&A));

116:   PetscCall(PetscFinalize());
117:   return 0;
118: }

120: /*TEST

122:    testset:
123:       args: -viewer_binary_mpiio 0
124:       output_file: output/empty.out
125:       test:
126:         suffix: stdio_1
127:         nsize: 1
128:       test:
129:         suffix: stdio_2
130:         nsize: 2
131:       test:
132:         suffix: stdio_3
133:         nsize: 3
134:       test:
135:         suffix: stdio_4
136:         nsize: 4
137:       test:
138:         suffix: stdio_5
139:         nsize: 4

141:    testset:
142:       requires: mpiio
143:       args: -viewer_binary_mpiio 1
144:       output_file: output/empty.out
145:       test:
146:         suffix: mpiio_1
147:         nsize: 1
148:       test:
149:         suffix: mpiio_2
150:         nsize: 2
151:       test:
152:         suffix: mpiio_3
153:         nsize: 3
154:       test:
155:         suffix: mpiio_4
156:         nsize: 4
157:       test:
158:         suffix: mpiio_5
159:         nsize: 5

161: TEST*/