Actual source code: ex50.c

  1: static char help[] = "Tests MatView()/MatLoad() with binary viewers for SBAIJ 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;
 18:   PetscBool   seqsbaij, mpisbaij, sbaij;

 20:   PetscFunctionBegin;
 21:   PetscCall(MatGetSize(A, &M, &N));
 22:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 23:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &seqsbaij));
 24:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISBAIJ, &mpisbaij));
 25:   sbaij = (seqsbaij || mpisbaij) ? PETSC_TRUE : PETSC_FALSE;
 26:   for (i = rstart; i < rend; i++) {
 27:     for (j = (sbaij ? i : 0); j < N; j++) {
 28:       PetscCall(MatGetValue(A, i, j, &val));
 29:       v = MakeValue(i, j, M);
 30:       w = PetscRealPart(val);
 31:       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);
 32:     }
 33:   }
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 35: }

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

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

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

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

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

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

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

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

121: /*TEST

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

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

162: TEST*/