1 static char help[] = "Tests MatView()/MatLoad() with binary viewers for AIJ 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 = 11, N = 13; 34 PetscInt rstart, rend, i, j; 35 PetscViewer view; 36 37 PetscFunctionBeginUser; 38 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 39 /* 40 Create a parallel AIJ matrix shared by all processors 41 */ 42 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &A)); 43 44 /* 45 Set values into the matrix 46 */ 47 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 48 for (i = rstart; i < rend; i++) { 49 for (j = 0; j < N; j++) { 50 PetscReal v = MakeValue(i, j, M); 51 if (PetscAbsReal(v) > 0) PetscCall(MatSetValue(A, i, j, v, INSERT_VALUES)); 52 } 53 } 54 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 55 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 56 PetscCall(MatViewFromOptions(A, NULL, "-mat_base_view")); 57 58 /* 59 Store the binary matrix to a file 60 */ 61 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view)); 62 for (i = 0; i < 3; i++) PetscCall(MatView(A, view)); 63 PetscCall(PetscViewerDestroy(&view)); 64 PetscCall(MatDestroy(&A)); 65 66 /* 67 Now reload the matrix and check its values 68 */ 69 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view)); 70 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 71 PetscCall(MatSetType(A, MATAIJ)); 72 for (i = 0; i < 3; i++) { 73 if (i > 0) PetscCall(MatZeroEntries(A)); 74 PetscCall(MatLoad(A, view)); 75 PetscCall(CheckValuesAIJ(A)); 76 } 77 PetscCall(PetscViewerDestroy(&view)); 78 PetscCall(MatViewFromOptions(A, NULL, "-mat_load_view")); 79 PetscCall(MatDestroy(&A)); 80 81 /* 82 Reload in SEQAIJ matrix and check its values 83 */ 84 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "matrix.dat", FILE_MODE_READ, &view)); 85 PetscCall(MatCreate(PETSC_COMM_SELF, &A)); 86 PetscCall(MatSetType(A, MATSEQAIJ)); 87 for (i = 0; i < 3; i++) { 88 if (i > 0) PetscCall(MatZeroEntries(A)); 89 PetscCall(MatLoad(A, view)); 90 PetscCall(CheckValuesAIJ(A)); 91 } 92 PetscCall(PetscViewerDestroy(&view)); 93 PetscCall(MatDestroy(&A)); 94 95 /* 96 Reload in MPIAIJ matrix and check its values 97 */ 98 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view)); 99 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 100 PetscCall(MatSetType(A, MATMPIAIJ)); 101 for (i = 0; i < 3; i++) { 102 if (i > 0) PetscCall(MatZeroEntries(A)); 103 PetscCall(MatLoad(A, view)); 104 PetscCall(CheckValuesAIJ(A)); 105 } 106 PetscCall(PetscViewerDestroy(&view)); 107 PetscCall(MatDestroy(&A)); 108 109 PetscCall(PetscFinalize()); 110 return 0; 111 } 112 113 /*TEST 114 115 testset: 116 args: -viewer_binary_mpiio 0 117 output_file: output/ex44.out 118 test: 119 suffix: stdio_1 120 nsize: 1 121 test: 122 suffix: stdio_2 123 nsize: 2 124 test: 125 suffix: stdio_3 126 nsize: 3 127 test: 128 suffix: stdio_4 129 nsize: 4 130 test: 131 suffix: stdio_15 132 nsize: 15 133 134 testset: 135 requires: mpiio 136 args: -viewer_binary_mpiio 1 137 output_file: output/ex44.out 138 test: 139 suffix: mpiio_1 140 nsize: 1 141 test: 142 suffix: mpiio_2 143 nsize: 2 144 test: 145 suffix: mpiio_3 146 nsize: 3 147 test: 148 suffix: mpiio_4 149 nsize: 4 150 test: 151 suffix: mpiio_15 152 nsize: 15 153 154 TEST*/ 155