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