1: static char help[] = "Tests MatView()/MatLoad() with binary viewers for AIJ 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 = 11, N = 13;
 36:   PetscInt    rstart, rend, i, j;
 37:   PetscViewer view;

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

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

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

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

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

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

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

115: /*TEST

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

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

156: TEST*/