1: static char help[] = "Tests MatView()/MatLoad() with binary viewers for SBAIJ 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:   PetscBool       seqsbaij,mpisbaij,sbaij;
 19:   PetscErrorCode  ierr;

 22:   MatGetSize(A,&M,&N);
 23:   MatGetOwnershipRange(A,&rstart,&rend);
 24:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&seqsbaij);
 25:   PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&mpisbaij);
 26:   sbaij = (seqsbaij||mpisbaij) ? PETSC_TRUE : PETSC_FALSE;
 27:   for (i=rstart; i<rend; i++) {
 28:     for (j=(sbaij?i:0); j<N; j++) {
 29:       MatGetValue(A,i,j,&val);
 30:       v = MakeValue(i,j,M); w = PetscRealPart(val);
 31:       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);
 32:     }
 33:   }
 34:   return(0);
 35: }

 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;
 43:   PetscViewer    view;

 45:   PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr;
 46:   /*
 47:       Create a parallel SBAIJ matrix shared by all processors
 48:   */
 49:   MatCreateSBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,M,N,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&A);

 51:   /*
 52:       Set values into the matrix
 53:   */
 54:   MatGetSize(A,&M,&N);
 55:   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) {
 60:         MatSetValue(A,i,j,v,INSERT_VALUES);
 61:       }
 62:     }
 63:   }
 64:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 65:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 66:   MatViewFromOptions(A,NULL,"-mat_base_view");

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

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

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

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

121:   PetscFinalize();
122:   return ierr;
123: }

125: /*TEST

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

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

166: TEST*/