Actual source code: ex246.cxx
1: static char help[] = "Tests MATHTOOL with a derived htool::IMatrix<PetscScalar> class\n\n";
3: #include <petscmat.h>
4: #include <htool/misc/petsc.hpp>
6: static PetscErrorCode GenEntries(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt *J,const PetscInt *K,PetscScalar *ptr,void *ctx)
7: {
8: PetscInt d,j,k;
9: PetscReal diff = 0.0,*coords = (PetscReal*)(ctx);
12: for (j = 0; j < M; j++) {
13: for (k = 0; k < N; k++) {
14: diff = 0.0;
15: for (d = 0; d < sdim; d++) diff += (coords[J[j]*sdim+d] - coords[K[k]*sdim+d]) * (coords[J[j]*sdim+d] - coords[K[k]*sdim+d]);
16: ptr[j+M*k] = 1.0/(1.0e-2 + PetscSqrtReal(diff));
17: }
18: }
19: return 0;
20: }
21: class MyIMatrix : public htool::VirtualGenerator<PetscScalar> {
22: private:
23: PetscReal *coords;
24: PetscInt sdim;
25: public:
26: MyIMatrix(PetscInt M,PetscInt N,PetscInt spacedim,PetscReal* gcoords) : htool::VirtualGenerator<PetscScalar>(M,N),coords(gcoords),sdim(spacedim) { }
28: void copy_submatrix(PetscInt M,PetscInt N,const PetscInt *J,const PetscInt *K,PetscScalar *ptr) const override
29: {
30: PetscReal diff = 0.0;
33: for (PetscInt j = 0; j < M; j++) /* could be optimized by the user how they see fit, e.g., vectorization */
34: for (PetscInt k = 0; k < N; k++) {
35: diff = 0.0;
36: for (PetscInt d = 0; d < sdim; d++) diff += (coords[J[j]*sdim+d] - coords[K[k]*sdim+d]) * (coords[J[j]*sdim+d] - coords[K[k]*sdim+d]);
37: ptr[j+M*k] = 1.0/(1.0e-2 + PetscSqrtReal(diff));
38: }
39: return;
40: }
41: };
43: int main(int argc,char **argv)
44: {
45: Mat A,B,P,R;
46: PetscInt m = 100,dim = 3,M,begin = 0;
47: PetscMPIInt size;
48: PetscReal *coords,*gcoords,norm,epsilon,relative;
49: PetscBool sym = PETSC_FALSE;
50: PetscRandom rdm;
51: MatHtoolKernel kernel = GenEntries;
52: MyIMatrix *imatrix;
54: PetscInitialize(&argc,&argv,(char*)NULL,help);
55: PetscOptionsGetInt(NULL,NULL,"-m_local",&m,NULL);
56: PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);
57: PetscOptionsGetBool(NULL,NULL,"-symmetric",&sym,NULL);
58: PetscOptionsGetReal(NULL,NULL,"-mat_htool_epsilon",&epsilon,NULL);
59: MPI_Comm_size(PETSC_COMM_WORLD,&size);
60: M = size*m;
61: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
62: PetscMalloc1(m*dim,&coords);
63: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
64: PetscRandomGetValuesReal(rdm,m*dim,coords);
65: PetscCalloc1(M*dim,&gcoords);
66: MPI_Exscan(&m,&begin,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
67: PetscArraycpy(gcoords+begin*dim,coords,m*dim);
68: MPIU_Allreduce(MPI_IN_PLACE,gcoords,M*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);
69: imatrix = new MyIMatrix(M,M,dim,gcoords);
70: MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,NULL,imatrix,&A); /* block-wise assembly using htool::IMatrix<PetscScalar>::copy_submatrix() */
71: MatSetOption(A,MAT_SYMMETRIC,sym);
72: MatSetFromOptions(A);
73: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
74: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
75: MatViewFromOptions(A,NULL,"-A_view");
76: MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&B); /* entry-wise assembly using GenEntries() */
77: MatSetOption(B,MAT_SYMMETRIC,sym);
78: MatSetFromOptions(B);
79: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
81: MatViewFromOptions(B,NULL,"-B_view");
82: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&P);
83: MatNorm(P,NORM_FROBENIUS,&relative);
84: MatConvert(B,MATDENSE,MAT_INITIAL_MATRIX,&R);
85: MatAXPY(R,-1.0,P,SAME_NONZERO_PATTERN);
86: MatNorm(R,NORM_INFINITY,&norm);
88: MatDestroy(&B);
89: MatDestroy(&R);
90: MatDestroy(&P);
91: PetscRandomDestroy(&rdm);
92: MatDestroy(&A);
93: PetscFree(gcoords);
94: PetscFree(coords);
95: delete imatrix;
96: PetscFinalize();
97: return 0;
98: }
100: /*TEST
102: build:
103: requires: htool
104: test:
105: requires: htool
106: suffix: 2
107: nsize: 4
108: args: -m_local 120 -mat_htool_epsilon 1.0e-2 -symmetric {{false true}shared output}
109: output_file: output/ex101.out
111: TEST*/