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