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;
 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;
 48:   Vec            right,left,perm;

 51:   PetscInitialize(&argc,&argv,(char*)NULL,help);if (ierr) return ierr;
 52:   PetscOptionsGetInt(NULL,NULL,"-m_local",&m,NULL);
 53:   PetscOptionsGetInt(NULL,NULL,"-n_local",&n,NULL);
 54:   PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);
 55:   PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL);
 56:   PetscOptionsGetBool(NULL,NULL,"-symmetric",&sym,NULL);
 57:   PetscOptionsGetReal(NULL,NULL,"-mat_htool_epsilon",&epsilon,NULL);
 58:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 59:   M = size*m;
 60:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 61:   PetscMalloc1(m*dim,&coords);
 62:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
 63:   PetscRandomGetValuesReal(rdm,m*dim,coords);
 64:   PetscCalloc1(M*dim,&gcoords);
 65:   MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,NULL,&B);
 66:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 67:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 68:   MatSetRandom(B,rdm);
 69:   MatGetOwnershipRange(B,&begin,NULL);
 70:   PetscArraycpy(gcoords+begin*dim,coords,m*dim);
 71:   MPIU_Allreduce(MPI_IN_PLACE,gcoords,M*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);
 72:   MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&A);
 73:   MatSetOption(A,MAT_SYMMETRIC,sym);
 74:   MatSetFromOptions(A);
 75:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 76:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 77:   MatViewFromOptions(A,NULL,"-A_view");
 78:   MatCreateVecs(A,&right,&left);
 79:   VecSetRandom(right,rdm);
 80:   MatMult(A,right,left);
 81:   MatHtoolGetPermutationSource(A,&iss);
 82:   MatHtoolGetPermutationTarget(A,&ist);
 83:   VecDuplicate(left,&perm);
 84:   VecCopy(left,perm);
 85:   VecPermute(perm,ist,PETSC_FALSE);
 86:   VecPermute(right,iss,PETSC_FALSE);
 87:   MatHtoolUsePermutation(A,PETSC_FALSE);
 88:   MatMult(A,right,left);
 89:   VecAXPY(left,-1.0,perm);
 90:   VecNorm(left,NORM_INFINITY,&norm);
 91:   if (PetscAbsReal(norm) > PETSC_SMALL) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||y(with permutation)-y(without permutation)|| = %g (> %g)",(double)PetscAbsReal(norm),(double)PETSC_SMALL);
 92:   MatHtoolUsePermutation(A,PETSC_TRUE);
 93:   VecDestroy(&perm);
 94:   VecDestroy(&left);
 95:   VecDestroy(&right);
 96:   ISDestroy(&ist);
 97:   ISDestroy(&iss);
 98:   if (PetscAbsReal(epsilon) >= PETSC_SMALL) { /* when there is compression, it is more difficult to check against MATDENSE, so just compare symmetric and nonsymmetric assemblies */
 99:     PetscReal relative;
100:     MatDestroy(&B);
101:     MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&B);
102:     MatSetOption(B,MAT_SYMMETRIC,(PetscBool)!sym);
103:     MatSetFromOptions(B);
104:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
105:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
106:     MatViewFromOptions(B,NULL,"-B_view");
107:     MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&P);
108:     MatNorm(P,NORM_FROBENIUS,&relative);
109:     MatConvert(B,MATDENSE,MAT_INITIAL_MATRIX,&R);
110:     MatAXPY(R,-1.0,P,SAME_NONZERO_PATTERN);
111:     MatNorm(R,NORM_INFINITY,&norm);
112:     if (PetscAbsReal(norm/relative) > epsilon) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||A(!symmetric)-A(symmetric)|| = %g (> %g)",(double)PetscAbsReal(norm/relative),(double)epsilon);
113:     MatDestroy(&B);
114:     MatDestroy(&R);
115:     MatDestroy(&P);
116:   } else {
117:     MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&D);
118:     MatViewFromOptions(D,NULL,"-D_view");
119:     MatMultEqual(A,D,10,&flg);
120:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Ax != Dx");
121:     MatMultAddEqual(A,D,10,&flg);
122:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"y+Ax != y+Dx");
123:     MatGetOwnershipRange(B,&begin,NULL);
124:     MatDenseGetArrayWrite(D,&ptr);
125:     for (PetscInt i = begin; i < m+begin; ++i)
126:       for (PetscInt j = 0; j < M; ++j) GenEntries(dim,1,1,&i,&j,ptr+i-begin+j*m,gcoords);
127:     MatDenseRestoreArrayWrite(D,&ptr);
128:     MatMultEqual(A,D,10,&flg);
129:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Ax != Dx");
130:     MatTranspose(D,MAT_INPLACE_MATRIX,&D);
131:     MatTranspose(A,MAT_INITIAL_MATRIX,&AT);
132:     MatMultEqual(AT,D,10,&flg);
133:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^Tx != D^Tx");
134:     MatTranspose(A,MAT_REUSE_MATRIX,&AT);
135:     MatMultEqual(AT,D,10,&flg);
136:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^Tx != D^Tx");
137:     MatAXPY(D,-1.0,AT,SAME_NONZERO_PATTERN);
138:     MatNorm(D,NORM_INFINITY,&norm);
139:     if (PetscAbsReal(norm) > PETSC_SMALL) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||A-D|| = %g (> %g)",(double)norm,(double)PETSC_SMALL);
140:     MatDestroy(&AT);
141:     MatDestroy(&D);
142:     MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&P);
143:     MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
144:     MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
145:     MatMatMultEqual(A,B,P,10,&flg);
146:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ABx != Px");
147:     MatDestroy(&B);
148:     MatDestroy(&P);
149:     if (n) {
150:       PetscMalloc1(n*dim,&scoords);
151:       PetscRandomGetValuesReal(rdm,n*dim,scoords);
152:       N = n;
153:       MPIU_Allreduce(MPI_IN_PLACE,&N,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
154:       PetscCalloc1(N*dim,&gscoords);
155:       MPI_Exscan(&n,&begin,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
156:       PetscArraycpy(gscoords+begin*dim,scoords,n*dim);
157:       MPIU_Allreduce(MPI_IN_PLACE,gscoords,N*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);
158:       kernel = GenEntriesRectangular;
159:       ctx[0] = gcoords;
160:       ctx[1] = gscoords;
161:       MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,n,M,N,dim,coords,scoords,kernel,ctx,&R);
162:       MatSetFromOptions(R);
163:       MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY);
164:       MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY);
165:       MatViewFromOptions(R,NULL,"-R_view");
166:       MatConvert(R,MATDENSE,MAT_INITIAL_MATRIX,&D);
167:       MatViewFromOptions(D,NULL,"-D_view");
168:       MatMultEqual(R,D,10,&flg);
169:       if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Rx != Dx");
170:       MatTranspose(D,MAT_INPLACE_MATRIX,&D);
171:       MatTranspose(R,MAT_INITIAL_MATRIX,&RT);
172:       MatMultEqual(RT,D,10,&flg);
173:       if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"R^Tx != D^Tx");
174:       MatTranspose(R,MAT_REUSE_MATRIX,&RT);
175:       MatMultEqual(RT,D,10,&flg);
176:       if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"R^Tx != D^Tx");
177:       MatDestroy(&RT);
178:       MatDestroy(&D);
179:       MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,PETSC_DETERMINE,K,NULL,&B);
180:       MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
181:       MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
182:       MatSetRandom(B,rdm);
183:       MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&P);
184:       MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
185:       MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
186:       MatMatMultEqual(R,B,P,10,&flg);
187:       if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"RBx != Px");
188:       MatDestroy(&B);
189:       MatDestroy(&P);
190:       MatCreateVecs(R,&right,&left);
191:       VecSetRandom(right,rdm);
192:       MatMult(R,right,left);
193:       MatHtoolGetPermutationSource(R,&iss);
194:       MatHtoolGetPermutationTarget(R,&ist);
195:       VecDuplicate(left,&perm);
196:       VecCopy(left,perm);
197:       VecPermute(perm,ist,PETSC_FALSE);
198:       VecPermute(right,iss,PETSC_FALSE);
199:       MatHtoolUsePermutation(R,PETSC_FALSE);
200:       MatMult(R,right,left);
201:       VecAXPY(left,-1.0,perm);
202:       VecNorm(left,NORM_INFINITY,&norm);
203:       if (PetscAbsReal(norm) > PETSC_SMALL) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||y(with permutation)-y(without permutation)|| = %g (> %g)",(double)PetscAbsReal(norm),(double)PETSC_SMALL);
204:       MatHtoolUsePermutation(R,PETSC_TRUE);
205:       VecDestroy(&perm);
206:       VecDestroy(&left);
207:       VecDestroy(&right);
208:       ISDestroy(&ist);
209:       ISDestroy(&iss);
210:       MatDestroy(&R);
211:       PetscFree(gscoords);
212:       PetscFree(scoords);
213:     }
214:   }
215:   PetscRandomDestroy(&rdm);
216:   MatDestroy(&A);
217:   PetscFree(gcoords);
218:   PetscFree(coords);
219:   PetscFinalize();
220:   return ierr;
221: }

223: /*TEST

225:    build:
226:       requires: htool

228:    test:
229:       requires: htool
230:       suffix: 1
231:       nsize: 4
232:       args: -m_local 80 -n_local 25 -mat_htool_epsilon 1.0e-11 -symmetric {{false true}shared output}
233:       output_file: output/ex101.out

235:    test:
236:       requires: htool
237:       suffix: 2
238:       nsize: 4
239:       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}
240:       output_file: output/ex101.out

242: TEST*/