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*/