1 #include <../src/mat/impls/htool/htool.hpp> /*I "petscmat.h" I*/ 2 #include <petscblaslapack.h> 3 #include <set> 4 5 const char *const MatHtoolCompressorTypes[] = { "sympartialACA", "fullACA", "SVD" }; 6 const char *const MatHtoolClusteringTypes[] = { "PCARegular", "PCAGeometric", "BoundingBox1Regular", "BoundingBox1Geometric" }; 7 const char HtoolCitation[] = "@article{marchand2020two,\n" 8 " Author = {Marchand, Pierre and Claeys, Xavier and Jolivet, Pierre and Nataf, Fr\\'ed\\'eric and Tournier, Pierre-Henri},\n" 9 " Title = {Two-level preconditioning for $h$-version boundary element approximation of hypersingular operator with {GenEO}},\n" 10 " Year = {2020},\n" 11 " Publisher = {Elsevier},\n" 12 " Journal = {Numerische Mathematik},\n" 13 " Volume = {146},\n" 14 " Pages = {597--628},\n" 15 " Url = {https://github.com/htool-ddm/htool}\n" 16 "}\n"; 17 static PetscBool HtoolCite = PETSC_FALSE; 18 19 static PetscErrorCode MatGetDiagonal_Htool(Mat A,Vec v) 20 { 21 Mat_Htool *a = (Mat_Htool*)A->data; 22 PetscScalar *x; 23 PetscBool flg; 24 25 PetscFunctionBegin; 26 PetscCall(MatHasCongruentLayouts(A,&flg)); 27 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only congruent layouts supported"); 28 PetscCall(VecGetArrayWrite(v,&x)); 29 a->hmatrix->copy_local_diagonal(x); 30 PetscCall(VecRestoreArrayWrite(v,&x)); 31 PetscCall(VecScale(v,a->s)); 32 PetscFunctionReturn(0); 33 } 34 35 static PetscErrorCode MatGetDiagonalBlock_Htool(Mat A,Mat *b) 36 { 37 Mat_Htool *a = (Mat_Htool*)A->data; 38 Mat B; 39 PetscScalar *ptr; 40 PetscBool flg; 41 42 PetscFunctionBegin; 43 PetscCall(MatHasCongruentLayouts(A,&flg)); 44 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only congruent layouts supported"); 45 PetscCall(PetscObjectQuery((PetscObject)A,"DiagonalBlock",(PetscObject*)&B)); /* same logic as in MatGetDiagonalBlock_MPIDense() */ 46 if (!B) { 47 PetscCall(MatCreateDense(PETSC_COMM_SELF,A->rmap->n,A->rmap->n,A->rmap->n,A->rmap->n,NULL,&B)); 48 PetscCall(MatDenseGetArrayWrite(B,&ptr)); 49 a->hmatrix->copy_local_diagonal_block(ptr); 50 PetscCall(MatDenseRestoreArrayWrite(B,&ptr)); 51 PetscCall(MatPropagateSymmetryOptions(A,B)); 52 PetscCall(MatScale(B,a->s)); 53 PetscCall(PetscObjectCompose((PetscObject)A,"DiagonalBlock",(PetscObject)B)); 54 *b = B; 55 PetscCall(MatDestroy(&B)); 56 } else *b = B; 57 PetscFunctionReturn(0); 58 } 59 60 static PetscErrorCode MatMult_Htool(Mat A,Vec x,Vec y) 61 { 62 Mat_Htool *a = (Mat_Htool*)A->data; 63 const PetscScalar *in; 64 PetscScalar *out; 65 66 PetscFunctionBegin; 67 PetscCall(VecGetArrayRead(x,&in)); 68 PetscCall(VecGetArrayWrite(y,&out)); 69 a->hmatrix->mvprod_local_to_local(in,out); 70 PetscCall(VecRestoreArrayRead(x,&in)); 71 PetscCall(VecRestoreArrayWrite(y,&out)); 72 PetscCall(VecScale(y,a->s)); 73 PetscFunctionReturn(0); 74 } 75 76 /* naive implementation of MatMultAdd() needed for FEM-BEM coupling via MATNEST */ 77 static PetscErrorCode MatMultAdd_Htool(Mat A,Vec v1,Vec v2,Vec v3) 78 { 79 Mat_Htool *a = (Mat_Htool*)A->data; 80 Vec tmp; 81 const PetscScalar scale = a->s; 82 83 PetscFunctionBegin; 84 PetscCall(VecDuplicate(v2,&tmp)); 85 PetscCall(VecCopy(v2,v3)); /* no-op in MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]) since VecCopy() checks for x == y */ 86 a->s = 1.0; /* set s to 1.0 since VecAXPY() may be used to scale the MatMult() output Vec */ 87 PetscCall(MatMult_Htool(A,v1,tmp)); 88 PetscCall(VecAXPY(v3,scale,tmp)); 89 PetscCall(VecDestroy(&tmp)); 90 a->s = scale; /* set s back to its original value */ 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode MatMultTranspose_Htool(Mat A,Vec x,Vec y) 95 { 96 Mat_Htool *a = (Mat_Htool*)A->data; 97 const PetscScalar *in; 98 PetscScalar *out; 99 100 PetscFunctionBegin; 101 PetscCall(VecGetArrayRead(x,&in)); 102 PetscCall(VecGetArrayWrite(y,&out)); 103 a->hmatrix->mvprod_transp_local_to_local(in,out); 104 PetscCall(VecRestoreArrayRead(x,&in)); 105 PetscCall(VecRestoreArrayWrite(y,&out)); 106 PetscCall(VecScale(y,a->s)); 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode MatIncreaseOverlap_Htool(Mat A,PetscInt is_max,IS is[],PetscInt ov) 111 { 112 std::set<PetscInt> set; 113 const PetscInt *idx; 114 PetscInt *oidx,size,bs[2]; 115 PetscMPIInt csize; 116 117 PetscFunctionBegin; 118 PetscCall(MatGetBlockSizes(A,bs,bs+1)); 119 if (bs[0] != bs[1]) bs[0] = 1; 120 for (PetscInt i=0; i<is_max; ++i) { 121 /* basic implementation that adds indices by shifting an IS by -ov, -ov+1..., -1, 1..., ov-1, ov */ 122 /* needed to avoid subdomain matrices to replicate A since it is dense */ 123 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is[i]),&csize)); 124 PetscCheck(csize == 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported parallel IS"); 125 PetscCall(ISGetSize(is[i],&size)); 126 PetscCall(ISGetIndices(is[i],&idx)); 127 for (PetscInt j=0; j<size; ++j) { 128 set.insert(idx[j]); 129 for (PetscInt k=1; k<=ov; ++k) { /* for each layer of overlap */ 130 if (idx[j] - k >= 0) set.insert(idx[j] - k); /* do not insert negative indices */ 131 if (idx[j] + k < A->rmap->N && idx[j] + k < A->cmap->N) set.insert(idx[j] + k); /* do not insert indices greater than the dimension of A */ 132 } 133 } 134 PetscCall(ISRestoreIndices(is[i],&idx)); 135 PetscCall(ISDestroy(is+i)); 136 if (bs[0] > 1) { 137 for (std::set<PetscInt>::iterator it=set.cbegin(); it!=set.cend(); it++) { 138 std::vector<PetscInt> block(bs[0]); 139 std::iota(block.begin(),block.end(),(*it/bs[0])*bs[0]); 140 set.insert(block.cbegin(),block.cend()); 141 } 142 } 143 size = set.size(); /* size with overlap */ 144 PetscCall(PetscMalloc1(size,&oidx)); 145 for (const PetscInt j : set) *oidx++ = j; 146 oidx -= size; 147 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,size,oidx,PETSC_OWN_POINTER,is+i)); 148 } 149 PetscFunctionReturn(0); 150 } 151 152 static PetscErrorCode MatCreateSubMatrices_Htool(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[]) 153 { 154 Mat_Htool *a = (Mat_Htool*)A->data; 155 Mat D,B,BT; 156 const PetscScalar *copy; 157 PetscScalar *ptr; 158 const PetscInt *idxr,*idxc,*it; 159 PetscInt nrow,m,i; 160 PetscBool flg; 161 162 PetscFunctionBegin; 163 if (scall != MAT_REUSE_MATRIX) { 164 PetscCall(PetscCalloc1(n,submat)); 165 } 166 for (i=0; i<n; ++i) { 167 PetscCall(ISGetLocalSize(irow[i],&nrow)); 168 PetscCall(ISGetLocalSize(icol[i],&m)); 169 PetscCall(ISGetIndices(irow[i],&idxr)); 170 PetscCall(ISGetIndices(icol[i],&idxc)); 171 if (scall != MAT_REUSE_MATRIX) { 172 PetscCall(MatCreateDense(PETSC_COMM_SELF,nrow,m,nrow,m,NULL,(*submat)+i)); 173 } 174 PetscCall(MatDenseGetArrayWrite((*submat)[i],&ptr)); 175 if (irow[i] == icol[i]) { /* same row and column IS? */ 176 PetscCall(MatHasCongruentLayouts(A,&flg)); 177 if (flg) { 178 PetscCall(ISSorted(irow[i],&flg)); 179 if (flg) { /* sorted IS? */ 180 it = std::lower_bound(idxr,idxr+nrow,A->rmap->rstart); 181 if (it != idxr+nrow && *it == A->rmap->rstart) { /* rmap->rstart in IS? */ 182 if (std::distance(idxr,it) + A->rmap->n <= nrow) { /* long enough IS to store the local diagonal block? */ 183 for (PetscInt j=0; j<A->rmap->n && flg; ++j) if (PetscUnlikely(it[j] != A->rmap->rstart+j)) flg = PETSC_FALSE; 184 if (flg) { /* complete local diagonal block in IS? */ 185 /* fast extraction when the local diagonal block is part of the submatrix, e.g., for PCASM or PCHPDDM 186 * [ B C E ] 187 * A = [ B D E ] 188 * [ B F E ] 189 */ 190 m = std::distance(idxr,it); /* shift of the coefficient (0,0) of block D from above */ 191 PetscCall(MatGetDiagonalBlock_Htool(A,&D)); 192 PetscCall(MatDenseGetArrayRead(D,©)); 193 for (PetscInt k=0; k<A->rmap->n; ++k) { 194 PetscCall(PetscArraycpy(ptr+(m+k)*nrow+m,copy+k*A->rmap->n,A->rmap->n)); /* block D from above */ 195 } 196 PetscCall(MatDenseRestoreArrayRead(D,©)); 197 if (m) { 198 a->wrapper->copy_submatrix(nrow,m,idxr,idxc,ptr); /* vertical block B from above */ 199 /* entry-wise assembly may be costly, so transpose already-computed entries when possible */ 200 if (A->symmetric || A->hermitian) { 201 PetscCall(MatCreateDense(PETSC_COMM_SELF,A->rmap->n,m,A->rmap->n,m,ptr+m,&B)); 202 PetscCall(MatDenseSetLDA(B,nrow)); 203 PetscCall(MatCreateDense(PETSC_COMM_SELF,m,A->rmap->n,m,A->rmap->n,ptr+m*nrow,&BT)); 204 PetscCall(MatDenseSetLDA(BT,nrow)); 205 if (A->hermitian && PetscDefined(USE_COMPLEX)) { 206 PetscCall(MatHermitianTranspose(B,MAT_REUSE_MATRIX,&BT)); 207 } else { 208 PetscCall(MatTranspose(B,MAT_REUSE_MATRIX,&BT)); 209 } 210 PetscCall(MatDestroy(&B)); 211 PetscCall(MatDestroy(&BT)); 212 } else { 213 for (PetscInt k=0; k<A->rmap->n; ++k) { /* block C from above */ 214 a->wrapper->copy_submatrix(m,1,idxr,idxc+m+k,ptr+(m+k)*nrow); 215 } 216 } 217 } 218 if (m+A->rmap->n != nrow) { 219 a->wrapper->copy_submatrix(nrow,std::distance(it+A->rmap->n,idxr+nrow),idxr,idxc+m+A->rmap->n,ptr+(m+A->rmap->n)*nrow); /* vertical block E from above */ 220 /* entry-wise assembly may be costly, so transpose already-computed entries when possible */ 221 if (A->symmetric || A->hermitian) { 222 PetscCall(MatCreateDense(PETSC_COMM_SELF,A->rmap->n,nrow-(m+A->rmap->n),A->rmap->n,nrow-(m+A->rmap->n),ptr+(m+A->rmap->n)*nrow+m,&B)); 223 PetscCall(MatDenseSetLDA(B,nrow)); 224 PetscCall(MatCreateDense(PETSC_COMM_SELF,nrow-(m+A->rmap->n),A->rmap->n,nrow-(m+A->rmap->n),A->rmap->n,ptr+m*nrow+m+A->rmap->n,&BT)); 225 PetscCall(MatDenseSetLDA(BT,nrow)); 226 if (A->hermitian && PetscDefined(USE_COMPLEX)) { 227 PetscCall(MatHermitianTranspose(B,MAT_REUSE_MATRIX,&BT)); 228 } else { 229 PetscCall(MatTranspose(B,MAT_REUSE_MATRIX,&BT)); 230 } 231 PetscCall(MatDestroy(&B)); 232 PetscCall(MatDestroy(&BT)); 233 } else { 234 for (PetscInt k=0; k<A->rmap->n; ++k) { /* block F from above */ 235 a->wrapper->copy_submatrix(std::distance(it+A->rmap->n,idxr+nrow),1,it+A->rmap->n,idxc+m+k,ptr+(m+k)*nrow+m+A->rmap->n); 236 } 237 } 238 } 239 } /* complete local diagonal block not in IS */ 240 } else flg = PETSC_FALSE; /* IS not long enough to store the local diagonal block */ 241 } else flg = PETSC_FALSE; /* rmap->rstart not in IS */ 242 } /* unsorted IS */ 243 } 244 } else flg = PETSC_FALSE; /* different row and column IS */ 245 if (!flg) a->wrapper->copy_submatrix(nrow,m,idxr,idxc,ptr); /* reassemble everything */ 246 PetscCall(ISRestoreIndices(irow[i],&idxr)); 247 PetscCall(ISRestoreIndices(icol[i],&idxc)); 248 PetscCall(MatDenseRestoreArrayWrite((*submat)[i],&ptr)); 249 PetscCall(MatScale((*submat)[i],a->s)); 250 } 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode MatDestroy_Htool(Mat A) 255 { 256 Mat_Htool *a = (Mat_Htool*)A->data; 257 PetscContainer container; 258 MatHtoolKernelTranspose *kernelt; 259 260 PetscFunctionBegin; 261 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 262 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_seqdense_C",NULL)); 263 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_mpidense_C",NULL)); 264 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_seqdense_C",NULL)); 265 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_mpidense_C",NULL)); 266 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetHierarchicalMat_C",NULL)); 267 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolSetKernel_C",NULL)); 268 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationSource_C",NULL)); 269 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationTarget_C",NULL)); 270 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolUsePermutation_C",NULL)); 271 PetscCall(PetscObjectQuery((PetscObject)A,"KernelTranspose",(PetscObject*)&container)); 272 if (container) { /* created in MatTranspose_Htool() */ 273 PetscCall(PetscContainerGetPointer(container,(void**)&kernelt)); 274 PetscCall(MatDestroy(&kernelt->A)); 275 PetscCall(PetscFree(kernelt)); 276 PetscCall(PetscContainerDestroy(&container)); 277 PetscCall(PetscObjectCompose((PetscObject)A,"KernelTranspose",NULL)); 278 } 279 if (a->gcoords_source != a->gcoords_target) { 280 PetscCall(PetscFree(a->gcoords_source)); 281 } 282 PetscCall(PetscFree(a->gcoords_target)); 283 PetscCall(PetscFree2(a->work_source,a->work_target)); 284 delete a->wrapper; 285 delete a->hmatrix; 286 PetscCall(PetscFree(A->data)); 287 PetscFunctionReturn(0); 288 } 289 290 static PetscErrorCode MatView_Htool(Mat A,PetscViewer pv) 291 { 292 Mat_Htool *a = (Mat_Htool*)A->data; 293 PetscBool flg; 294 295 PetscFunctionBegin; 296 PetscCall(PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&flg)); 297 if (flg) { 298 PetscCall(PetscViewerASCIIPrintf(pv,"symmetry: %c\n",a->hmatrix->get_symmetry_type())); 299 if (PetscAbsScalar(a->s-1.0) > PETSC_MACHINE_EPSILON) { 300 #if defined(PETSC_USE_COMPLEX) 301 PetscCall(PetscViewerASCIIPrintf(pv,"scaling: %g+%gi\n",(double)PetscRealPart(a->s),(double)PetscImaginaryPart(a->s))); 302 #else 303 PetscCall(PetscViewerASCIIPrintf(pv,"scaling: %g\n",(double)a->s)); 304 #endif 305 } 306 PetscCall(PetscViewerASCIIPrintf(pv,"minimum cluster size: %" PetscInt_FMT "\n",a->bs[0])); 307 PetscCall(PetscViewerASCIIPrintf(pv,"maximum block size: %" PetscInt_FMT "\n",a->bs[1])); 308 PetscCall(PetscViewerASCIIPrintf(pv,"epsilon: %g\n",(double)a->epsilon)); 309 PetscCall(PetscViewerASCIIPrintf(pv,"eta: %g\n",(double)a->eta)); 310 PetscCall(PetscViewerASCIIPrintf(pv,"minimum target depth: %" PetscInt_FMT "\n",a->depth[0])); 311 PetscCall(PetscViewerASCIIPrintf(pv,"minimum source depth: %" PetscInt_FMT "\n",a->depth[1])); 312 PetscCall(PetscViewerASCIIPrintf(pv,"compressor: %s\n",MatHtoolCompressorTypes[a->compressor])); 313 PetscCall(PetscViewerASCIIPrintf(pv,"clustering: %s\n",MatHtoolClusteringTypes[a->clustering])); 314 PetscCall(PetscViewerASCIIPrintf(pv,"compression ratio: %s\n",a->hmatrix->get_infos("Compression_ratio").c_str())); 315 PetscCall(PetscViewerASCIIPrintf(pv,"space saving: %s\n",a->hmatrix->get_infos("Space_saving").c_str())); 316 PetscCall(PetscViewerASCIIPrintf(pv,"number of dense (resp. low rank) matrices: %s (resp. %s)\n",a->hmatrix->get_infos("Number_of_dmat").c_str(),a->hmatrix->get_infos("Number_of_lrmat").c_str())); 317 PetscCall(PetscViewerASCIIPrintf(pv,"(minimum, mean, maximum) dense block sizes: (%s, %s, %s)\n",a->hmatrix->get_infos("Dense_block_size_min").c_str(),a->hmatrix->get_infos("Dense_block_size_mean").c_str(),a->hmatrix->get_infos("Dense_block_size_max").c_str())); 318 PetscCall(PetscViewerASCIIPrintf(pv,"(minimum, mean, maximum) low rank block sizes: (%s, %s, %s)\n",a->hmatrix->get_infos("Low_rank_block_size_min").c_str(),a->hmatrix->get_infos("Low_rank_block_size_mean").c_str(),a->hmatrix->get_infos("Low_rank_block_size_max").c_str())); 319 PetscCall(PetscViewerASCIIPrintf(pv,"(minimum, mean, maximum) ranks: (%s, %s, %s)\n",a->hmatrix->get_infos("Rank_min").c_str(),a->hmatrix->get_infos("Rank_mean").c_str(),a->hmatrix->get_infos("Rank_max").c_str())); 320 } 321 PetscFunctionReturn(0); 322 } 323 324 static PetscErrorCode MatScale_Htool(Mat A,PetscScalar s) 325 { 326 Mat_Htool *a = (Mat_Htool*)A->data; 327 328 PetscFunctionBegin; 329 a->s *= s; 330 PetscFunctionReturn(0); 331 } 332 333 /* naive implementation of MatGetRow() needed for MatConvert_Nest_AIJ() */ 334 static PetscErrorCode MatGetRow_Htool(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 335 { 336 Mat_Htool *a = (Mat_Htool*)A->data; 337 PetscInt *idxc; 338 PetscBLASInt one = 1,bn; 339 340 PetscFunctionBegin; 341 if (nz) *nz = A->cmap->N; 342 if (idx || v) { /* even if !idx, need to set idxc for htool::copy_submatrix() */ 343 PetscCall(PetscMalloc1(A->cmap->N,&idxc)); 344 for (PetscInt i=0; i<A->cmap->N; ++i) idxc[i] = i; 345 } 346 if (idx) *idx = idxc; 347 if (v) { 348 PetscCall(PetscMalloc1(A->cmap->N,v)); 349 if (a->wrapper) a->wrapper->copy_submatrix(1,A->cmap->N,&row,idxc,*v); 350 else reinterpret_cast<htool::VirtualGenerator<PetscScalar>*>(a->kernelctx)->copy_submatrix(1,A->cmap->N,&row,idxc,*v); 351 PetscCall(PetscBLASIntCast(A->cmap->N,&bn)); 352 PetscStackCallBLAS("BLASscal",BLASscal_(&bn,&a->s,*v,&one)); 353 } 354 if (!idx) { 355 PetscCall(PetscFree(idxc)); 356 } 357 PetscFunctionReturn(0); 358 } 359 360 static PetscErrorCode MatRestoreRow_Htool(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 361 { 362 PetscFunctionBegin; 363 if (nz) *nz = 0; 364 if (idx) { 365 PetscCall(PetscFree(*idx)); 366 } 367 if (v) { 368 PetscCall(PetscFree(*v)); 369 } 370 PetscFunctionReturn(0); 371 } 372 373 static PetscErrorCode MatSetFromOptions_Htool(PetscOptionItems *PetscOptionsObject,Mat A) 374 { 375 Mat_Htool *a = (Mat_Htool*)A->data; 376 PetscInt n; 377 PetscBool flg; 378 379 PetscFunctionBegin; 380 PetscOptionsHeadBegin(PetscOptionsObject,"Htool options"); 381 PetscCall(PetscOptionsInt("-mat_htool_min_cluster_size","Minimal leaf size in cluster tree",NULL,a->bs[0],a->bs,NULL)); 382 PetscCall(PetscOptionsInt("-mat_htool_max_block_size","Maximal number of coefficients in a dense block",NULL,a->bs[1],a->bs + 1,NULL)); 383 PetscCall(PetscOptionsReal("-mat_htool_epsilon","Relative error in Frobenius norm when approximating a block",NULL,a->epsilon,&a->epsilon,NULL)); 384 PetscCall(PetscOptionsReal("-mat_htool_eta","Admissibility condition tolerance",NULL,a->eta,&a->eta,NULL)); 385 PetscCall(PetscOptionsInt("-mat_htool_min_target_depth","Minimal cluster tree depth associated with the rows",NULL,a->depth[0],a->depth,NULL)); 386 PetscCall(PetscOptionsInt("-mat_htool_min_source_depth","Minimal cluster tree depth associated with the columns",NULL,a->depth[1],a->depth + 1,NULL)); 387 n = 0; 388 PetscCall(PetscOptionsEList("-mat_htool_compressor","Type of compression","MatHtoolCompressorType",MatHtoolCompressorTypes,PETSC_STATIC_ARRAY_LENGTH(MatHtoolCompressorTypes),MatHtoolCompressorTypes[MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA],&n,&flg)); 389 if (flg) a->compressor = MatHtoolCompressorType(n); 390 n = 0; 391 PetscCall(PetscOptionsEList("-mat_htool_clustering","Type of clustering","MatHtoolClusteringType",MatHtoolClusteringTypes,PETSC_STATIC_ARRAY_LENGTH(MatHtoolClusteringTypes),MatHtoolClusteringTypes[MAT_HTOOL_CLUSTERING_PCA_REGULAR],&n,&flg)); 392 if (flg) a->clustering = MatHtoolClusteringType(n); 393 PetscOptionsHeadEnd(); 394 PetscFunctionReturn(0); 395 } 396 397 static PetscErrorCode MatAssemblyEnd_Htool(Mat A,MatAssemblyType type) 398 { 399 Mat_Htool *a = (Mat_Htool*)A->data; 400 const PetscInt *ranges; 401 PetscInt *offset; 402 PetscMPIInt size; 403 char S = PetscDefined(USE_COMPLEX) && A->hermitian ? 'H' : (A->symmetric ? 'S' : 'N'),uplo = S == 'N' ? 'N' : 'U'; 404 htool::VirtualGenerator<PetscScalar> *generator = nullptr; 405 std::shared_ptr<htool::VirtualCluster> t,s = nullptr; 406 std::shared_ptr<htool::VirtualLowRankGenerator<PetscScalar>> compressor = nullptr; 407 408 PetscFunctionBegin; 409 PetscCall(PetscCitationsRegister(HtoolCitation,&HtoolCite)); 410 delete a->wrapper; 411 delete a->hmatrix; 412 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A),&size)); 413 PetscCall(PetscMalloc1(2*size,&offset)); 414 PetscCall(MatGetOwnershipRanges(A,&ranges)); 415 for (PetscInt i=0; i<size; ++i) { 416 offset[2*i] = ranges[i]; 417 offset[2*i+1] = ranges[i+1] - ranges[i]; 418 } 419 switch (a->clustering) { 420 case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC: 421 t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim); 422 break; 423 case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC: 424 t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim); 425 break; 426 case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR: 427 t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim); 428 break; 429 default: 430 t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim); 431 } 432 t->set_minclustersize(a->bs[0]); 433 t->build(A->rmap->N,a->gcoords_target,offset); 434 if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N,A->cmap->N,a->dim,a->kernel,a->kernelctx); 435 else { 436 a->wrapper = NULL; 437 generator = reinterpret_cast<htool::VirtualGenerator<PetscScalar>*>(a->kernelctx); 438 } 439 if (a->gcoords_target != a->gcoords_source) { 440 PetscCall(MatGetOwnershipRangesColumn(A,&ranges)); 441 for (PetscInt i=0; i<size; ++i) { 442 offset[2*i] = ranges[i]; 443 offset[2*i+1] = ranges[i+1] - ranges[i]; 444 } 445 switch (a->clustering) { 446 case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC: 447 s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim); 448 break; 449 case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC: 450 s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim); 451 break; 452 case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR: 453 s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim); 454 break; 455 default: 456 s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim); 457 } 458 s->set_minclustersize(a->bs[0]); 459 s->build(A->cmap->N,a->gcoords_source,offset); 460 S = uplo = 'N'; 461 } 462 PetscCall(PetscFree(offset)); 463 switch (a->compressor) { 464 case MAT_HTOOL_COMPRESSOR_FULL_ACA: 465 compressor = std::make_shared<htool::fullACA<PetscScalar>>(); 466 break; 467 case MAT_HTOOL_COMPRESSOR_SVD: 468 compressor = std::make_shared<htool::SVD<PetscScalar>>(); 469 break; 470 default: 471 compressor = std::make_shared<htool::sympartialACA<PetscScalar>>(); 472 } 473 a->hmatrix = dynamic_cast<htool::VirtualHMatrix<PetscScalar>*>(new htool::HMatrix<PetscScalar>(t,s ? s : t,a->epsilon,a->eta,S,uplo)); 474 a->hmatrix->set_compression(compressor); 475 a->hmatrix->set_maxblocksize(a->bs[1]); 476 a->hmatrix->set_mintargetdepth(a->depth[0]); 477 a->hmatrix->set_minsourcedepth(a->depth[1]); 478 if (s) a->hmatrix->build(a->wrapper ? *a->wrapper : *generator,a->gcoords_target,a->gcoords_source); 479 else a->hmatrix->build(a->wrapper ? *a->wrapper : *generator,a->gcoords_target); 480 PetscFunctionReturn(0); 481 } 482 483 static PetscErrorCode MatProductNumeric_Htool(Mat C) 484 { 485 Mat_Product *product = C->product; 486 Mat_Htool *a = (Mat_Htool*)product->A->data; 487 const PetscScalar *in; 488 PetscScalar *out; 489 PetscInt N,lda; 490 491 PetscFunctionBegin; 492 MatCheckProduct(C,1); 493 PetscCall(MatGetSize(C,NULL,&N)); 494 PetscCall(MatDenseGetLDA(C,&lda)); 495 PetscCheck(lda == C->rmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")",lda,C->rmap->n); 496 PetscCall(MatDenseGetArrayRead(product->B,&in)); 497 PetscCall(MatDenseGetArrayWrite(C,&out)); 498 switch (product->type) { 499 case MATPRODUCT_AB: 500 a->hmatrix->mvprod_local_to_local(in,out,N); 501 break; 502 case MATPRODUCT_AtB: 503 a->hmatrix->mvprod_transp_local_to_local(in,out,N); 504 break; 505 default: 506 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductType %s is not supported",MatProductTypes[product->type]); 507 } 508 PetscCall(MatDenseRestoreArrayWrite(C,&out)); 509 PetscCall(MatDenseRestoreArrayRead(product->B,&in)); 510 PetscCall(MatScale(C,a->s)); 511 PetscFunctionReturn(0); 512 } 513 514 static PetscErrorCode MatProductSymbolic_Htool(Mat C) 515 { 516 Mat_Product *product = C->product; 517 Mat A,B; 518 PetscBool flg; 519 520 PetscFunctionBegin; 521 MatCheckProduct(C,1); 522 A = product->A; 523 B = product->B; 524 PetscCall(PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,"")); 525 PetscCheck(flg,PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"MatProduct_AB not supported for %s",((PetscObject)product->B)->type_name); 526 switch (product->type) { 527 case MATPRODUCT_AB: 528 if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) { 529 PetscCall(MatSetSizes(C,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N)); 530 } 531 break; 532 case MATPRODUCT_AtB: 533 if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) { 534 PetscCall(MatSetSizes(C,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N)); 535 } 536 break; 537 default: 538 SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[product->type]); 539 } 540 PetscCall(MatSetType(C,MATDENSE)); 541 PetscCall(MatSetUp(C)); 542 PetscCall(MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE)); 543 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 544 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 545 C->ops->productsymbolic = NULL; 546 C->ops->productnumeric = MatProductNumeric_Htool; 547 PetscFunctionReturn(0); 548 } 549 550 static PetscErrorCode MatProductSetFromOptions_Htool(Mat C) 551 { 552 PetscFunctionBegin; 553 MatCheckProduct(C,1); 554 if (C->product->type == MATPRODUCT_AB || C->product->type == MATPRODUCT_AtB) C->ops->productsymbolic = MatProductSymbolic_Htool; 555 PetscFunctionReturn(0); 556 } 557 558 static PetscErrorCode MatHtoolGetHierarchicalMat_Htool(Mat A,const htool::VirtualHMatrix<PetscScalar> **hmatrix) 559 { 560 Mat_Htool *a = (Mat_Htool*)A->data; 561 562 PetscFunctionBegin; 563 *hmatrix = a->hmatrix; 564 PetscFunctionReturn(0); 565 } 566 567 /*@C 568 MatHtoolGetHierarchicalMat - Retrieves the opaque pointer to a Htool virtual matrix stored in a MATHTOOL. 569 570 Input Parameter: 571 . A - hierarchical matrix 572 573 Output Parameter: 574 . hmatrix - opaque pointer to a Htool virtual matrix 575 576 Level: advanced 577 578 .seealso: `MATHTOOL` 579 @*/ 580 PETSC_EXTERN PetscErrorCode MatHtoolGetHierarchicalMat(Mat A,const htool::VirtualHMatrix<PetscScalar> **hmatrix) 581 { 582 PetscFunctionBegin; 583 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 584 PetscValidPointer(hmatrix,2); 585 PetscTryMethod(A,"MatHtoolGetHierarchicalMat_C",(Mat,const htool::VirtualHMatrix<PetscScalar>**),(A,hmatrix)); 586 PetscFunctionReturn(0); 587 } 588 589 static PetscErrorCode MatHtoolSetKernel_Htool(Mat A,MatHtoolKernel kernel,void *kernelctx) 590 { 591 Mat_Htool *a = (Mat_Htool*)A->data; 592 593 PetscFunctionBegin; 594 a->kernel = kernel; 595 a->kernelctx = kernelctx; 596 delete a->wrapper; 597 if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N,A->cmap->N,a->dim,a->kernel,a->kernelctx); 598 PetscFunctionReturn(0); 599 } 600 601 /*@C 602 MatHtoolSetKernel - Sets the kernel and context used for the assembly of a MATHTOOL. 603 604 Input Parameters: 605 + A - hierarchical matrix 606 . kernel - computational kernel (or NULL) 607 - kernelctx - kernel context (if kernel is NULL, the pointer must be of type htool::VirtualGenerator<PetscScalar>*) 608 609 Level: advanced 610 611 .seealso: `MATHTOOL`, `MatCreateHtoolFromKernel()` 612 @*/ 613 PETSC_EXTERN PetscErrorCode MatHtoolSetKernel(Mat A,MatHtoolKernel kernel,void *kernelctx) 614 { 615 PetscFunctionBegin; 616 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 617 if (!kernelctx) PetscValidFunction(kernel,2); 618 if (!kernel) PetscValidPointer(kernelctx,3); 619 PetscTryMethod(A,"MatHtoolSetKernel_C",(Mat,MatHtoolKernel,void*),(A,kernel,kernelctx)); 620 PetscFunctionReturn(0); 621 } 622 623 static PetscErrorCode MatHtoolGetPermutationSource_Htool(Mat A,IS* is) 624 { 625 Mat_Htool *a = (Mat_Htool*)A->data; 626 std::vector<PetscInt> source; 627 628 PetscFunctionBegin; 629 source = a->hmatrix->get_source_cluster()->get_local_perm(); 630 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),source.size(),source.data(),PETSC_COPY_VALUES,is)); 631 PetscCall(ISSetPermutation(*is)); 632 PetscFunctionReturn(0); 633 } 634 635 /*@C 636 MatHtoolGetPermutationSource - Gets the permutation associated to the source cluster. 637 638 Input Parameter: 639 . A - hierarchical matrix 640 641 Output Parameter: 642 . is - permutation 643 644 Level: advanced 645 646 .seealso: `MATHTOOL`, `MatHtoolGetPermutationTarget()`, `MatHtoolUsePermutation()` 647 @*/ 648 PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationSource(Mat A,IS* is) 649 { 650 PetscFunctionBegin; 651 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 652 if (!is) PetscValidPointer(is,2); 653 PetscTryMethod(A,"MatHtoolGetPermutationSource_C",(Mat,IS*),(A,is)); 654 PetscFunctionReturn(0); 655 } 656 657 static PetscErrorCode MatHtoolGetPermutationTarget_Htool(Mat A,IS* is) 658 { 659 Mat_Htool *a = (Mat_Htool*)A->data; 660 std::vector<PetscInt> target; 661 662 PetscFunctionBegin; 663 target = a->hmatrix->get_target_cluster()->get_local_perm(); 664 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),target.size(),target.data(),PETSC_COPY_VALUES,is)); 665 PetscCall(ISSetPermutation(*is)); 666 PetscFunctionReturn(0); 667 } 668 669 /*@C 670 MatHtoolGetPermutationTarget - Gets the permutation associated to the target cluster. 671 672 Input Parameter: 673 . A - hierarchical matrix 674 675 Output Parameter: 676 . is - permutation 677 678 Level: advanced 679 680 .seealso: `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolUsePermutation()` 681 @*/ 682 PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationTarget(Mat A,IS* is) 683 { 684 PetscFunctionBegin; 685 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 686 if (!is) PetscValidPointer(is,2); 687 PetscTryMethod(A,"MatHtoolGetPermutationTarget_C",(Mat,IS*),(A,is)); 688 PetscFunctionReturn(0); 689 } 690 691 static PetscErrorCode MatHtoolUsePermutation_Htool(Mat A,PetscBool use) 692 { 693 Mat_Htool *a = (Mat_Htool*)A->data; 694 695 PetscFunctionBegin; 696 a->hmatrix->set_use_permutation(use); 697 PetscFunctionReturn(0); 698 } 699 700 /*@C 701 MatHtoolUsePermutation - Sets whether MATHTOOL should permute input (resp. output) vectors following its internal source (resp. target) permutation. 702 703 Input Parameters: 704 + A - hierarchical matrix 705 - use - Boolean value 706 707 Level: advanced 708 709 .seealso: `MATHTOOL`, `MatHtoolGetPermutationSource()`, `MatHtoolGetPermutationTarget()` 710 @*/ 711 PETSC_EXTERN PetscErrorCode MatHtoolUsePermutation(Mat A,PetscBool use) 712 { 713 PetscFunctionBegin; 714 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 715 PetscValidLogicalCollectiveBool(A,use,2); 716 PetscTryMethod(A,"MatHtoolUsePermutation_C",(Mat,PetscBool),(A,use)); 717 PetscFunctionReturn(0); 718 } 719 720 static PetscErrorCode MatConvert_Htool_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *B) 721 { 722 Mat C; 723 Mat_Htool *a = (Mat_Htool*)A->data; 724 PetscInt lda; 725 PetscScalar *array; 726 727 PetscFunctionBegin; 728 if (reuse == MAT_REUSE_MATRIX) { 729 C = *B; 730 PetscCheck(C->rmap->n == A->rmap->n && C->cmap->N == A->cmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible dimensions"); 731 PetscCall(MatDenseGetLDA(C,&lda)); 732 PetscCheck(lda == C->rmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported leading dimension (%" PetscInt_FMT " != %" PetscInt_FMT ")",lda,C->rmap->n); 733 } else { 734 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 735 PetscCall(MatSetSizes(C,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 736 PetscCall(MatSetType(C,MATDENSE)); 737 PetscCall(MatSetUp(C)); 738 PetscCall(MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE)); 739 } 740 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 741 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 742 PetscCall(MatDenseGetArrayWrite(C,&array)); 743 a->hmatrix->copy_local_dense_perm(array); 744 PetscCall(MatDenseRestoreArrayWrite(C,&array)); 745 PetscCall(MatScale(C,a->s)); 746 if (reuse == MAT_INPLACE_MATRIX) { 747 PetscCall(MatHeaderReplace(A,&C)); 748 } else *B = C; 749 PetscFunctionReturn(0); 750 } 751 752 static PetscErrorCode GenEntriesTranspose(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt *rows,const PetscInt *cols,PetscScalar *ptr,void *ctx) 753 { 754 MatHtoolKernelTranspose *generator = (MatHtoolKernelTranspose*)ctx; 755 PetscScalar *tmp; 756 757 PetscFunctionBegin; 758 generator->kernel(sdim,N,M,cols,rows,ptr,generator->kernelctx); 759 PetscCall(PetscMalloc1(M*N,&tmp)); 760 PetscCall(PetscArraycpy(tmp,ptr,M*N)); 761 for (PetscInt i=0; i<M; ++i) { 762 for (PetscInt j=0; j<N; ++j) ptr[i+j*M] = tmp[j+i*N]; 763 } 764 PetscCall(PetscFree(tmp)); 765 PetscFunctionReturn(0); 766 } 767 768 /* naive implementation which keeps a reference to the original Mat */ 769 static PetscErrorCode MatTranspose_Htool(Mat A,MatReuse reuse,Mat *B) 770 { 771 Mat C; 772 Mat_Htool *a = (Mat_Htool*)A->data,*c; 773 PetscInt M = A->rmap->N,N = A->cmap->N,m = A->rmap->n,n = A->cmap->n; 774 PetscContainer container; 775 MatHtoolKernelTranspose *kernelt; 776 777 PetscFunctionBegin; 778 PetscCheck(reuse != MAT_INPLACE_MATRIX,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatTranspose() with MAT_INPLACE_MATRIX not supported"); 779 if (reuse == MAT_INITIAL_MATRIX) { 780 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 781 PetscCall(MatSetSizes(C,n,m,N,M)); 782 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 783 PetscCall(MatSetUp(C)); 784 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)C),&container)); 785 PetscCall(PetscNew(&kernelt)); 786 PetscCall(PetscContainerSetPointer(container,kernelt)); 787 PetscCall(PetscObjectCompose((PetscObject)C,"KernelTranspose",(PetscObject)container)); 788 } else { 789 C = *B; 790 PetscCall(PetscObjectQuery((PetscObject)C,"KernelTranspose",(PetscObject*)&container)); 791 PetscCheck(container,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatTranspose() with MAT_INITIAL_MATRIX first"); 792 PetscCall(PetscContainerGetPointer(container,(void**)&kernelt)); 793 } 794 c = (Mat_Htool*)C->data; 795 c->dim = a->dim; 796 c->s = a->s; 797 c->kernel = GenEntriesTranspose; 798 if (kernelt->A != A) { 799 PetscCall(MatDestroy(&kernelt->A)); 800 kernelt->A = A; 801 PetscCall(PetscObjectReference((PetscObject)A)); 802 } 803 kernelt->kernel = a->kernel; 804 kernelt->kernelctx = a->kernelctx; 805 c->kernelctx = kernelt; 806 if (reuse == MAT_INITIAL_MATRIX) { 807 PetscCall(PetscMalloc1(N*c->dim,&c->gcoords_target)); 808 PetscCall(PetscArraycpy(c->gcoords_target,a->gcoords_source,N*c->dim)); 809 if (a->gcoords_target != a->gcoords_source) { 810 PetscCall(PetscMalloc1(M*c->dim,&c->gcoords_source)); 811 PetscCall(PetscArraycpy(c->gcoords_source,a->gcoords_target,M*c->dim)); 812 } else c->gcoords_source = c->gcoords_target; 813 PetscCall(PetscCalloc2(M,&c->work_source,N,&c->work_target)); 814 } 815 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 816 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 817 if (reuse == MAT_INITIAL_MATRIX) *B = C; 818 PetscFunctionReturn(0); 819 } 820 821 /*@C 822 MatCreateHtoolFromKernel - Creates a MATHTOOL from a user-supplied kernel. 823 824 Input Parameters: 825 + comm - MPI communicator 826 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 827 . n - number of local columns (or PETSC_DECIDE to have calculated if N is given) 828 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 829 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 830 . spacedim - dimension of the space coordinates 831 . coords_target - coordinates of the target 832 . coords_source - coordinates of the source 833 . kernel - computational kernel (or NULL) 834 - kernelctx - kernel context (if kernel is NULL, the pointer must be of type htool::VirtualGenerator<PetscScalar>*) 835 836 Output Parameter: 837 . B - matrix 838 839 Options Database Keys: 840 + -mat_htool_min_cluster_size <PetscInt> - minimal leaf size in cluster tree 841 . -mat_htool_max_block_size <PetscInt> - maximal number of coefficients in a dense block 842 . -mat_htool_epsilon <PetscReal> - relative error in Frobenius norm when approximating a block 843 . -mat_htool_eta <PetscReal> - admissibility condition tolerance 844 . -mat_htool_min_target_depth <PetscInt> - minimal cluster tree depth associated with the rows 845 . -mat_htool_min_source_depth <PetscInt> - minimal cluster tree depth associated with the columns 846 . -mat_htool_compressor <sympartialACA, fullACA, SVD> - type of compression 847 - -mat_htool_clustering <PCARegular, PCAGeometric, BounbingBox1Regular, BoundingBox1Geometric> - type of clustering 848 849 Level: intermediate 850 851 .seealso: `MatCreate()`, `MATHTOOL`, `PCSetCoordinates()`, `MatHtoolSetKernel()`, `MatHtoolCompressorType`, `MATH2OPUS`, `MatCreateH2OpusFromKernel()` 852 @*/ 853 PetscErrorCode MatCreateHtoolFromKernel(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt spacedim,const PetscReal coords_target[],const PetscReal coords_source[],MatHtoolKernel kernel,void *kernelctx,Mat *B) 854 { 855 Mat A; 856 Mat_Htool *a; 857 858 PetscFunctionBegin; 859 PetscCall(MatCreate(comm,&A)); 860 PetscValidLogicalCollectiveInt(A,spacedim,6); 861 PetscValidRealPointer(coords_target,7); 862 PetscValidRealPointer(coords_source,8); 863 if (!kernelctx) PetscValidFunction(kernel,9); 864 if (!kernel) PetscValidPointer(kernelctx,10); 865 PetscCall(MatSetSizes(A,m,n,M,N)); 866 PetscCall(MatSetType(A,MATHTOOL)); 867 PetscCall(MatSetUp(A)); 868 a = (Mat_Htool*)A->data; 869 a->dim = spacedim; 870 a->s = 1.0; 871 a->kernel = kernel; 872 a->kernelctx = kernelctx; 873 PetscCall(PetscCalloc1(A->rmap->N*spacedim,&a->gcoords_target)); 874 PetscCall(PetscArraycpy(a->gcoords_target+A->rmap->rstart*spacedim,coords_target,A->rmap->n*spacedim)); 875 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,a->gcoords_target,A->rmap->N*spacedim,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)A))); /* global target coordinates */ 876 if (coords_target != coords_source) { 877 PetscCall(PetscCalloc1(A->cmap->N*spacedim,&a->gcoords_source)); 878 PetscCall(PetscArraycpy(a->gcoords_source+A->cmap->rstart*spacedim,coords_source,A->cmap->n*spacedim)); 879 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,a->gcoords_source,A->cmap->N*spacedim,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)A))); /* global source coordinates */ 880 } else a->gcoords_source = a->gcoords_target; 881 PetscCall(PetscCalloc2(A->cmap->N,&a->work_source,A->rmap->N,&a->work_target)); 882 *B = A; 883 PetscFunctionReturn(0); 884 } 885 886 /*MC 887 MATHTOOL = "htool" - A matrix type for hierarchical matrices using the Htool package. 888 889 Use ./configure --download-htool to install PETSc to use Htool. 890 891 Options Database Keys: 892 . -mat_type htool - matrix type to "htool" during a call to MatSetFromOptions() 893 894 Level: beginner 895 896 .seealso: `MATH2OPUS`, `MATDENSE`, `MatCreateHtoolFromKernel()`, `MatHtoolSetKernel()` 897 M*/ 898 PETSC_EXTERN PetscErrorCode MatCreate_Htool(Mat A) 899 { 900 Mat_Htool *a; 901 902 PetscFunctionBegin; 903 PetscCall(PetscNewLog(A,&a)); 904 A->data = (void*)a; 905 PetscCall(PetscObjectChangeTypeName((PetscObject)A,MATHTOOL)); 906 PetscCall(PetscMemzero(A->ops,sizeof(struct _MatOps))); 907 A->ops->getdiagonal = MatGetDiagonal_Htool; 908 A->ops->getdiagonalblock = MatGetDiagonalBlock_Htool; 909 A->ops->mult = MatMult_Htool; 910 A->ops->multadd = MatMultAdd_Htool; 911 A->ops->multtranspose = MatMultTranspose_Htool; 912 if (!PetscDefined(USE_COMPLEX)) A->ops->multhermitiantranspose = MatMultTranspose_Htool; 913 A->ops->increaseoverlap = MatIncreaseOverlap_Htool; 914 A->ops->createsubmatrices = MatCreateSubMatrices_Htool; 915 A->ops->transpose = MatTranspose_Htool; 916 A->ops->destroy = MatDestroy_Htool; 917 A->ops->view = MatView_Htool; 918 A->ops->setfromoptions = MatSetFromOptions_Htool; 919 A->ops->scale = MatScale_Htool; 920 A->ops->getrow = MatGetRow_Htool; 921 A->ops->restorerow = MatRestoreRow_Htool; 922 A->ops->assemblyend = MatAssemblyEnd_Htool; 923 a->dim = 0; 924 a->gcoords_target = NULL; 925 a->gcoords_source = NULL; 926 a->s = 1.0; 927 a->bs[0] = 10; 928 a->bs[1] = 1000000; 929 a->epsilon = PetscSqrtReal(PETSC_SMALL); 930 a->eta = 10.0; 931 a->depth[0] = 0; 932 a->depth[1] = 0; 933 a->compressor = MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA; 934 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_seqdense_C",MatProductSetFromOptions_Htool)); 935 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_mpidense_C",MatProductSetFromOptions_Htool)); 936 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_seqdense_C",MatConvert_Htool_Dense)); 937 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_mpidense_C",MatConvert_Htool_Dense)); 938 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetHierarchicalMat_C",MatHtoolGetHierarchicalMat_Htool)); 939 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolSetKernel_C",MatHtoolSetKernel_Htool)); 940 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationSource_C",MatHtoolGetPermutationSource_Htool)); 941 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationTarget_C",MatHtoolGetPermutationTarget_Htool)); 942 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatHtoolUsePermutation_C",MatHtoolUsePermutation_Htool)); 943 PetscFunctionReturn(0); 944 } 945