Actual source code: ex241.c
1: static char help[] = "Tests MATHTOOL\n\n";
3: #include <petscmat.h>
5: static PetscErrorCode GenEntries(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt *J,const PetscInt *K,PetscScalar *ptr,void *ctx)
6: {
7: PetscInt d,j,k;
8: PetscReal diff = 0.0,*coords = (PetscReal*)(ctx);
11: for (j = 0; j < M; j++) {
12: for (k = 0; k < N; k++) {
13: diff = 0.0;
14: 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]);
15: ptr[j+M*k] = 1.0/(1.0e-2 + PetscSqrtReal(diff));
16: }
17: }
18: return 0;
19: }
21: static PetscErrorCode GenEntriesRectangular(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt *J,const PetscInt *K,PetscScalar *ptr,void *ctx)
22: {
23: PetscInt d,j,k;
24: PetscReal diff = 0.0,**coords = (PetscReal**)(ctx);
27: for (j = 0; j < M; j++) {
28: for (k = 0; k < N; k++) {
29: diff = 0.0;
30: for (d = 0; d < sdim; d++) diff += (coords[0][J[j]*sdim+d] - coords[1][K[k]*sdim+d]) * (coords[0][J[j]*sdim+d] - coords[1][K[k]*sdim+d]);
31: ptr[j+M*k] = 1.0/(1.0e-2 + PetscSqrtReal(diff));
32: }
33: }
34: return 0;
35: }
37: int main(int argc,char **argv)
38: {
39: Mat A,AT,D,B,P,R,RT;
40: PetscInt m = 100,dim = 3,M,K = 10,begin,n = 0,N,bs;
41: PetscMPIInt size;
42: PetscScalar *ptr;
43: PetscReal *coords,*gcoords,*scoords,*gscoords,*(ctx[2]),norm,epsilon;
44: MatHtoolKernel kernel = GenEntries;
45: PetscBool flg,sym = PETSC_FALSE;
46: PetscRandom rdm;
47: IS iss,ist,is[2];
48: Vec right,left,perm;
50: PetscInitialize(&argc,&argv,(char*)NULL,help);
51: PetscOptionsGetInt(NULL,NULL,"-m_local",&m,NULL);
52: PetscOptionsGetInt(NULL,NULL,"-n_local",&n,NULL);
53: PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);
54: PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL);
55: PetscOptionsGetBool(NULL,NULL,"-symmetric",&sym,NULL);
56: PetscOptionsGetReal(NULL,NULL,"-mat_htool_epsilon",&epsilon,NULL);
57: MPI_Comm_size(PETSC_COMM_WORLD,&size);
58: M = size*m;
59: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
60: PetscMalloc1(m*dim,&coords);
61: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
62: PetscRandomGetValuesReal(rdm,m*dim,coords);
63: PetscCalloc1(M*dim,&gcoords);
64: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,NULL,&B);
65: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
67: MatSetRandom(B,rdm);
68: MatGetOwnershipRange(B,&begin,NULL);
69: PetscArraycpy(gcoords+begin*dim,coords,m*dim);
70: MPIU_Allreduce(MPI_IN_PLACE,gcoords,M*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);
71: MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&A);
72: MatSetOption(A,MAT_SYMMETRIC,sym);
73: MatSetFromOptions(A);
74: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
76: MatViewFromOptions(A,NULL,"-A_view");
77: MatGetOwnershipIS(A,is,NULL);
78: ISDuplicate(is[0],is+1);
79: MatIncreaseOverlap(A,1,is,2);
80: MatSetBlockSize(A,2);
81: MatIncreaseOverlap(A,1,is+1,1);
82: ISGetBlockSize(is[1],&bs);
84: MatSetBlockSize(A,1);
85: ISEqual(is[0],is[1],&flg);
87: ISDestroy(is);
88: ISDestroy(is+1);
89: MatCreateVecs(A,&right,&left);
90: VecSetRandom(right,rdm);
91: MatMult(A,right,left);
92: MatHtoolGetPermutationSource(A,&iss);
93: MatHtoolGetPermutationTarget(A,&ist);
94: VecDuplicate(left,&perm);
95: VecCopy(left,perm);
96: VecPermute(perm,ist,PETSC_FALSE);
97: VecPermute(right,iss,PETSC_FALSE);
98: MatHtoolUsePermutation(A,PETSC_FALSE);
99: MatMult(A,right,left);
100: VecAXPY(left,-1.0,perm);
101: VecNorm(left,NORM_INFINITY,&norm);
103: MatHtoolUsePermutation(A,PETSC_TRUE);
104: VecDestroy(&perm);
105: VecDestroy(&left);
106: VecDestroy(&right);
107: ISDestroy(&ist);
108: ISDestroy(&iss);
109: if (PetscAbsReal(epsilon) >= PETSC_SMALL) { /* when there is compression, it is more difficult to check against MATDENSE, so just compare symmetric and nonsymmetric assemblies */
110: PetscReal relative;
111: MatDestroy(&B);
112: MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&B);
113: MatSetOption(B,MAT_SYMMETRIC,(PetscBool)!sym);
114: MatSetFromOptions(B);
115: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
116: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
117: MatViewFromOptions(B,NULL,"-B_view");
118: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&P);
119: MatNorm(P,NORM_FROBENIUS,&relative);
120: MatConvert(B,MATDENSE,MAT_INITIAL_MATRIX,&R);
121: MatAXPY(R,-1.0,P,SAME_NONZERO_PATTERN);
122: MatNorm(R,NORM_INFINITY,&norm);
124: MatDestroy(&B);
125: MatDestroy(&R);
126: MatDestroy(&P);
127: } else {
128: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&D);
129: MatViewFromOptions(D,NULL,"-D_view");
130: MatMultEqual(A,D,10,&flg);
132: MatMultTransposeEqual(A,D,10,&flg);
134: MatMultAddEqual(A,D,10,&flg);
136: MatGetOwnershipRange(B,&begin,NULL);
137: MatDenseGetArrayWrite(D,&ptr);
138: for (PetscInt i = begin; i < m+begin; ++i)
139: for (PetscInt j = 0; j < M; ++j) GenEntries(dim,1,1,&i,&j,ptr+i-begin+j*m,gcoords);
140: MatDenseRestoreArrayWrite(D,&ptr);
141: MatMultEqual(A,D,10,&flg);
143: MatTranspose(D,MAT_INPLACE_MATRIX,&D);
144: MatTranspose(A,MAT_INITIAL_MATRIX,&AT);
145: MatMultEqual(AT,D,10,&flg);
147: MatTranspose(A,MAT_REUSE_MATRIX,&AT);
148: MatMultEqual(AT,D,10,&flg);
150: MatAXPY(D,-1.0,AT,SAME_NONZERO_PATTERN);
151: MatNorm(D,NORM_INFINITY,&norm);
153: MatDestroy(&AT);
154: MatDestroy(&D);
155: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&P);
156: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
157: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
158: MatMatMultEqual(A,B,P,10,&flg);
160: MatTransposeMatMultEqual(A,B,P,10,&flg);
162: MatDestroy(&B);
163: MatDestroy(&P);
164: if (n) {
165: PetscMalloc1(n*dim,&scoords);
166: PetscRandomGetValuesReal(rdm,n*dim,scoords);
167: N = n;
168: MPIU_Allreduce(MPI_IN_PLACE,&N,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
169: PetscCalloc1(N*dim,&gscoords);
170: MPI_Exscan(&n,&begin,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
171: PetscArraycpy(gscoords+begin*dim,scoords,n*dim);
172: MPIU_Allreduce(MPI_IN_PLACE,gscoords,N*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);
173: kernel = GenEntriesRectangular;
174: ctx[0] = gcoords;
175: ctx[1] = gscoords;
176: MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,n,M,N,dim,coords,scoords,kernel,ctx,&R);
177: MatSetFromOptions(R);
178: MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY);
179: MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY);
180: MatViewFromOptions(R,NULL,"-R_view");
181: MatConvert(R,MATDENSE,MAT_INITIAL_MATRIX,&D);
182: MatViewFromOptions(D,NULL,"-D_view");
183: MatMultEqual(R,D,10,&flg);
185: MatTranspose(D,MAT_INPLACE_MATRIX,&D);
186: MatTranspose(R,MAT_INITIAL_MATRIX,&RT);
187: MatMultEqual(RT,D,10,&flg);
189: MatTranspose(R,MAT_REUSE_MATRIX,&RT);
190: MatMultEqual(RT,D,10,&flg);
192: MatDestroy(&RT);
193: MatDestroy(&D);
194: MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,PETSC_DETERMINE,K,NULL,&B);
195: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
196: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
197: MatSetRandom(B,rdm);
198: MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&P);
199: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
200: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
201: MatMatMultEqual(R,B,P,10,&flg);
203: MatDestroy(&B);
204: MatDestroy(&P);
205: MatCreateVecs(R,&right,&left);
206: VecSetRandom(right,rdm);
207: MatMult(R,right,left);
208: MatHtoolGetPermutationSource(R,&iss);
209: MatHtoolGetPermutationTarget(R,&ist);
210: VecDuplicate(left,&perm);
211: VecCopy(left,perm);
212: VecPermute(perm,ist,PETSC_FALSE);
213: VecPermute(right,iss,PETSC_FALSE);
214: MatHtoolUsePermutation(R,PETSC_FALSE);
215: MatMult(R,right,left);
216: VecAXPY(left,-1.0,perm);
217: VecNorm(left,NORM_INFINITY,&norm);
219: MatHtoolUsePermutation(R,PETSC_TRUE);
220: VecDestroy(&perm);
221: VecDestroy(&left);
222: VecDestroy(&right);
223: ISDestroy(&ist);
224: ISDestroy(&iss);
225: MatDestroy(&R);
226: PetscFree(gscoords);
227: PetscFree(scoords);
228: }
229: }
230: PetscRandomDestroy(&rdm);
231: MatDestroy(&A);
232: PetscFree(gcoords);
233: PetscFree(coords);
234: PetscFinalize();
235: return 0;
236: }
238: /*TEST
240: build:
241: requires: htool
243: test:
244: requires: htool
245: suffix: 1
246: nsize: 4
247: args: -m_local 80 -n_local 25 -mat_htool_epsilon 1.0e-11 -symmetric {{false true}shared output}
248: output_file: output/ex101.out
250: test:
251: requires: htool
252: suffix: 2
253: nsize: 4
254: args: -m_local 120 -mat_htool_epsilon 1.0e-2 -mat_htool_compressor {{sympartialACA fullACA SVD}shared output} -mat_htool_clustering {{PCARegular PCAGeometric BoundingBox1Regular BoundingBox1Geometric}shared output}
255: output_file: output/ex101.out
257: TEST*/