Actual source code: ex45.c

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

  3: #include <petscmat.h>
  4: #include <petscviewer.h>

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

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

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

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

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

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

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

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

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

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

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

117: /*TEST

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

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

158: TEST*/