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