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