1: static char help[] = "Tests MatView()/MatLoad() with binary viewers for BAIJ matrices.\n\n";

  3: #include <petscmat.h>
  4: #include <petscviewer.h>

  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: }

 13: static PetscErrorCode CheckValuesAIJ(Mat A)
 14: {
 15:   PetscInt        M,N,rstart,rend,i,j;
 16:   PetscReal       v,w;
 17:   PetscScalar     val;
 18:   PetscErrorCode  ierr;

 21:   MatGetSize(A,&M,&N);
 22:   MatGetOwnershipRange(A,&rstart,&rend);
 23:   for (i=rstart; i<rend; i++) {
 24:     for (j=0; j<N; j++) {
 25:       MatGetValue(A,i,j,&val);
 26:       v = MakeValue(i,j,M); w = PetscRealPart(val);
 27:       if (PetscAbsReal(v-w) > 0) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Matrix entry (%D,%D) should be %g, got %g",i,j,(double)v,(double)w);
 28:     }
 29:   }
 30:   return(0);
 31: }

 33: int main(int argc,char **args)
 34: {
 35:   Mat            A;
 36:   PetscInt       M = 24,N = 48,bs = 2;
 37:   PetscInt       rstart,rend,i,j;
 39:   PetscViewer    view;

 41:   PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr;
 42:   /*
 43:       Create a parallel BAIJ matrix shared by all processors
 44:   */
 45:   MatCreateBAIJ(PETSC_COMM_WORLD,
 46:                        bs,
 47:                        PETSC_DECIDE,PETSC_DECIDE,
 48:                        M,N,
 49:                        PETSC_DECIDE,NULL,
 50:                        PETSC_DECIDE,NULL,
 51:                        &A);

 53:   /*
 54:       Set values into the matrix
 55:   */
 56:   MatGetSize(A,&M,&N);
 57:   MatGetOwnershipRange(A,&rstart,&rend);
 58:   for (i=rstart; i<rend; i++) {
 59:     for (j=0; j<N; j++) {
 60:       PetscReal v = MakeValue(i,j,M);
 61:       if (PetscAbsReal(v) > 0) {
 62:         MatSetValue(A,i,j,v,INSERT_VALUES);
 63:       }
 64:     }
 65:   }
 66:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 67:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 68:   MatViewFromOptions(A,NULL,"-mat_base_view");

 70:   /*
 71:       Store the binary matrix to a file
 72:   */
 73:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view);
 74:   for (i=0; i<3; i++) {
 75:     MatView(A,view);
 76:   }
 77:   PetscViewerDestroy(&view);
 78:   MatDestroy(&A);

 80:   /*
 81:       Now reload the matrix and check its values
 82:   */
 83:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"matrix.dat",FILE_MODE_READ,&view);
 84:   MatCreate(PETSC_COMM_WORLD,&A);
 85:   MatSetType(A,MATBAIJ);
 86:   for (i=0; i<3; i++) {
 87:     if (i > 0) {MatZeroEntries(A);}
 88:     MatLoad(A,view);
 89:     CheckValuesAIJ(A);
 90:   }
 91:   PetscViewerDestroy(&view);
 92:   MatViewFromOptions(A,NULL,"-mat_load_view");
 93:   MatDestroy(&A);

 95:   /*
 96:       Reload in SEQBAIJ matrix and check its values
 97:   */
 98:   PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_READ,&view);
 99:   MatCreate(PETSC_COMM_SELF,&A);
100:   MatSetType(A,MATSEQBAIJ);
101:   for (i=0; i<3; i++) {
102:     if (i > 0) {MatZeroEntries(A);}
103:     MatLoad(A,view);
104:     CheckValuesAIJ(A);
105:   }
106:   PetscViewerDestroy(&view);
107:   MatDestroy(&A);

109:   /*
110:      Reload in MPIBAIJ matrix and check its values
111:   */
112:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"matrix.dat",FILE_MODE_READ,&view);
113:   MatCreate(PETSC_COMM_WORLD,&A);
114:   MatSetType(A,MATMPIBAIJ);
115:   for (i=0; i<3; i++) {
116:     if (i > 0) {MatZeroEntries(A);}
117:     MatLoad(A,view);
118:     CheckValuesAIJ(A);
119:   }
120:   PetscViewerDestroy(&view);
121:   MatDestroy(&A);

123:   PetscFinalize();
124:   return ierr;
125: }

127: /*TEST

129:    testset:
130:       args: -viewer_binary_mpiio 0
131:       output_file: output/ex45.out
132:       test:
133:         suffix: stdio_1
134:         nsize: 1
135:       test:
136:         suffix: stdio_2
137:         nsize: 2
138:       test:
139:         suffix: stdio_3
140:         nsize: 3
141:       test:
142:         suffix: stdio_4
143:         nsize: 4
144:       test:
145:         suffix: stdio_5
146:         nsize: 4

148:    testset:
149:       requires: mpiio
150:       args: -viewer_binary_mpiio 1
151:       output_file: output/ex45.out
152:       test:
153:         suffix: mpiio_1
154:         nsize: 1
155:       test:
156:         suffix: mpiio_2
157:         nsize: 2
158:       test:
159:         suffix: mpiio_3
160:         nsize: 3
161:       test:
162:         suffix: mpiio_4
163:         nsize: 4
164:       test:
165:         suffix: mpiio_5
166:         nsize: 5

168: TEST*/