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 { 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 PetscBool seqsbaij, mpisbaij, sbaij; 19 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 } 36 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; 43 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)); 50 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")); 65 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)); 73 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)); 88 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)); 102 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)); 116 117 PetscCall(PetscFinalize()); 118 return 0; 119 } 120 121 /*TEST 122 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 141 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 161 162 TEST*/ 163