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