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