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