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); w = PetscRealPart(val); 26 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); 27 } 28 } 29 PetscFunctionReturn(0); 30 } 31 32 int main(int argc,char **args) 33 { 34 Mat A; 35 PetscInt M = 11,N = 13; 36 PetscInt rstart,rend,i,j; 37 PetscViewer view; 38 39 PetscCall(PetscInitialize(&argc,&args,NULL,help)); 40 /* 41 Create a parallel AIJ matrix shared by all processors 42 */ 43 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,M,N,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&A)); 44 45 /* 46 Set values into the matrix 47 */ 48 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 49 for (i=rstart; i<rend; i++) { 50 for (j=0; j<N; j++) { 51 PetscReal v = MakeValue(i,j,M); 52 if (PetscAbsReal(v) > 0) { 53 PetscCall(MatSetValue(A,i,j,v,INSERT_VALUES)); 54 } 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++) { 66 PetscCall(MatView(A,view)); 67 } 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,MATAIJ)); 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 SEQAIJ 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,MATSEQAIJ)); 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 MPIAIJ 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,MATMPIAIJ)); 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/ex44.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_15 137 nsize: 15 138 139 testset: 140 requires: mpiio 141 args: -viewer_binary_mpiio 1 142 output_file: output/ex44.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_15 157 nsize: 15 158 159 TEST*/ 160