1 static char help[] = "Tests MATHTOOL\n\n"; 2 3 #include <petscmat.h> 4 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); 9 10 PetscFunctionBeginUser; 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 PetscFunctionReturn(0); 19 } 20 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); 25 26 PetscFunctionBeginUser; 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 PetscFunctionReturn(0); 35 } 36 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; 49 PetscErrorCode ierr; 50 51 ierr = PetscInitialize(&argc,&argv,(char*)NULL,help);if (ierr) return ierr; 52 ierr = PetscOptionsGetInt(NULL,NULL,"-m_local",&m,NULL);CHKERRQ(ierr); 53 ierr = PetscOptionsGetInt(NULL,NULL,"-n_local",&n,NULL);CHKERRQ(ierr); 54 ierr = PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);CHKERRQ(ierr); 55 ierr = PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL);CHKERRQ(ierr); 56 ierr = PetscOptionsGetBool(NULL,NULL,"-symmetric",&sym,NULL);CHKERRQ(ierr); 57 ierr = PetscOptionsGetReal(NULL,NULL,"-mat_htool_epsilon",&epsilon,NULL);CHKERRQ(ierr); 58 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 59 M = size*m; 60 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 61 ierr = PetscMalloc1(m*dim,&coords);CHKERRQ(ierr); 62 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr); 63 ierr = PetscRandomGetValuesReal(rdm,m*dim,coords);CHKERRQ(ierr); 64 ierr = PetscCalloc1(M*dim,&gcoords);CHKERRQ(ierr); 65 ierr = MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,NULL,&B);CHKERRQ(ierr); 66 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 67 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 68 ierr = MatSetRandom(B,rdm);CHKERRQ(ierr); 69 ierr = MatGetOwnershipRange(B,&begin,NULL);CHKERRQ(ierr); 70 ierr = PetscArraycpy(gcoords+begin*dim,coords,m*dim);CHKERRQ(ierr); 71 ierr = MPIU_Allreduce(MPI_IN_PLACE,gcoords,M*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);CHKERRMPI(ierr); 72 ierr = MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&A);CHKERRQ(ierr); 73 ierr = MatSetOption(A,MAT_SYMMETRIC,sym);CHKERRQ(ierr); 74 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 75 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 76 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 77 ierr = MatViewFromOptions(A,NULL,"-A_view");CHKERRQ(ierr); 78 ierr = MatGetOwnershipIS(A,is,NULL);CHKERRQ(ierr); 79 ierr = ISDuplicate(is[0],is+1);CHKERRQ(ierr); 80 ierr = MatIncreaseOverlap(A,1,is,2);CHKERRQ(ierr); 81 ierr = MatSetBlockSize(A,2);CHKERRQ(ierr); 82 ierr = MatIncreaseOverlap(A,1,is+1,1);CHKERRQ(ierr); 83 ierr = ISGetBlockSize(is[1],&bs);CHKERRQ(ierr); 84 PetscCheckFalse(bs != 2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrect block size %" PetscInt_FMT " != 2",bs); 85 ierr = MatSetBlockSize(A,1);CHKERRQ(ierr); 86 ierr = ISEqual(is[0],is[1],&flg);CHKERRQ(ierr); 87 PetscCheckFalse(!flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unequal index sets"); 88 ierr = ISDestroy(is);CHKERRQ(ierr); 89 ierr = ISDestroy(is+1);CHKERRQ(ierr); 90 ierr = MatCreateVecs(A,&right,&left);CHKERRQ(ierr); 91 ierr = VecSetRandom(right,rdm);CHKERRQ(ierr); 92 ierr = MatMult(A,right,left);CHKERRQ(ierr); 93 ierr = MatHtoolGetPermutationSource(A,&iss);CHKERRQ(ierr); 94 ierr = MatHtoolGetPermutationTarget(A,&ist);CHKERRQ(ierr); 95 ierr = VecDuplicate(left,&perm);CHKERRQ(ierr); 96 ierr = VecCopy(left,perm);CHKERRQ(ierr); 97 ierr = VecPermute(perm,ist,PETSC_FALSE);CHKERRQ(ierr); 98 ierr = VecPermute(right,iss,PETSC_FALSE);CHKERRQ(ierr); 99 ierr = MatHtoolUsePermutation(A,PETSC_FALSE);CHKERRQ(ierr); 100 ierr = MatMult(A,right,left);CHKERRQ(ierr); 101 ierr = VecAXPY(left,-1.0,perm);CHKERRQ(ierr); 102 ierr = VecNorm(left,NORM_INFINITY,&norm);CHKERRQ(ierr); 103 PetscCheckFalse(PetscAbsReal(norm) > PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||y(with permutation)-y(without permutation)|| = %g (> %g)",(double)PetscAbsReal(norm),(double)PETSC_SMALL); 104 ierr = MatHtoolUsePermutation(A,PETSC_TRUE);CHKERRQ(ierr); 105 ierr = VecDestroy(&perm);CHKERRQ(ierr); 106 ierr = VecDestroy(&left);CHKERRQ(ierr); 107 ierr = VecDestroy(&right);CHKERRQ(ierr); 108 ierr = ISDestroy(&ist);CHKERRQ(ierr); 109 ierr = ISDestroy(&iss);CHKERRQ(ierr); 110 if (PetscAbsReal(epsilon) >= PETSC_SMALL) { /* when there is compression, it is more difficult to check against MATDENSE, so just compare symmetric and nonsymmetric assemblies */ 111 PetscReal relative; 112 ierr = MatDestroy(&B);CHKERRQ(ierr); 113 ierr = MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&B);CHKERRQ(ierr); 114 ierr = MatSetOption(B,MAT_SYMMETRIC,(PetscBool)!sym);CHKERRQ(ierr); 115 ierr = MatSetFromOptions(B);CHKERRQ(ierr); 116 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 117 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 118 ierr = MatViewFromOptions(B,NULL,"-B_view");CHKERRQ(ierr); 119 ierr = MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&P);CHKERRQ(ierr); 120 ierr = MatNorm(P,NORM_FROBENIUS,&relative);CHKERRQ(ierr); 121 ierr = MatConvert(B,MATDENSE,MAT_INITIAL_MATRIX,&R);CHKERRQ(ierr); 122 ierr = MatAXPY(R,-1.0,P,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 123 ierr = MatNorm(R,NORM_INFINITY,&norm);CHKERRQ(ierr); 124 PetscCheckFalse(PetscAbsReal(norm/relative) > epsilon,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||A(!symmetric)-A(symmetric)|| = %g (> %g)",(double)PetscAbsReal(norm/relative),(double)epsilon); 125 ierr = MatDestroy(&B);CHKERRQ(ierr); 126 ierr = MatDestroy(&R);CHKERRQ(ierr); 127 ierr = MatDestroy(&P);CHKERRQ(ierr); 128 } else { 129 ierr = MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&D);CHKERRQ(ierr); 130 ierr = MatViewFromOptions(D,NULL,"-D_view");CHKERRQ(ierr); 131 ierr = MatMultEqual(A,D,10,&flg);CHKERRQ(ierr); 132 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Ax != Dx"); 133 ierr = MatMultTransposeEqual(A,D,10,&flg);CHKERRQ(ierr); 134 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^Tx != D^Tx"); 135 ierr = MatMultAddEqual(A,D,10,&flg);CHKERRQ(ierr); 136 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"y+Ax != y+Dx"); 137 ierr = MatGetOwnershipRange(B,&begin,NULL);CHKERRQ(ierr); 138 ierr = MatDenseGetArrayWrite(D,&ptr);CHKERRQ(ierr); 139 for (PetscInt i = begin; i < m+begin; ++i) 140 for (PetscInt j = 0; j < M; ++j) GenEntries(dim,1,1,&i,&j,ptr+i-begin+j*m,gcoords); 141 ierr = MatDenseRestoreArrayWrite(D,&ptr);CHKERRQ(ierr); 142 ierr = MatMultEqual(A,D,10,&flg);CHKERRQ(ierr); 143 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Ax != Dx"); 144 ierr = MatTranspose(D,MAT_INPLACE_MATRIX,&D);CHKERRQ(ierr); 145 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&AT);CHKERRQ(ierr); 146 ierr = MatMultEqual(AT,D,10,&flg);CHKERRQ(ierr); 147 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^Tx != D^Tx"); 148 ierr = MatTranspose(A,MAT_REUSE_MATRIX,&AT);CHKERRQ(ierr); 149 ierr = MatMultEqual(AT,D,10,&flg);CHKERRQ(ierr); 150 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^Tx != D^Tx"); 151 ierr = MatAXPY(D,-1.0,AT,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 152 ierr = MatNorm(D,NORM_INFINITY,&norm);CHKERRQ(ierr); 153 PetscCheckFalse(PetscAbsReal(norm) > PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||A-D|| = %g (> %g)",(double)norm,(double)PETSC_SMALL); 154 ierr = MatDestroy(&AT);CHKERRQ(ierr); 155 ierr = MatDestroy(&D);CHKERRQ(ierr); 156 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&P);CHKERRQ(ierr); 157 ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 158 ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 159 ierr = MatMatMultEqual(A,B,P,10,&flg);CHKERRQ(ierr); 160 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ABx != Px"); 161 ierr = MatTransposeMatMultEqual(A,B,P,10,&flg);CHKERRQ(ierr); 162 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^TBx != P^Tx"); 163 ierr = MatDestroy(&B);CHKERRQ(ierr); 164 ierr = MatDestroy(&P);CHKERRQ(ierr); 165 if (n) { 166 ierr = PetscMalloc1(n*dim,&scoords);CHKERRQ(ierr); 167 ierr = PetscRandomGetValuesReal(rdm,n*dim,scoords);CHKERRQ(ierr); 168 N = n; 169 ierr = MPIU_Allreduce(MPI_IN_PLACE,&N,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRMPI(ierr); 170 ierr = PetscCalloc1(N*dim,&gscoords);CHKERRQ(ierr); 171 ierr = MPI_Exscan(&n,&begin,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRMPI(ierr); 172 ierr = PetscArraycpy(gscoords+begin*dim,scoords,n*dim);CHKERRQ(ierr); 173 ierr = MPIU_Allreduce(MPI_IN_PLACE,gscoords,N*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD);CHKERRMPI(ierr); 174 kernel = GenEntriesRectangular; 175 ctx[0] = gcoords; 176 ctx[1] = gscoords; 177 ierr = MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,n,M,N,dim,coords,scoords,kernel,ctx,&R);CHKERRQ(ierr); 178 ierr = MatSetFromOptions(R);CHKERRQ(ierr); 179 ierr = MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 180 ierr = MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 181 ierr = MatViewFromOptions(R,NULL,"-R_view");CHKERRQ(ierr); 182 ierr = MatConvert(R,MATDENSE,MAT_INITIAL_MATRIX,&D);CHKERRQ(ierr); 183 ierr = MatViewFromOptions(D,NULL,"-D_view");CHKERRQ(ierr); 184 ierr = MatMultEqual(R,D,10,&flg);CHKERRQ(ierr); 185 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Rx != Dx"); 186 ierr = MatTranspose(D,MAT_INPLACE_MATRIX,&D);CHKERRQ(ierr); 187 ierr = MatTranspose(R,MAT_INITIAL_MATRIX,&RT);CHKERRQ(ierr); 188 ierr = MatMultEqual(RT,D,10,&flg);CHKERRQ(ierr); 189 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"R^Tx != D^Tx"); 190 ierr = MatTranspose(R,MAT_REUSE_MATRIX,&RT);CHKERRQ(ierr); 191 ierr = MatMultEqual(RT,D,10,&flg);CHKERRQ(ierr); 192 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"R^Tx != D^Tx"); 193 ierr = MatDestroy(&RT);CHKERRQ(ierr); 194 ierr = MatDestroy(&D);CHKERRQ(ierr); 195 ierr = MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,PETSC_DETERMINE,K,NULL,&B);CHKERRQ(ierr); 196 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 197 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 198 ierr = MatSetRandom(B,rdm);CHKERRQ(ierr); 199 ierr = MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&P);CHKERRQ(ierr); 200 ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 201 ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 202 ierr = MatMatMultEqual(R,B,P,10,&flg);CHKERRQ(ierr); 203 PetscCheckFalse(!flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"RBx != Px"); 204 ierr = MatDestroy(&B);CHKERRQ(ierr); 205 ierr = MatDestroy(&P);CHKERRQ(ierr); 206 ierr = MatCreateVecs(R,&right,&left);CHKERRQ(ierr); 207 ierr = VecSetRandom(right,rdm);CHKERRQ(ierr); 208 ierr = MatMult(R,right,left);CHKERRQ(ierr); 209 ierr = MatHtoolGetPermutationSource(R,&iss);CHKERRQ(ierr); 210 ierr = MatHtoolGetPermutationTarget(R,&ist);CHKERRQ(ierr); 211 ierr = VecDuplicate(left,&perm);CHKERRQ(ierr); 212 ierr = VecCopy(left,perm);CHKERRQ(ierr); 213 ierr = VecPermute(perm,ist,PETSC_FALSE);CHKERRQ(ierr); 214 ierr = VecPermute(right,iss,PETSC_FALSE);CHKERRQ(ierr); 215 ierr = MatHtoolUsePermutation(R,PETSC_FALSE);CHKERRQ(ierr); 216 ierr = MatMult(R,right,left);CHKERRQ(ierr); 217 ierr = VecAXPY(left,-1.0,perm);CHKERRQ(ierr); 218 ierr = VecNorm(left,NORM_INFINITY,&norm);CHKERRQ(ierr); 219 PetscCheckFalse(PetscAbsReal(norm) > PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||y(with permutation)-y(without permutation)|| = %g (> %g)",(double)PetscAbsReal(norm),(double)PETSC_SMALL); 220 ierr = MatHtoolUsePermutation(R,PETSC_TRUE);CHKERRQ(ierr); 221 ierr = VecDestroy(&perm);CHKERRQ(ierr); 222 ierr = VecDestroy(&left);CHKERRQ(ierr); 223 ierr = VecDestroy(&right);CHKERRQ(ierr); 224 ierr = ISDestroy(&ist);CHKERRQ(ierr); 225 ierr = ISDestroy(&iss);CHKERRQ(ierr); 226 ierr = MatDestroy(&R);CHKERRQ(ierr); 227 ierr = PetscFree(gscoords);CHKERRQ(ierr); 228 ierr = PetscFree(scoords);CHKERRQ(ierr); 229 } 230 } 231 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 232 ierr = MatDestroy(&A);CHKERRQ(ierr); 233 ierr = PetscFree(gcoords);CHKERRQ(ierr); 234 ierr = PetscFree(coords);CHKERRQ(ierr); 235 ierr = PetscFinalize(); 236 return ierr; 237 } 238 239 /*TEST 240 241 build: 242 requires: htool 243 244 test: 245 requires: htool 246 suffix: 1 247 nsize: 4 248 args: -m_local 80 -n_local 25 -mat_htool_epsilon 1.0e-11 -symmetric {{false true}shared output} 249 output_file: output/ex101.out 250 251 test: 252 requires: htool 253 suffix: 2 254 nsize: 4 255 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} 256 output_file: output/ex101.out 257 258 TEST*/ 259