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: %s\n",a->hmatrix->get_infos("Compression").c_str());CHKERRQ(ierr); 303 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); 304 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); 305 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); 306 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); 307 } 308 PetscFunctionReturn(0); 309 } 310 311 static PetscErrorCode MatScale_Htool(Mat A,PetscScalar s) 312 { 313 Mat_Htool *a = (Mat_Htool*)A->data; 314 315 PetscFunctionBegin; 316 a->s *= s; 317 PetscFunctionReturn(0); 318 } 319 320 /* naive implementation of MatGetRow() needed for MatConvert_Nest_AIJ() */ 321 static PetscErrorCode MatGetRow_Htool(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 322 { 323 Mat_Htool *a = (Mat_Htool*)A->data; 324 PetscInt *idxc; 325 PetscBLASInt one = 1,bn; 326 PetscErrorCode ierr; 327 328 PetscFunctionBegin; 329 if (nz) *nz = A->cmap->N; 330 if (idx || v) { /* even if !idx, need to set idxc for htool::copy_submatrix() */ 331 ierr = PetscMalloc1(A->cmap->N,&idxc);CHKERRQ(ierr); 332 for (PetscInt i=0; i<A->cmap->N; ++i) idxc[i] = i; 333 } 334 if (idx) *idx = idxc; 335 if (v) { 336 ierr = PetscMalloc1(A->cmap->N,v);CHKERRQ(ierr); 337 if (a->wrapper) a->wrapper->copy_submatrix(1,A->cmap->N,&row,idxc,*v); 338 else reinterpret_cast<htool::IMatrix<PetscScalar>*>(a->kernelctx)->copy_submatrix(1,A->cmap->N,&row,idxc,*v); 339 ierr = PetscBLASIntCast(A->cmap->N,&bn);CHKERRQ(ierr); 340 PetscStackCallBLAS("BLASscal",BLASscal_(&bn,&a->s,*v,&one)); 341 } 342 if (!idx) { 343 ierr = PetscFree(idxc);CHKERRQ(ierr); 344 } 345 PetscFunctionReturn(0); 346 } 347 348 static PetscErrorCode MatRestoreRow_Htool(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 if (nz) *nz = 0; 354 if (idx) { 355 ierr = PetscFree(*idx);CHKERRQ(ierr); 356 } 357 if (v) { 358 ierr = PetscFree(*v);CHKERRQ(ierr); 359 } 360 PetscFunctionReturn(0); 361 } 362 363 static PetscErrorCode MatSetFromOptions_Htool(PetscOptionItems *PetscOptionsObject,Mat A) 364 { 365 Mat_Htool *a = (Mat_Htool*)A->data; 366 PetscInt n; 367 PetscBool flg; 368 PetscErrorCode ierr; 369 370 PetscFunctionBegin; 371 ierr = PetscOptionsHead(PetscOptionsObject,"Htool options");CHKERRQ(ierr); 372 ierr = PetscOptionsInt("-mat_htool_min_cluster_size","Minimal leaf size in cluster tree",NULL,a->bs[0],a->bs,NULL);CHKERRQ(ierr); 373 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); 374 ierr = PetscOptionsReal("-mat_htool_epsilon","Relative error in Frobenius norm when approximating a block",NULL,a->epsilon,&a->epsilon,NULL);CHKERRQ(ierr); 375 ierr = PetscOptionsReal("-mat_htool_eta","Admissibility condition tolerance",NULL,a->eta,&a->eta,NULL);CHKERRQ(ierr); 376 ierr = PetscOptionsInt("-mat_htool_min_target_depth","Minimal cluster tree depth associated with the rows",NULL,a->depth[0],a->depth,NULL);CHKERRQ(ierr); 377 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); 378 n = 0; 379 ierr = PetscOptionsEList("-mat_htool_compressor","Type of compression","MatHtoolCompressorType",MatHtoolCompressorTypes,ALEN(MatHtoolCompressorTypes),MatHtoolCompressorTypes[MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA],&n,&flg);CHKERRQ(ierr); 380 if (flg) a->compressor = MatHtoolCompressorType(n); 381 n = 0; 382 ierr = PetscOptionsEList("-mat_htool_clustering","Type of clustering","MatHtoolClusteringType",MatHtoolClusteringTypes,ALEN(MatHtoolClusteringTypes),MatHtoolClusteringTypes[MAT_HTOOL_CLUSTERING_PCA_REGULAR],&n,&flg);CHKERRQ(ierr); 383 if (flg) a->clustering = MatHtoolClusteringType(n); 384 ierr = PetscOptionsTail();CHKERRQ(ierr); 385 PetscFunctionReturn(0); 386 } 387 388 static PetscErrorCode MatAssemblyEnd_Htool(Mat A,MatAssemblyType type) 389 { 390 Mat_Htool *a = (Mat_Htool*)A->data; 391 const PetscInt *ranges; 392 PetscInt *offset; 393 PetscMPIInt size; 394 char S = PetscDefined(USE_COMPLEX) && A->hermitian ? 'H' : (A->symmetric ? 'S' : 'N'),uplo = S == 'N' ? 'N' : 'U'; 395 htool::IMatrix<PetscScalar> *generator = nullptr; 396 std::shared_ptr<htool::VirtualCluster> t,s = nullptr; 397 PetscErrorCode ierr; 398 399 PetscFunctionBegin; 400 ierr = PetscCitationsRegister(HtoolCitation,&HtoolCite);CHKERRQ(ierr); 401 delete a->wrapper; 402 delete a->hmatrix; 403 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRMPI(ierr); 404 ierr = PetscMalloc1(2*size,&offset);CHKERRQ(ierr); 405 ierr = MatGetOwnershipRanges(A,&ranges);CHKERRQ(ierr); 406 for (PetscInt i=0; i<size; ++i) { 407 offset[2*i] = ranges[i]; 408 offset[2*i+1] = ranges[i+1] - ranges[i]; 409 } 410 switch (a->clustering) { 411 case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC: 412 t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim); 413 break; 414 case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC: 415 t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim); 416 break; 417 case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR: 418 t = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim); 419 break; 420 default: 421 t = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim); 422 } 423 t->set_minclustersize(a->bs[0]); 424 t->build(A->rmap->N,a->gcoords_target,offset); 425 if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N,A->cmap->N,a->dim,a->kernel,a->kernelctx); 426 else { 427 a->wrapper = NULL; 428 generator = reinterpret_cast<htool::IMatrix<PetscScalar>*>(a->kernelctx); 429 } 430 if (a->gcoords_target != a->gcoords_source) { 431 ierr = MatGetOwnershipRangesColumn(A,&ranges);CHKERRQ(ierr); 432 for (PetscInt i=0; i<size; ++i) { 433 offset[2*i] = ranges[i]; 434 offset[2*i+1] = ranges[i+1] - ranges[i]; 435 } 436 switch (a->clustering) { 437 case MAT_HTOOL_CLUSTERING_PCA_GEOMETRIC: 438 s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::GeometricSplitting>>>(a->dim); 439 break; 440 case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_GEOMETRIC: 441 s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::GeometricSplitting>>>(a->dim); 442 break; 443 case MAT_HTOOL_CLUSTERING_BOUNDING_BOX_1_REGULAR: 444 s = std::make_shared<htool::Cluster<htool::BoundingBox1<htool::SplittingTypes::RegularSplitting>>>(a->dim); 445 break; 446 default: 447 s = std::make_shared<htool::Cluster<htool::PCA<htool::SplittingTypes::RegularSplitting>>>(a->dim); 448 } 449 s->set_minclustersize(a->bs[0]); 450 s->build(A->cmap->N,a->gcoords_source,offset); 451 S = uplo = 'N'; 452 } 453 ierr = PetscFree(offset);CHKERRQ(ierr); 454 switch (a->compressor) { 455 case MAT_HTOOL_COMPRESSOR_FULL_ACA: 456 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)); 457 break; 458 case MAT_HTOOL_COMPRESSOR_SVD: 459 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)); 460 break; 461 default: 462 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)); 463 } 464 a->hmatrix->set_maxblocksize(a->bs[1]); 465 a->hmatrix->set_mintargetdepth(a->depth[0]); 466 a->hmatrix->set_minsourcedepth(a->depth[1]); 467 if (s) a->hmatrix->build_auto(a->wrapper ? *a->wrapper : *generator,a->gcoords_target,a->gcoords_source); 468 else a->hmatrix->build_auto_sym(a->wrapper ? *a->wrapper : *generator,a->gcoords_target); 469 PetscFunctionReturn(0); 470 } 471 472 static PetscErrorCode MatProductNumeric_Htool(Mat C) 473 { 474 Mat_Product *product = C->product; 475 Mat_Htool *a = (Mat_Htool*)product->A->data; 476 const PetscScalar *in; 477 PetscScalar *out; 478 PetscInt lda; 479 PetscErrorCode ierr; 480 481 PetscFunctionBegin; 482 MatCheckProduct(C,1); 483 switch (product->type) { 484 case MATPRODUCT_AB: 485 PetscInt N; 486 ierr = MatGetSize(C,NULL,&N);CHKERRQ(ierr); 487 ierr = MatDenseGetLDA(C,&lda);CHKERRQ(ierr); 488 if (lda != C->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported leading dimension (%D != %D)",lda,C->rmap->n); 489 ierr = MatDenseGetArrayRead(product->B,&in);CHKERRQ(ierr); 490 ierr = MatDenseGetArrayWrite(C,&out);CHKERRQ(ierr); 491 a->hmatrix->mvprod_local_to_local(in,out,N); 492 ierr = MatDenseRestoreArrayWrite(C,&out);CHKERRQ(ierr); 493 ierr = MatDenseRestoreArrayRead(product->B,&in);CHKERRQ(ierr); 494 ierr = MatScale(C,a->s);CHKERRQ(ierr); 495 break; 496 default: 497 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductType %s is not supported",MatProductTypes[product->type]); 498 } 499 PetscFunctionReturn(0); 500 } 501 502 static PetscErrorCode MatProductSymbolic_Htool(Mat C) 503 { 504 Mat_Product *product = C->product; 505 Mat A,B; 506 PetscBool flg; 507 PetscErrorCode ierr; 508 509 PetscFunctionBegin; 510 MatCheckProduct(C,1); 511 A = product->A; 512 B = product->B; 513 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 514 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"MatProduct_AB not supported for %s",((PetscObject)product->B)->type_name); 515 switch (product->type) { 516 case MATPRODUCT_AB: 517 if (C->rmap->n == PETSC_DECIDE || C->cmap->n == PETSC_DECIDE || C->rmap->N == PETSC_DECIDE || C->cmap->N == PETSC_DECIDE) { 518 ierr = MatSetSizes(C,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);CHKERRQ(ierr); 519 } 520 ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr); 521 ierr = MatSetUp(C);CHKERRQ(ierr); 522 ierr = MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 523 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 524 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 525 break; 526 default: 527 SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[product->type]); 528 } 529 C->ops->productsymbolic = NULL; 530 C->ops->productnumeric = MatProductNumeric_Htool; 531 PetscFunctionReturn(0); 532 } 533 534 static PetscErrorCode MatProductSetFromOptions_Htool(Mat C) 535 { 536 PetscFunctionBegin; 537 MatCheckProduct(C,1); 538 if (C->product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_Htool; 539 PetscFunctionReturn(0); 540 } 541 542 static PetscErrorCode MatHtoolGetHierarchicalMat_Htool(Mat A,const htool::VirtualHMatrix<PetscScalar> **hmatrix) 543 { 544 Mat_Htool *a = (Mat_Htool*)A->data; 545 546 PetscFunctionBegin; 547 *hmatrix = a->hmatrix; 548 PetscFunctionReturn(0); 549 } 550 551 /*@C 552 MatHtoolGetHierarchicalMat - Retrieves the opaque pointer to a Htool virtual matrix stored in a MATHTOOL. 553 554 Input Parameter: 555 . A - hierarchical matrix 556 557 Output Parameter: 558 . hmatrix - opaque pointer to a Htool virtual matrix 559 560 Level: advanced 561 562 .seealso: MATHTOOL 563 @*/ 564 PETSC_EXTERN PetscErrorCode MatHtoolGetHierarchicalMat(Mat A,const htool::VirtualHMatrix<PetscScalar> **hmatrix) 565 { 566 PetscErrorCode ierr; 567 568 PetscFunctionBegin; 569 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 570 PetscValidPointer(hmatrix,2); 571 ierr = PetscTryMethod(A,"MatHtoolGetHierarchicalMat_C",(Mat,const htool::VirtualHMatrix<PetscScalar>**),(A,hmatrix));CHKERRQ(ierr); 572 PetscFunctionReturn(0); 573 } 574 575 static PetscErrorCode MatHtoolSetKernel_Htool(Mat A,MatHtoolKernel kernel,void *kernelctx) 576 { 577 Mat_Htool *a = (Mat_Htool*)A->data; 578 579 PetscFunctionBegin; 580 a->kernel = kernel; 581 a->kernelctx = kernelctx; 582 delete a->wrapper; 583 if (a->kernel) a->wrapper = new WrapperHtool(A->rmap->N,A->cmap->N,a->dim,a->kernel,a->kernelctx); 584 PetscFunctionReturn(0); 585 } 586 587 /*@C 588 MatHtoolSetKernel - Sets the kernel and context used for the assembly of a MATHTOOL. 589 590 Input Parameters: 591 + A - hierarchical matrix 592 . kernel - computational kernel (or NULL) 593 - kernelctx - kernel context (if kernel is NULL, the pointer must be of type htool::IMatrix<PetscScalar>*) 594 595 Level: advanced 596 597 .seealso: MATHTOOL, MatCreateHtoolFromKernel() 598 @*/ 599 PETSC_EXTERN PetscErrorCode MatHtoolSetKernel(Mat A,MatHtoolKernel kernel,void *kernelctx) 600 { 601 PetscErrorCode ierr; 602 603 PetscFunctionBegin; 604 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 605 if (!kernelctx) PetscValidFunction(kernel,2); 606 if (!kernel) PetscValidPointer(kernelctx,3); 607 ierr = PetscTryMethod(A,"MatHtoolSetKernel_C",(Mat,MatHtoolKernel,void*),(A,kernel,kernelctx));CHKERRQ(ierr); 608 PetscFunctionReturn(0); 609 } 610 611 static PetscErrorCode MatHtoolGetPermutationSource_Htool(Mat A,IS* is) 612 { 613 Mat_Htool *a = (Mat_Htool*)A->data; 614 std::vector<PetscInt> source; 615 PetscErrorCode ierr; 616 617 PetscFunctionBegin; 618 source = a->hmatrix->get_local_perm_source(); 619 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),source.size(),source.data(),PETSC_COPY_VALUES,is);CHKERRQ(ierr); 620 ierr = ISSetPermutation(*is);CHKERRQ(ierr); 621 PetscFunctionReturn(0); 622 } 623 624 /*@C 625 MatHtoolGetPermutationSource - Gets the permutation associated to the source cluster. 626 627 Input Parameter: 628 . A - hierarchical matrix 629 630 Output Parameter: 631 . is - permutation 632 633 Level: advanced 634 635 .seealso: MATHTOOL, MatHtoolGetPermutationTarget(), MatHtoolUsePermutation() 636 @*/ 637 PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationSource(Mat A,IS* is) 638 { 639 PetscErrorCode ierr; 640 641 PetscFunctionBegin; 642 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 643 if (!is) PetscValidPointer(is,2); 644 ierr = PetscTryMethod(A,"MatHtoolGetPermutationSource_C",(Mat,IS*),(A,is));CHKERRQ(ierr); 645 PetscFunctionReturn(0); 646 } 647 648 static PetscErrorCode MatHtoolGetPermutationTarget_Htool(Mat A,IS* is) 649 { 650 Mat_Htool *a = (Mat_Htool*)A->data; 651 std::vector<PetscInt> target; 652 PetscErrorCode ierr; 653 654 PetscFunctionBegin; 655 target = a->hmatrix->get_local_perm_target(); 656 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),target.size(),target.data(),PETSC_COPY_VALUES,is);CHKERRQ(ierr); 657 ierr = ISSetPermutation(*is);CHKERRQ(ierr); 658 PetscFunctionReturn(0); 659 } 660 661 /*@C 662 MatHtoolGetPermutationTarget - Gets the permutation associated to the target cluster. 663 664 Input Parameter: 665 . A - hierarchical matrix 666 667 Output Parameter: 668 . is - permutation 669 670 Level: advanced 671 672 .seealso: MATHTOOL, MatHtoolGetPermutationSource(), MatHtoolUsePermutation() 673 @*/ 674 PETSC_EXTERN PetscErrorCode MatHtoolGetPermutationTarget(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,"MatHtoolGetPermutationTarget_C",(Mat,IS*),(A,is));CHKERRQ(ierr); 682 PetscFunctionReturn(0); 683 } 684 685 static PetscErrorCode MatHtoolUsePermutation_Htool(Mat A,PetscBool use) 686 { 687 Mat_Htool *a = (Mat_Htool*)A->data; 688 689 PetscFunctionBegin; 690 a->hmatrix->set_use_permutation(use); 691 PetscFunctionReturn(0); 692 } 693 694 /*@C 695 MatHtoolUsePermutation - Sets whether MATHTOOL should permute input (resp. output) vectors following its internal source (resp. target) permutation. 696 697 Input Parameters: 698 + A - hierarchical matrix 699 - use - Boolean value 700 701 Level: advanced 702 703 .seealso: MATHTOOL, MatHtoolGetPermutationSource(), MatHtoolGetPermutationTarget() 704 @*/ 705 PETSC_EXTERN PetscErrorCode MatHtoolUsePermutation(Mat A,PetscBool use) 706 { 707 PetscErrorCode ierr; 708 709 PetscFunctionBegin; 710 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 711 PetscValidLogicalCollectiveBool(A,use,2); 712 ierr = PetscTryMethod(A,"MatHtoolUsePermutation_C",(Mat,PetscBool),(A,use));CHKERRQ(ierr); 713 PetscFunctionReturn(0); 714 } 715 716 static PetscErrorCode MatConvert_Htool_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *B) 717 { 718 Mat C; 719 Mat_Htool *a = (Mat_Htool*)A->data; 720 PetscInt lda; 721 PetscScalar *array; 722 PetscErrorCode ierr; 723 724 PetscFunctionBegin; 725 if (reuse == MAT_REUSE_MATRIX) { 726 C = *B; 727 if (C->rmap->n != A->rmap->n || C->cmap->N != A->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible dimensions"); 728 ierr = MatDenseGetLDA(C,&lda);CHKERRQ(ierr); 729 if (lda != C->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported leading dimension (%D != %D)",lda,C->rmap->n); 730 } else { 731 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 732 ierr = MatSetSizes(C,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 733 ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr); 734 ierr = MatSetUp(C);CHKERRQ(ierr); 735 ierr = MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 736 } 737 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 738 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 739 ierr = MatDenseGetArrayWrite(C,&array);CHKERRQ(ierr); 740 a->hmatrix->copy_local_dense_perm(array); 741 ierr = MatDenseRestoreArrayWrite(C,&array);CHKERRQ(ierr); 742 ierr = MatScale(C,a->s);CHKERRQ(ierr); 743 if (reuse == MAT_INPLACE_MATRIX) { 744 ierr = MatHeaderReplace(A,&C);CHKERRQ(ierr); 745 } else *B = C; 746 PetscFunctionReturn(0); 747 } 748 749 static PetscErrorCode GenEntriesTranspose(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt *rows,const PetscInt *cols,PetscScalar *ptr,void *ctx) 750 { 751 MatHtoolKernelTranspose *generator = (MatHtoolKernelTranspose*)ctx; 752 PetscScalar *tmp; 753 PetscErrorCode ierr; 754 755 PetscFunctionBegin; 756 generator->kernel(sdim,N,M,cols,rows,ptr,generator->kernelctx); 757 ierr = PetscMalloc1(M*N,&tmp);CHKERRQ(ierr); 758 ierr = PetscArraycpy(tmp,ptr,M*N);CHKERRQ(ierr); 759 for (PetscInt i=0; i<M; ++i) { 760 for (PetscInt j=0; j<N; ++j) ptr[i+j*M] = tmp[j+i*N]; 761 } 762 ierr = PetscFree(tmp);CHKERRQ(ierr); 763 PetscFunctionReturn(0); 764 } 765 766 /* naive implementation which keeps a reference to the original Mat */ 767 static PetscErrorCode MatTranspose_Htool(Mat A,MatReuse reuse,Mat *B) 768 { 769 Mat C; 770 Mat_Htool *a = (Mat_Htool*)A->data,*c; 771 PetscInt M = A->rmap->N,N = A->cmap->N,m = A->rmap->n,n = A->cmap->n; 772 PetscContainer container; 773 MatHtoolKernelTranspose *kernelt; 774 PetscErrorCode ierr; 775 776 PetscFunctionBegin; 777 if (reuse == MAT_INPLACE_MATRIX) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatTranspose() with MAT_INPLACE_MATRIX not supported"); 778 if (reuse == MAT_INITIAL_MATRIX) { 779 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 780 ierr = MatSetSizes(C,n,m,N,M);CHKERRQ(ierr); 781 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 782 ierr = MatSetUp(C);CHKERRQ(ierr); 783 ierr = PetscContainerCreate(PetscObjectComm((PetscObject)C),&container);CHKERRQ(ierr); 784 ierr = PetscNew(&kernelt);CHKERRQ(ierr); 785 ierr = PetscContainerSetPointer(container,kernelt);CHKERRQ(ierr); 786 ierr = PetscObjectCompose((PetscObject)C,"KernelTranspose",(PetscObject)container);CHKERRQ(ierr); 787 } else { 788 C = *B; 789 ierr = PetscObjectQuery((PetscObject)C,"KernelTranspose",(PetscObject*)&container);CHKERRQ(ierr); 790 if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatTranspose() with MAT_INITIAL_MATRIX first"); 791 ierr = PetscContainerGetPointer(container,(void**)&kernelt);CHKERRQ(ierr); 792 } 793 c = (Mat_Htool*)C->data; 794 c->dim = a->dim; 795 c->s = a->s; 796 c->kernel = GenEntriesTranspose; 797 if (kernelt->A != A) { 798 ierr = MatDestroy(&kernelt->A);CHKERRQ(ierr); 799 kernelt->A = A; 800 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 801 } 802 kernelt->kernel = a->kernel; 803 kernelt->kernelctx = a->kernelctx; 804 c->kernelctx = kernelt; 805 if (reuse == MAT_INITIAL_MATRIX) { 806 ierr = PetscMalloc1(N*c->dim,&c->gcoords_target);CHKERRQ(ierr); 807 ierr = PetscArraycpy(c->gcoords_target,a->gcoords_source,N*c->dim);CHKERRQ(ierr); 808 if (a->gcoords_target != a->gcoords_source) { 809 ierr = PetscMalloc1(M*c->dim,&c->gcoords_source);CHKERRQ(ierr); 810 ierr = PetscArraycpy(c->gcoords_source,a->gcoords_target,M*c->dim);CHKERRQ(ierr); 811 } else c->gcoords_source = c->gcoords_target; 812 ierr = PetscCalloc2(M,&c->work_source,N,&c->work_target);CHKERRQ(ierr); 813 } 814 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 815 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 816 if (reuse == MAT_INITIAL_MATRIX) *B = C; 817 PetscFunctionReturn(0); 818 } 819 820 /*@C 821 MatCreateHtoolFromKernel - Creates a MATHTOOL from a user-supplied kernel. 822 823 Input Parameters: 824 + comm - MPI communicator 825 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 826 . n - number of local columns (or PETSC_DECIDE to have calculated if N is given) 827 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 828 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 829 . spacedim - dimension of the space coordinates 830 . coords_target - coordinates of the target 831 . coords_source - coordinates of the source 832 . kernel - computational kernel (or NULL) 833 - kernelctx - kernel context (if kernel is NULL, the pointer must be of type htool::IMatrix<PetscScalar>*) 834 835 Output Parameter: 836 . B - matrix 837 838 Options Database Keys: 839 + -mat_htool_min_cluster_size <PetscInt> - minimal leaf size in cluster tree 840 . -mat_htool_max_block_size <PetscInt> - maximal number of coefficients in a dense block 841 . -mat_htool_epsilon <PetscReal> - relative error in Frobenius norm when approximating a block 842 . -mat_htool_eta <PetscReal> - admissibility condition tolerance 843 . -mat_htool_min_target_depth <PetscInt> - minimal cluster tree depth associated with the rows 844 . -mat_htool_min_source_depth <PetscInt> - minimal cluster tree depth associated with the columns 845 . -mat_htool_compressor <sympartialACA, fullACA, SVD> - type of compression 846 - -mat_htool_clustering <PCARegular, PCAGeometric, BounbingBox1Regular, BoundingBox1Geometric> - type of clustering 847 848 Level: intermediate 849 850 .seealso: MatCreate(), MATHTOOL, PCSetCoordinates(), MatHtoolSetKernel(), MatHtoolCompressorType, MATHARA, MatCreateHaraFromKernel() 851 @*/ 852 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) 853 { 854 Mat A; 855 Mat_Htool *a; 856 PetscErrorCode ierr; 857 858 PetscFunctionBegin; 859 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 860 PetscValidLogicalCollectiveInt(A,spacedim,6); 861 PetscValidRealPointer(coords_target,7); 862 PetscValidRealPointer(coords_source,8); 863 if (!kernelctx) PetscValidFunction(kernel,9); 864 if (!kernel) PetscValidPointer(kernelctx,10); 865 ierr = MatSetSizes(A,m,n,M,N);CHKERRQ(ierr); 866 ierr = MatSetType(A,MATHTOOL);CHKERRQ(ierr); 867 ierr = MatSetUp(A);CHKERRQ(ierr); 868 a = (Mat_Htool*)A->data; 869 a->dim = spacedim; 870 a->s = 1.0; 871 a->kernel = kernel; 872 a->kernelctx = kernelctx; 873 ierr = PetscCalloc1(A->rmap->N*spacedim,&a->gcoords_target);CHKERRQ(ierr); 874 ierr = PetscArraycpy(a->gcoords_target+A->rmap->rstart*spacedim,coords_target,A->rmap->n*spacedim);CHKERRQ(ierr); 875 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 */ 876 if (coords_target != coords_source) { 877 ierr = PetscCalloc1(A->cmap->N*spacedim,&a->gcoords_source);CHKERRQ(ierr); 878 ierr = PetscArraycpy(a->gcoords_source+A->cmap->rstart*spacedim,coords_source,A->cmap->n*spacedim);CHKERRQ(ierr); 879 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 */ 880 } else a->gcoords_source = a->gcoords_target; 881 ierr = PetscCalloc2(A->cmap->N,&a->work_source,A->rmap->N,&a->work_target);CHKERRQ(ierr); 882 *B = A; 883 PetscFunctionReturn(0); 884 } 885 886 /*MC 887 MATHTOOL = "htool" - A matrix type for hierarchical matrices using the Htool package. 888 889 Use ./configure --download-htool to install PETSc to use Htool. 890 891 Options Database Keys: 892 . -mat_type htool - matrix type to "htool" during a call to MatSetFromOptions() 893 894 Level: beginner 895 896 .seealso: MATHARA, MATDENSE, MatCreateHtoolFromKernel(), MatHtoolSetKernel() 897 M*/ 898 PETSC_EXTERN PetscErrorCode MatCreate_Htool(Mat A) 899 { 900 Mat_Htool *a; 901 PetscErrorCode ierr; 902 903 PetscFunctionBegin; 904 ierr = PetscNewLog(A,&a);CHKERRQ(ierr); 905 A->data = (void*)a; 906 ierr = PetscObjectChangeTypeName((PetscObject)A,MATHTOOL);CHKERRQ(ierr); 907 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 908 A->ops->getdiagonal = MatGetDiagonal_Htool; 909 A->ops->getdiagonalblock = MatGetDiagonalBlock_Htool; 910 A->ops->mult = MatMult_Htool; 911 A->ops->multadd = MatMultAdd_Htool; 912 A->ops->increaseoverlap = MatIncreaseOverlap_Htool; 913 A->ops->createsubmatrices = MatCreateSubMatrices_Htool; 914 A->ops->transpose = MatTranspose_Htool; 915 A->ops->destroy = MatDestroy_Htool; 916 A->ops->view = MatView_Htool; 917 A->ops->setfromoptions = MatSetFromOptions_Htool; 918 A->ops->scale = MatScale_Htool; 919 A->ops->getrow = MatGetRow_Htool; 920 A->ops->restorerow = MatRestoreRow_Htool; 921 A->ops->assemblyend = MatAssemblyEnd_Htool; 922 a->dim = 0; 923 a->gcoords_target = NULL; 924 a->gcoords_source = NULL; 925 a->s = 1.0; 926 a->bs[0] = 10; 927 a->bs[1] = 1000000; 928 a->epsilon = PetscSqrtReal(PETSC_SMALL); 929 a->eta = 10.0; 930 a->depth[0] = 0; 931 a->depth[1] = 0; 932 a->compressor = MAT_HTOOL_COMPRESSOR_SYMPARTIAL_ACA; 933 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_seqdense_C",MatProductSetFromOptions_Htool);CHKERRQ(ierr); 934 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_htool_mpidense_C",MatProductSetFromOptions_Htool);CHKERRQ(ierr); 935 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_seqdense_C",MatConvert_Htool_Dense);CHKERRQ(ierr); 936 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_htool_mpidense_C",MatConvert_Htool_Dense);CHKERRQ(ierr); 937 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetHierarchicalMat_C",MatHtoolGetHierarchicalMat_Htool);CHKERRQ(ierr); 938 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHtoolSetKernel_C",MatHtoolSetKernel_Htool);CHKERRQ(ierr); 939 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationSource_C",MatHtoolGetPermutationSource_Htool);CHKERRQ(ierr); 940 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHtoolGetPermutationTarget_C",MatHtoolGetPermutationTarget_Htool);CHKERRQ(ierr); 941 ierr = PetscObjectComposeFunction((PetscObject)A,"MatHtoolUsePermutation_C",MatHtoolUsePermutation_Htool);CHKERRQ(ierr); 942 PetscFunctionReturn(0); 943 } 944