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 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m_local",&m,NULL)); 53 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n_local",&n,NULL)); 54 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL)); 55 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL)); 56 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-symmetric",&sym,NULL)); 57 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mat_htool_epsilon",&epsilon,NULL)); 58 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 59 M = size*m; 60 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 61 CHKERRQ(PetscMalloc1(m*dim,&coords)); 62 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&rdm)); 63 CHKERRQ(PetscRandomGetValuesReal(rdm,m*dim,coords)); 64 CHKERRQ(PetscCalloc1(M*dim,&gcoords)); 65 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,NULL,&B)); 66 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 67 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 68 CHKERRQ(MatSetRandom(B,rdm)); 69 CHKERRQ(MatGetOwnershipRange(B,&begin,NULL)); 70 CHKERRQ(PetscArraycpy(gcoords+begin*dim,coords,m*dim)); 71 CHKERRMPI(MPIU_Allreduce(MPI_IN_PLACE,gcoords,M*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD)); 72 CHKERRQ(MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&A)); 73 CHKERRQ(MatSetOption(A,MAT_SYMMETRIC,sym)); 74 CHKERRQ(MatSetFromOptions(A)); 75 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 76 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 77 CHKERRQ(MatViewFromOptions(A,NULL,"-A_view")); 78 CHKERRQ(MatGetOwnershipIS(A,is,NULL)); 79 CHKERRQ(ISDuplicate(is[0],is+1)); 80 CHKERRQ(MatIncreaseOverlap(A,1,is,2)); 81 CHKERRQ(MatSetBlockSize(A,2)); 82 CHKERRQ(MatIncreaseOverlap(A,1,is+1,1)); 83 CHKERRQ(ISGetBlockSize(is[1],&bs)); 84 PetscCheckFalse(bs != 2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrect block size %" PetscInt_FMT " != 2",bs); 85 CHKERRQ(MatSetBlockSize(A,1)); 86 CHKERRQ(ISEqual(is[0],is[1],&flg)); 87 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unequal index sets"); 88 CHKERRQ(ISDestroy(is)); 89 CHKERRQ(ISDestroy(is+1)); 90 CHKERRQ(MatCreateVecs(A,&right,&left)); 91 CHKERRQ(VecSetRandom(right,rdm)); 92 CHKERRQ(MatMult(A,right,left)); 93 CHKERRQ(MatHtoolGetPermutationSource(A,&iss)); 94 CHKERRQ(MatHtoolGetPermutationTarget(A,&ist)); 95 CHKERRQ(VecDuplicate(left,&perm)); 96 CHKERRQ(VecCopy(left,perm)); 97 CHKERRQ(VecPermute(perm,ist,PETSC_FALSE)); 98 CHKERRQ(VecPermute(right,iss,PETSC_FALSE)); 99 CHKERRQ(MatHtoolUsePermutation(A,PETSC_FALSE)); 100 CHKERRQ(MatMult(A,right,left)); 101 CHKERRQ(VecAXPY(left,-1.0,perm)); 102 CHKERRQ(VecNorm(left,NORM_INFINITY,&norm)); 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 CHKERRQ(MatHtoolUsePermutation(A,PETSC_TRUE)); 105 CHKERRQ(VecDestroy(&perm)); 106 CHKERRQ(VecDestroy(&left)); 107 CHKERRQ(VecDestroy(&right)); 108 CHKERRQ(ISDestroy(&ist)); 109 CHKERRQ(ISDestroy(&iss)); 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 CHKERRQ(MatDestroy(&B)); 113 CHKERRQ(MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,m,M,M,dim,coords,coords,kernel,gcoords,&B)); 114 CHKERRQ(MatSetOption(B,MAT_SYMMETRIC,(PetscBool)!sym)); 115 CHKERRQ(MatSetFromOptions(B)); 116 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 117 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 118 CHKERRQ(MatViewFromOptions(B,NULL,"-B_view")); 119 CHKERRQ(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&P)); 120 CHKERRQ(MatNorm(P,NORM_FROBENIUS,&relative)); 121 CHKERRQ(MatConvert(B,MATDENSE,MAT_INITIAL_MATRIX,&R)); 122 CHKERRQ(MatAXPY(R,-1.0,P,SAME_NONZERO_PATTERN)); 123 CHKERRQ(MatNorm(R,NORM_INFINITY,&norm)); 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 CHKERRQ(MatDestroy(&B)); 126 CHKERRQ(MatDestroy(&R)); 127 CHKERRQ(MatDestroy(&P)); 128 } else { 129 CHKERRQ(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&D)); 130 CHKERRQ(MatViewFromOptions(D,NULL,"-D_view")); 131 CHKERRQ(MatMultEqual(A,D,10,&flg)); 132 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Ax != Dx"); 133 CHKERRQ(MatMultTransposeEqual(A,D,10,&flg)); 134 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^Tx != D^Tx"); 135 CHKERRQ(MatMultAddEqual(A,D,10,&flg)); 136 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"y+Ax != y+Dx"); 137 CHKERRQ(MatGetOwnershipRange(B,&begin,NULL)); 138 CHKERRQ(MatDenseGetArrayWrite(D,&ptr)); 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 CHKERRQ(MatDenseRestoreArrayWrite(D,&ptr)); 142 CHKERRQ(MatMultEqual(A,D,10,&flg)); 143 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Ax != Dx"); 144 CHKERRQ(MatTranspose(D,MAT_INPLACE_MATRIX,&D)); 145 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&AT)); 146 CHKERRQ(MatMultEqual(AT,D,10,&flg)); 147 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^Tx != D^Tx"); 148 CHKERRQ(MatTranspose(A,MAT_REUSE_MATRIX,&AT)); 149 CHKERRQ(MatMultEqual(AT,D,10,&flg)); 150 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^Tx != D^Tx"); 151 CHKERRQ(MatAXPY(D,-1.0,AT,SAME_NONZERO_PATTERN)); 152 CHKERRQ(MatNorm(D,NORM_INFINITY,&norm)); 153 PetscCheckFalse(PetscAbsReal(norm) > PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"||A-D|| = %g (> %g)",(double)norm,(double)PETSC_SMALL); 154 CHKERRQ(MatDestroy(&AT)); 155 CHKERRQ(MatDestroy(&D)); 156 CHKERRQ(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&P)); 157 CHKERRQ(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY)); 158 CHKERRQ(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY)); 159 CHKERRQ(MatMatMultEqual(A,B,P,10,&flg)); 160 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ABx != Px"); 161 CHKERRQ(MatTransposeMatMultEqual(A,B,P,10,&flg)); 162 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A^TBx != P^Tx"); 163 CHKERRQ(MatDestroy(&B)); 164 CHKERRQ(MatDestroy(&P)); 165 if (n) { 166 CHKERRQ(PetscMalloc1(n*dim,&scoords)); 167 CHKERRQ(PetscRandomGetValuesReal(rdm,n*dim,scoords)); 168 N = n; 169 CHKERRMPI(MPIU_Allreduce(MPI_IN_PLACE,&N,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD)); 170 CHKERRQ(PetscCalloc1(N*dim,&gscoords)); 171 CHKERRMPI(MPI_Exscan(&n,&begin,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD)); 172 CHKERRQ(PetscArraycpy(gscoords+begin*dim,scoords,n*dim)); 173 CHKERRMPI(MPIU_Allreduce(MPI_IN_PLACE,gscoords,N*dim,MPIU_REAL,MPI_SUM,PETSC_COMM_WORLD)); 174 kernel = GenEntriesRectangular; 175 ctx[0] = gcoords; 176 ctx[1] = gscoords; 177 CHKERRQ(MatCreateHtoolFromKernel(PETSC_COMM_WORLD,m,n,M,N,dim,coords,scoords,kernel,ctx,&R)); 178 CHKERRQ(MatSetFromOptions(R)); 179 CHKERRQ(MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY)); 180 CHKERRQ(MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY)); 181 CHKERRQ(MatViewFromOptions(R,NULL,"-R_view")); 182 CHKERRQ(MatConvert(R,MATDENSE,MAT_INITIAL_MATRIX,&D)); 183 CHKERRQ(MatViewFromOptions(D,NULL,"-D_view")); 184 CHKERRQ(MatMultEqual(R,D,10,&flg)); 185 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Rx != Dx"); 186 CHKERRQ(MatTranspose(D,MAT_INPLACE_MATRIX,&D)); 187 CHKERRQ(MatTranspose(R,MAT_INITIAL_MATRIX,&RT)); 188 CHKERRQ(MatMultEqual(RT,D,10,&flg)); 189 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"R^Tx != D^Tx"); 190 CHKERRQ(MatTranspose(R,MAT_REUSE_MATRIX,&RT)); 191 CHKERRQ(MatMultEqual(RT,D,10,&flg)); 192 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"R^Tx != D^Tx"); 193 CHKERRQ(MatDestroy(&RT)); 194 CHKERRQ(MatDestroy(&D)); 195 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,PETSC_DETERMINE,K,NULL,&B)); 196 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 197 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 198 CHKERRQ(MatSetRandom(B,rdm)); 199 CHKERRQ(MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&P)); 200 CHKERRQ(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY)); 201 CHKERRQ(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY)); 202 CHKERRQ(MatMatMultEqual(R,B,P,10,&flg)); 203 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"RBx != Px"); 204 CHKERRQ(MatDestroy(&B)); 205 CHKERRQ(MatDestroy(&P)); 206 CHKERRQ(MatCreateVecs(R,&right,&left)); 207 CHKERRQ(VecSetRandom(right,rdm)); 208 CHKERRQ(MatMult(R,right,left)); 209 CHKERRQ(MatHtoolGetPermutationSource(R,&iss)); 210 CHKERRQ(MatHtoolGetPermutationTarget(R,&ist)); 211 CHKERRQ(VecDuplicate(left,&perm)); 212 CHKERRQ(VecCopy(left,perm)); 213 CHKERRQ(VecPermute(perm,ist,PETSC_FALSE)); 214 CHKERRQ(VecPermute(right,iss,PETSC_FALSE)); 215 CHKERRQ(MatHtoolUsePermutation(R,PETSC_FALSE)); 216 CHKERRQ(MatMult(R,right,left)); 217 CHKERRQ(VecAXPY(left,-1.0,perm)); 218 CHKERRQ(VecNorm(left,NORM_INFINITY,&norm)); 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 CHKERRQ(MatHtoolUsePermutation(R,PETSC_TRUE)); 221 CHKERRQ(VecDestroy(&perm)); 222 CHKERRQ(VecDestroy(&left)); 223 CHKERRQ(VecDestroy(&right)); 224 CHKERRQ(ISDestroy(&ist)); 225 CHKERRQ(ISDestroy(&iss)); 226 CHKERRQ(MatDestroy(&R)); 227 CHKERRQ(PetscFree(gscoords)); 228 CHKERRQ(PetscFree(scoords)); 229 } 230 } 231 CHKERRQ(PetscRandomDestroy(&rdm)); 232 CHKERRQ(MatDestroy(&A)); 233 CHKERRQ(PetscFree(gcoords)); 234 CHKERRQ(PetscFree(coords)); 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