1: static char help[] = "Tests MatView()/MatLoad() with binary viewers for BAIJ 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;

 18:   PetscFunctionBegin;
 19:   PetscCall(MatGetSize(A, &M, &N));
 20:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 21:   for (i = rstart; i < rend; i++) {
 22:     for (j = 0; j < N; j++) {
 23:       PetscCall(MatGetValue(A, i, j, &val));
 24:       v = MakeValue(i, j, M);
 25:       w = PetscRealPart(val);
 26:       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);
 27:     }
 28:   }
 29:   PetscFunctionReturn(PETSC_SUCCESS);
 30: }

 32: int main(int argc, char **args)
 33: {
 34:   Mat         A;
 35:   PetscInt    M = 24, N = 48, bs = 2;
 36:   PetscInt    rstart, rend, i, j;
 37:   PetscViewer view;

 39:   PetscFunctionBeginUser;
 40:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 41:   /*
 42:       Create a parallel BAIJ matrix shared by all processors
 43:   */
 44:   PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &A));

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

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

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

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

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

112:   PetscCall(PetscFinalize());
113:   return 0;
114: }

116: /*TEST

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

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

157: TEST*/