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