1 static char help[] = "Tests MatView()/MatLoad() with binary viewers for BAIJ matrices.\n\n"; 2 3 #include <petscmat.h> 4 #include <petscviewer.h> 5 6 #include <petsc/private/hashtable.h> 7 static PetscReal MakeValue(PetscInt i, PetscInt j, PetscInt M) { 8 PetscHash_t h = PetscHashCombine(PetscHashInt(i), PetscHashInt(j)); 9 return (PetscReal)((h % 5 == 0) ? (1 + i + j * M) : 0); 10 } 11 12 static PetscErrorCode CheckValuesAIJ(Mat A) { 13 PetscInt M, N, rstart, rend, i, j; 14 PetscReal v, w; 15 PetscScalar val; 16 17 PetscFunctionBegin; 18 PetscCall(MatGetSize(A, &M, &N)); 19 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 20 for (i = rstart; i < rend; i++) { 21 for (j = 0; j < N; j++) { 22 PetscCall(MatGetValue(A, i, j, &val)); 23 v = MakeValue(i, j, M); 24 w = PetscRealPart(val); 25 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); 26 } 27 } 28 PetscFunctionReturn(0); 29 } 30 31 int main(int argc, char **args) { 32 Mat A; 33 PetscInt M = 24, N = 48, bs = 2; 34 PetscInt rstart, rend, i, j; 35 PetscViewer view; 36 37 PetscFunctionBeginUser; 38 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 39 /* 40 Create a parallel BAIJ matrix shared by all processors 41 */ 42 PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &A)); 43 44 /* 45 Set values into the matrix 46 */ 47 PetscCall(MatGetSize(A, &M, &N)); 48 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 49 for (i = rstart; i < rend; i++) { 50 for (j = 0; j < N; j++) { 51 PetscReal v = MakeValue(i, j, M); 52 if (PetscAbsReal(v) > 0) PetscCall(MatSetValue(A, i, j, v, INSERT_VALUES)); 53 } 54 } 55 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 56 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 57 PetscCall(MatViewFromOptions(A, NULL, "-mat_base_view")); 58 59 /* 60 Store the binary matrix to a file 61 */ 62 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view)); 63 for (i = 0; i < 3; i++) PetscCall(MatView(A, view)); 64 PetscCall(PetscViewerDestroy(&view)); 65 PetscCall(MatDestroy(&A)); 66 67 /* 68 Now reload the matrix and check its values 69 */ 70 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view)); 71 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 72 PetscCall(MatSetType(A, MATBAIJ)); 73 for (i = 0; i < 3; i++) { 74 if (i > 0) PetscCall(MatZeroEntries(A)); 75 PetscCall(MatLoad(A, view)); 76 PetscCall(CheckValuesAIJ(A)); 77 } 78 PetscCall(PetscViewerDestroy(&view)); 79 PetscCall(MatViewFromOptions(A, NULL, "-mat_load_view")); 80 PetscCall(MatDestroy(&A)); 81 82 /* 83 Reload in SEQBAIJ matrix and check its values 84 */ 85 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "matrix.dat", FILE_MODE_READ, &view)); 86 PetscCall(MatCreate(PETSC_COMM_SELF, &A)); 87 PetscCall(MatSetType(A, MATSEQBAIJ)); 88 for (i = 0; i < 3; i++) { 89 if (i > 0) PetscCall(MatZeroEntries(A)); 90 PetscCall(MatLoad(A, view)); 91 PetscCall(CheckValuesAIJ(A)); 92 } 93 PetscCall(PetscViewerDestroy(&view)); 94 PetscCall(MatDestroy(&A)); 95 96 /* 97 Reload in MPIBAIJ matrix and check its values 98 */ 99 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view)); 100 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 101 PetscCall(MatSetType(A, MATMPIBAIJ)); 102 for (i = 0; i < 3; i++) { 103 if (i > 0) PetscCall(MatZeroEntries(A)); 104 PetscCall(MatLoad(A, view)); 105 PetscCall(CheckValuesAIJ(A)); 106 } 107 PetscCall(PetscViewerDestroy(&view)); 108 PetscCall(MatDestroy(&A)); 109 110 PetscCall(PetscFinalize()); 111 return 0; 112 } 113 114 /*TEST 115 116 testset: 117 args: -viewer_binary_mpiio 0 118 output_file: output/ex45.out 119 test: 120 suffix: stdio_1 121 nsize: 1 122 test: 123 suffix: stdio_2 124 nsize: 2 125 test: 126 suffix: stdio_3 127 nsize: 3 128 test: 129 suffix: stdio_4 130 nsize: 4 131 test: 132 suffix: stdio_5 133 nsize: 4 134 135 testset: 136 requires: mpiio 137 args: -viewer_binary_mpiio 1 138 output_file: output/ex45.out 139 test: 140 suffix: mpiio_1 141 nsize: 1 142 test: 143 suffix: mpiio_2 144 nsize: 2 145 test: 146 suffix: mpiio_3 147 nsize: 3 148 test: 149 suffix: mpiio_4 150 nsize: 4 151 test: 152 suffix: mpiio_5 153 nsize: 5 154 155 TEST*/ 156