Lines Matching +full:- +full:m
13 Mat M;
16 /* sampler for Newton-Schultz */
31 /* Newton-Schultz customization */
54 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
59 if (pch2opus->sdim) PetscFunctionReturn(PETSC_SUCCESS);
61 if (!dm) PetscCall(MatGetDM(pc->useAmat ? pc->mat : pc->pmat, &dm));
80 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
83 pch2opus->sdim = 0;
84 pch2opus->nlocc = 0;
85 PetscCall(PetscFree(pch2opus->coords));
91 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
95 if (pch2opus->sdim && sdim == pch2opus->sdim && nlocc == pch2opus->nlocc) {
96 PetscCall(PetscArraycmp(pch2opus->coords, coords, sdim * nlocc, &reset));
102 PetscCall(PetscMalloc1(sdim * nlocc, &pch2opus->coords));
103 PetscCall(PetscArraycpy(pch2opus->coords, coords, sdim * nlocc));
104 pch2opus->sdim = sdim;
105 pch2opus->nlocc = nlocc;
112 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
116 PetscCall(MatDestroy(&pch2opus->A));
117 PetscCall(MatDestroy(&pch2opus->M));
118 PetscCall(MatDestroy(&pch2opus->T));
119 PetscCall(VecDestroy(&pch2opus->w));
120 PetscCall(MatDestroy(&pch2opus->S));
121 PetscCall(VecDestroy(&pch2opus->wns[0]));
122 PetscCall(VecDestroy(&pch2opus->wns[1]));
123 PetscCall(VecDestroy(&pch2opus->wns[2]));
124 PetscCall(VecDestroy(&pch2opus->wns[3]));
125 PetscCall(MatDestroy(&pch2opus->wnsmat[0]));
126 PetscCall(MatDestroy(&pch2opus->wnsmat[1]));
127 PetscCall(MatDestroy(&pch2opus->wnsmat[2]));
128 PetscCall(MatDestroy(&pch2opus->wnsmat[3]));
130 PetscCall(PetscFree(pc->data));
136 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
140 PetscCall(PetscOptionsInt("-pc_h2opus_maxits", "Maximum number of iterations for Newton-Schultz", NULL, pch2opus->maxits, &pch2opus->maxits, NULL));
141 PetscCall(PetscOptionsBool("-pc_h2opus_monitor", "Monitor Newton-Schultz convergence", NULL, pch2opus->monitor, &pch2opus->monitor, NULL));
142 PetscCall(PetscOptionsReal("-pc_h2opus_atol", "Absolute tolerance", NULL, pch2opus->atol, &pch2opus->atol, NULL));
143 PetscCall(PetscOptionsReal("-pc_h2opus_rtol", "Relative tolerance", NULL, pch2opus->rtol, &pch2opus->rtol, NULL));
144 PetscCall(PetscOptionsEnum("-pc_h2opus_norm_type", "Norm type for convergence monitoring", NULL, NormTypes, (PetscEnum)pch2opus->normtype, (PetscEnum *)&pch2opus->normtype, NULL));
145 PetscCall(PetscOptionsInt("-pc_h2opus_hyperorder", "Hyper power order of sampling", NULL, pch2opus->hyperorder, &pch2opus->hyperorder, NULL));
146 PetscCall(PetscOptionsInt("-pc_h2opus_leafsize", "Leaf size when constructed from kernel", NULL, pch2opus->leafsize, &pch2opus->leafsize, NULL));
147 PetscCall(PetscOptionsReal("-pc_h2opus_eta", "Admissibility condition tolerance", NULL, pch2opus->eta, &pch2opus->eta, NULL));
148 PetscCall(PetscOptionsInt("-pc_h2opus_maxrank", "Maximum rank when constructed from matvecs", NULL, pch2opus->max_rank, &pch2opus->max_rank, NULL));
149 PetscCall(PetscOptionsInt("-pc_h2opus_samples", "Number of samples to be taken concurrently when constructing from matvecs", NULL, pch2opus->bs, &pch2opus->bs, NULL));
150 PetscCall(PetscOptionsReal("-pc_h2opus_mrtol", "Relative tolerance for construction from sampling", NULL, pch2opus->mrtol, &pch2opus->mrtol, NULL));
151 PetscCall(PetscOptionsBool("-pc_h2opus_forcecpu", "Force construction on CPU", NULL, pch2opus->forcecpu, &pch2opus->forcecpu, NULL));
158 Mat M;
168 /* PetscCall(MatMultTranspose(aat->M,x,aat->w)); */
169 PetscCall(MatMultTranspose(aat->A, x, aat->w));
170 PetscCall(MatMult(aat->A, aat->w, y));
176 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
177 Mat A = pc->useAmat ? pc->mat : pc->pmat, AAt;
178 PetscInt M, m;
185 aat.M = pch2opus->M; /* unused so far */
187 PetscCall(MatGetSize(A, &M, NULL));
188 PetscCall(MatGetLocalSize(A, &m, NULL));
189 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A), m, m, M, M, &aat, &AAt));
190 PetscCall(MatBindToCPU(AAt, pch2opus->boundtocpu));
197 pch2opus->s0 = 1. / n;
205 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
208 if (t) PetscCall(MatMultTranspose(pch2opus->M, x, y));
209 else PetscCall(MatMult(pch2opus->M, x, y));
215 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
218 if (t) PetscCall(MatTransposeMatMult(pch2opus->M, X, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Y));
219 else PetscCall(MatMatMult(pch2opus->M, X, MAT_REUSE_MATRIX, PETSC_CURRENT, &Y));
251 /* used to test the norm of (M^-1 A - I) */
252 static PetscErrorCode MatMultKernel_MAmI(Mat M, Vec x, Vec y, PetscBool t)
260 PetscCall(MatShellGetContext(M, &pc));
261 pch2opus = (PC_H2OPUS *)pc->data;
262 if (!pch2opus->w) PetscCall(MatCreateVecs(pch2opus->M, &pch2opus->w, NULL));
263 A = pch2opus->A;
264 PetscCall(VecBindToCPU(pch2opus->w, pch2opus->boundtocpu));
267 PetscCall(PCApplyTranspose_H2OPUS(pc, x, pch2opus->w));
268 PetscCall(MatMultTranspose(A, pch2opus->w, y));
270 PetscCall(MatMultTranspose(A, x, pch2opus->w));
271 PetscCall(PCApplyTranspose_H2OPUS(pc, pch2opus->w, y));
275 PetscCall(MatMult(A, x, pch2opus->w));
276 PetscCall(PCApply_H2OPUS(pc, pch2opus->w, y));
278 PetscCall(PCApply_H2OPUS(pc, x, pch2opus->w));
279 PetscCall(MatMult(A, pch2opus->w, y));
282 PetscCall(VecAXPY(y, -1.0, x));
302 for i = 1 . . . l - 1 do
303 R = (I - A * Xk) * R
307 static PetscErrorCode MatMultKernel_Hyper(Mat M, Vec x, Vec y, PetscBool t)
314 PetscCall(MatShellGetContext(M, &pc));
315 pch2opus = (PC_H2OPUS *)pc->data;
316 A = pch2opus->A;
317 PetscCall(MatCreateVecs(pch2opus->M, pch2opus->wns[0] ? NULL : &pch2opus->wns[0], pch2opus->wns[1] ? NULL : &pch2opus->wns[1]));
318 PetscCall(MatCreateVecs(pch2opus->M, pch2opus->wns[2] ? NULL : &pch2opus->wns[2], pch2opus->wns[3] ? NULL : &pch2opus->wns[3]));
319 PetscCall(VecBindToCPU(pch2opus->wns[0], pch2opus->boundtocpu));
320 PetscCall(VecBindToCPU(pch2opus->wns[1], pch2opus->boundtocpu));
321 PetscCall(VecBindToCPU(pch2opus->wns[2], pch2opus->boundtocpu));
322 PetscCall(VecBindToCPU(pch2opus->wns[3], pch2opus->boundtocpu));
323 PetscCall(VecCopy(x, pch2opus->wns[0]));
324 PetscCall(VecCopy(x, pch2opus->wns[3]));
326 for (PetscInt i = 0; i < pch2opus->hyperorder - 1; i++) {
327 PetscCall(MatMultTranspose(A, pch2opus->wns[0], pch2opus->wns[1]));
328 PetscCall(PCApplyTranspose_H2OPUS(pc, pch2opus->wns[1], pch2opus->wns[2]));
329 PetscCall(VecAXPY(pch2opus->wns[0], -1., pch2opus->wns[2]));
330 PetscCall(VecAXPY(pch2opus->wns[3], 1., pch2opus->wns[0]));
332 PetscCall(PCApplyTranspose_H2OPUS(pc, pch2opus->wns[3], y));
334 for (PetscInt i = 0; i < pch2opus->hyperorder - 1; i++) {
335 PetscCall(PCApply_H2OPUS(pc, pch2opus->wns[0], pch2opus->wns[1]));
336 PetscCall(MatMult(A, pch2opus->wns[1], pch2opus->wns[2]));
337 PetscCall(VecAXPY(pch2opus->wns[0], -1., pch2opus->wns[2]));
338 PetscCall(VecAXPY(pch2opus->wns[3], 1., pch2opus->wns[0]));
340 PetscCall(PCApply_H2OPUS(pc, pch2opus->wns[3], y));
345 static PetscErrorCode MatMult_Hyper(Mat M, Vec x, Vec y)
348 PetscCall(MatMultKernel_Hyper(M, x, y, PETSC_FALSE));
352 static PetscErrorCode MatMultTranspose_Hyper(Mat M, Vec x, Vec y)
355 PetscCall(MatMultKernel_Hyper(M, x, y, PETSC_TRUE));
360 static PetscErrorCode MatMatMultKernel_Hyper(Mat M, Mat X, Mat Y, PetscBool t)
368 PetscCall(MatShellGetContext(M, &pc));
369 pch2opus = (PC_H2OPUS *)pc->data;
370 A = pch2opus->A;
371 if (pch2opus->wnsmat[0] && pch2opus->wnsmat[0]->cmap->N != X->cmap->N) {
372 PetscCall(MatDestroy(&pch2opus->wnsmat[0]));
373 PetscCall(MatDestroy(&pch2opus->wnsmat[1]));
375 if (!pch2opus->wnsmat[0]) {
376 PetscCall(MatDuplicate(X, MAT_SHARE_NONZERO_PATTERN, &pch2opus->wnsmat[0]));
377 PetscCall(MatDuplicate(Y, MAT_SHARE_NONZERO_PATTERN, &pch2opus->wnsmat[1]));
379 if (pch2opus->wnsmat[2] && pch2opus->wnsmat[2]->cmap->N != X->cmap->N) {
380 PetscCall(MatDestroy(&pch2opus->wnsmat[2]));
381 PetscCall(MatDestroy(&pch2opus->wnsmat[3]));
383 if (!pch2opus->wnsmat[2]) {
384 PetscCall(MatDuplicate(X, MAT_SHARE_NONZERO_PATTERN, &pch2opus->wnsmat[2]));
385 PetscCall(MatDuplicate(Y, MAT_SHARE_NONZERO_PATTERN, &pch2opus->wnsmat[3]));
387 PetscCall(MatCopy(X, pch2opus->wnsmat[0], SAME_NONZERO_PATTERN));
388 PetscCall(MatCopy(X, pch2opus->wnsmat[3], SAME_NONZERO_PATTERN));
390 for (i = 0; i < pch2opus->hyperorder - 1; i++) {
391 PetscCall(MatTransposeMatMult(A, pch2opus->wnsmat[0], MAT_REUSE_MATRIX, PETSC_CURRENT, &pch2opus->wnsmat[1]));
392 PetscCall(PCApplyTransposeMat_H2OPUS(pc, pch2opus->wnsmat[1], pch2opus->wnsmat[2]));
393 PetscCall(MatAXPY(pch2opus->wnsmat[0], -1., pch2opus->wnsmat[2], SAME_NONZERO_PATTERN));
394 PetscCall(MatAXPY(pch2opus->wnsmat[3], 1., pch2opus->wnsmat[0], SAME_NONZERO_PATTERN));
396 PetscCall(PCApplyTransposeMat_H2OPUS(pc, pch2opus->wnsmat[3], Y));
398 for (i = 0; i < pch2opus->hyperorder - 1; i++) {
399 PetscCall(PCApplyMat_H2OPUS(pc, pch2opus->wnsmat[0], pch2opus->wnsmat[1]));
400 PetscCall(MatMatMult(A, pch2opus->wnsmat[1], MAT_REUSE_MATRIX, PETSC_CURRENT, &pch2opus->wnsmat[2]));
401 PetscCall(MatAXPY(pch2opus->wnsmat[0], -1., pch2opus->wnsmat[2], SAME_NONZERO_PATTERN));
402 PetscCall(MatAXPY(pch2opus->wnsmat[3], 1., pch2opus->wnsmat[0], SAME_NONZERO_PATTERN));
404 PetscCall(PCApplyMat_H2OPUS(pc, pch2opus->wnsmat[3], Y));
409 static PetscErrorCode MatMatMultNumeric_Hyper(Mat M, Mat X, Mat Y, PetscCtx ctx)
412 PetscCall(MatMatMultKernel_Hyper(M, X, Y, PETSC_FALSE));
416 /* Basic Newton-Schultz sampler: (2 * I - M * A)*M */
417 static PetscErrorCode MatMultKernel_NS(Mat M, Vec x, Vec y, PetscBool t)
424 PetscCall(MatShellGetContext(M, &pc));
425 pch2opus = (PC_H2OPUS *)pc->data;
426 A = pch2opus->A;
427 PetscCall(MatCreateVecs(pch2opus->M, pch2opus->wns[0] ? NULL : &pch2opus->wns[0], pch2opus->wns[1] ? NULL : &pch2opus->wns[1]));
428 PetscCall(VecBindToCPU(pch2opus->wns[0], pch2opus->boundtocpu));
429 PetscCall(VecBindToCPU(pch2opus->wns[1], pch2opus->boundtocpu));
432 PetscCall(MatMultTranspose(A, y, pch2opus->wns[1]));
433 PetscCall(PCApplyTranspose_H2OPUS(pc, pch2opus->wns[1], pch2opus->wns[0]));
434 PetscCall(VecAXPBY(y, -1., 2., pch2opus->wns[0]));
437 PetscCall(MatMult(A, y, pch2opus->wns[0]));
438 PetscCall(PCApply_H2OPUS(pc, pch2opus->wns[0], pch2opus->wns[1]));
439 PetscCall(VecAXPBY(y, -1., 2., pch2opus->wns[1]));
444 static PetscErrorCode MatMult_NS(Mat M, Vec x, Vec y)
447 PetscCall(MatMultKernel_NS(M, x, y, PETSC_FALSE));
451 static PetscErrorCode MatMultTranspose_NS(Mat M, Vec x, Vec y)
454 PetscCall(MatMultKernel_NS(M, x, y, PETSC_TRUE));
458 /* Basic Newton-Schultz sampler: (2 * I - M * A)*M, MatMat version */
459 static PetscErrorCode MatMatMultKernel_NS(Mat M, Mat X, Mat Y, PetscBool t)
466 PetscCall(MatShellGetContext(M, &pc));
467 pch2opus = (PC_H2OPUS *)pc->data;
468 A = pch2opus->A;
469 if (pch2opus->wnsmat[0] && pch2opus->wnsmat[0]->cmap->N != X->cmap->N) {
470 PetscCall(MatDestroy(&pch2opus->wnsmat[0]));
471 PetscCall(MatDestroy(&pch2opus->wnsmat[1]));
473 if (!pch2opus->wnsmat[0]) {
474 PetscCall(MatDuplicate(X, MAT_SHARE_NONZERO_PATTERN, &pch2opus->wnsmat[0]));
475 PetscCall(MatDuplicate(Y, MAT_SHARE_NONZERO_PATTERN, &pch2opus->wnsmat[1]));
479 PetscCall(MatTransposeMatMult(A, Y, MAT_REUSE_MATRIX, PETSC_CURRENT, &pch2opus->wnsmat[1]));
480 PetscCall(PCApplyTransposeMat_H2OPUS(pc, pch2opus->wnsmat[1], pch2opus->wnsmat[0]));
482 PetscCall(MatAXPY(Y, -1., pch2opus->wnsmat[0], SAME_NONZERO_PATTERN));
485 PetscCall(MatMatMult(A, Y, MAT_REUSE_MATRIX, PETSC_CURRENT, &pch2opus->wnsmat[0]));
486 PetscCall(PCApplyMat_H2OPUS(pc, pch2opus->wnsmat[0], pch2opus->wnsmat[1]));
488 PetscCall(MatAXPY(Y, -1., pch2opus->wnsmat[1], SAME_NONZERO_PATTERN));
493 static PetscErrorCode MatMatMultNumeric_NS(Mat M, Mat X, Mat Y, PetscCtx ctx)
496 PetscCall(MatMatMultKernel_NS(M, X, Y, PETSC_FALSE));
502 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
503 Mat A = pc->useAmat ? pc->mat : pc->pmat;
506 if (!pch2opus->S) {
507 PetscInt M, N, m, n;
509 PetscCall(MatGetSize(A, &M, &N));
510 PetscCall(MatGetLocalSize(A, &m, &n));
511 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A), m, n, M, N, pc, &pch2opus->S));
512 PetscCall(MatSetBlockSizesFromMats(pch2opus->S, A, A));
514 PetscCall(MatShellSetVecType(pch2opus->S, VECCUDA));
517 if (pch2opus->hyperorder >= 2) {
518 PetscCall(MatShellSetOperation(pch2opus->S, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Hyper));
519 PetscCall(MatShellSetOperation(pch2opus->S, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Hyper));
520 PetscCall(MatShellSetMatProductOperation(pch2opus->S, MATPRODUCT_AB, NULL, MatMatMultNumeric_Hyper, NULL, MATDENSE, MATDENSE));
521 PetscCall(MatShellSetMatProductOperation(pch2opus->S, MATPRODUCT_AB, NULL, MatMatMultNumeric_Hyper, NULL, MATDENSECUDA, MATDENSECUDA));
523 PetscCall(MatShellSetOperation(pch2opus->S, MATOP_MULT, (PetscErrorCodeFn *)MatMult_NS));
524 PetscCall(MatShellSetOperation(pch2opus->S, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_NS));
525 PetscCall(MatShellSetMatProductOperation(pch2opus->S, MATPRODUCT_AB, NULL, MatMatMultNumeric_NS, NULL, MATDENSE, MATDENSE));
526 PetscCall(MatShellSetMatProductOperation(pch2opus->S, MATPRODUCT_AB, NULL, MatMatMultNumeric_NS, NULL, MATDENSECUDA, MATDENSECUDA));
528 PetscCall(MatPropagateSymmetryOptions(A, pch2opus->S));
529 PetscCall(MatBindToCPU(pch2opus->S, pch2opus->boundtocpu));
531 PetscCall(MatSetOption(pch2opus->S, MAT_SYMMETRIC, PETSC_TRUE));
537 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
538 Mat A = pc->useAmat ? pc->mat : pc->pmat;
539 NormType norm = pch2opus->normtype;
544 if (!pch2opus->T) {
545 PetscInt M, N, m, n;
549 PetscCall(MatGetSize(A, &M, &N));
550 PetscCall(MatGetLocalSize(A, &m, &n));
551 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc->pmat), m, n, M, N, pc, &pch2opus->T));
552 PetscCall(MatSetBlockSizesFromMats(pch2opus->T, A, A));
553 PetscCall(MatShellSetOperation(pch2opus->T, MATOP_MULT, (PetscErrorCodeFn *)MatMult_MAmI));
554 PetscCall(MatShellSetOperation(pch2opus->T, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_MAmI));
555 PetscCall(MatShellSetOperation(pch2opus->T, MATOP_NORM, (PetscErrorCodeFn *)MatNorm_H2OPUS));
557 PetscCall(MatShellSetVecType(pch2opus->T, VECCUDA));
559 PetscCall(MatSetOptionsPrefix(pch2opus->T, prefix));
560 PetscCall(MatAppendOptionsPrefix(pch2opus->T, "pc_h2opus_est_"));
562 PetscCall(MatDestroy(&pch2opus->A));
566 pch2opus->A = A;
571 if (!pch2opus->sdim) PetscCall(PCH2OpusInferCoordinates_Private(pc));
572 PetscCall(MatCreateH2OpusFromMat(A, pch2opus->sdim, pch2opus->coords, PETSC_FALSE, pch2opus->eta, pch2opus->leafsize, pch2opus->max_rank, pch2opus->bs, pch2opus->mrtol, &pch2opus->A));
574 PetscCall(MatSetOption(pch2opus->A, MAT_SYMMETRIC, PETSC_TRUE));
576 PetscCall(MatSetOptionsPrefix(pch2opus->A, prefix));
577 PetscCall(MatAppendOptionsPrefix(pch2opus->A, "pc_h2opus_init_"));
578 PetscCall(MatSetFromOptions(pch2opus->A));
579 PetscCall(MatAssemblyBegin(pch2opus->A, MAT_FINAL_ASSEMBLY));
580 PetscCall(MatAssemblyEnd(pch2opus->A, MAT_FINAL_ASSEMBLY));
582 PetscCall(MatSetOption(pch2opus->A, MAT_SYMMETRIC, PETSC_TRUE));
585 PetscCall(MatBindToCPU(pch2opus->A, pch2opus->forcecpu));
588 pch2opus->boundtocpu = pch2opus->forcecpu ? PETSC_TRUE : pch2opus->A->boundtocpu;
590 PetscCall(MatBindToCPU(pch2opus->T, pch2opus->boundtocpu));
591 if (pch2opus->M) { /* see if we can reuse M as initial guess */
594 PetscCall(MatBindToCPU(pch2opus->M, pch2opus->boundtocpu));
595 PetscCall(MatNorm(pch2opus->T, norm, &reuse));
596 if (reuse >= 1.0) PetscCall(MatDestroy(&pch2opus->M));
598 if (!pch2opus->M) {
600 PetscCall(MatDuplicate(pch2opus->A, MAT_COPY_VALUES, &pch2opus->M));
602 PetscCall(MatSetOptionsPrefix(pch2opus->M, prefix));
603 PetscCall(MatAppendOptionsPrefix(pch2opus->M, "pc_h2opus_inv_"));
604 PetscCall(MatSetFromOptions(pch2opus->M));
606 PetscCall(MatScale(pch2opus->M, pch2opus->s0));
608 /* A and M have the same h2 matrix nonzero structure, save on reordering routines */
609 PetscCall(MatH2OpusSetNativeMult(pch2opus->A, PETSC_TRUE));
610 PetscCall(MatH2OpusSetNativeMult(pch2opus->M, PETSC_TRUE));
611 if (norm == NORM_1 || norm == NORM_2 || norm == NORM_INFINITY) PetscCall(MatNorm(pch2opus->T, norm, &initerr));
612 if (PetscIsInfOrNanReal(initerr)) pc->failedreason = PC_SETUP_ERROR;
614 if (pch2opus->monitor) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)pc), "%" PetscInt_FMT ": ||M*A - I|| NORM%s abs %g rel %g\n", 0, NormTypes[norm], (double)err, (double)(err / initerr)));
615 if (initerr > pch2opus->atol && !pc->failedreason) {
619 for (i = 0; i < pch2opus->maxits; i++) {
620 Mat M;
623 PetscCall(MatDuplicate(pch2opus->M, MAT_SHARE_NONZERO_PATTERN, &M));
624 PetscCall(MatGetOptionsPrefix(pch2opus->M, &prefix));
625 PetscCall(MatSetOptionsPrefix(M, prefix));
626 PetscCall(MatH2OpusSetSamplingMat(M, pch2opus->S, PETSC_DECIDE, PETSC_DECIDE));
627 PetscCall(MatSetFromOptions(M));
628 PetscCall(MatH2OpusSetNativeMult(M, PETSC_TRUE));
629 PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
630 PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
631 PetscCall(MatDestroy(&pch2opus->M));
632 pch2opus->M = M;
633 if (norm == NORM_1 || norm == NORM_2 || norm == NORM_INFINITY) PetscCall(MatNorm(pch2opus->T, norm, &err));
634 if (pch2opus->monitor) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)pc), "%" PetscInt_FMT ": ||M*A - I|| NORM%s abs %g rel %g\n", i + 1, NormTypes[norm], (double)err, (double)(err / initerr)));
635 if (PetscIsInfOrNanReal(err)) pc->failedreason = PC_SETUP_ERROR;
636 if (err < pch2opus->atol || err < pch2opus->rtol * initerr || pc->failedreason) break;
640 PetscCall(MatH2OpusSetNativeMult(pch2opus->A, PETSC_FALSE));
641 PetscCall(MatH2OpusSetNativeMult(pch2opus->M, PETSC_FALSE));
642 PetscCall(VecDestroy(&pch2opus->wns[0]));
643 PetscCall(VecDestroy(&pch2opus->wns[1]));
644 PetscCall(VecDestroy(&pch2opus->wns[2]));
645 PetscCall(VecDestroy(&pch2opus->wns[3]));
646 PetscCall(MatDestroy(&pch2opus->wnsmat[0]));
647 PetscCall(MatDestroy(&pch2opus->wnsmat[1]));
648 PetscCall(MatDestroy(&pch2opus->wnsmat[2]));
649 PetscCall(MatDestroy(&pch2opus->wnsmat[3]));
655 PC_H2OPUS *pch2opus = (PC_H2OPUS *)pc->data;
661 if (pch2opus->A && pch2opus->A != pc->mat && pch2opus->A != pc->pmat) {
665 PetscCall(MatView(pch2opus->A, viewer));
669 if (pch2opus->M) {
673 PetscCall(MatView(pch2opus->M, viewer));
677 PetscCall(PetscViewerASCIIPrintf(viewer, "Initial scaling: %g\n", pch2opus->s0));
683 PCH2OPUS = "h2opus" - A preconditioner type for, `MATH2OPUS`, hierarchical matrices using the H2Opus package.
686 + -pc_type h2opus - pc type to "h2opus" during a call to `PCSetFromOptions()`
687 . -pc_h2opus_maxits - maximum number of iterations for Newton-Schultz
688 . -pc_h2opus_monitor - monitor Newton-Schultz convergence
689 . -pc_h2opus_atol - absolute tolerance
690 . -pc_h2opus_rtol - relative tolerance
691 . -pc_h2opus_norm_type - normtype
692 . -pc_h2opus_hyperorder - Hyper power order of sampling
693 . -pc_h2opus_leafsize - leaf size when constructed from kernel
694 . -pc_h2opus_eta - admissibility condition tolerance
695 . -pc_h2opus_maxrank - maximum rank when constructed from matvecs
696 . -pc_h2opus_samples - number of samples to be taken concurrently when constructing from matvecs
697 . -pc_h2opus_mrtol - relative tolerance for construction from sampling
698 - -pc_h2opus_forcecpu - force construction of preconditioner on CPU
703 M*/
711 pc->data = (void *)pch2opus;
713 pch2opus->atol = 1.e-2;
714 pch2opus->rtol = 1.e-6;
715 pch2opus->maxits = 50;
716 pch2opus->hyperorder = 1; /* defaults to basic Newton-Schultz */
717 pch2opus->normtype = NORM_2;
720 pch2opus->eta = PETSC_DECIDE;
721 pch2opus->leafsize = PETSC_DECIDE;
722 pch2opus->max_rank = PETSC_DECIDE;
723 pch2opus->bs = PETSC_DECIDE;
724 pch2opus->mrtol = PETSC_DECIDE;
725 pch2opus->boundtocpu = PetscDefined(H2OPUS_USE_GPU) ? PETSC_FALSE : PETSC_TRUE;
726 pc->ops->destroy = PCDestroy_H2OPUS;
727 pc->ops->setup = PCSetUp_H2OPUS;
728 pc->ops->apply = PCApply_H2OPUS;
729 pc->ops->matapply = PCApplyMat_H2OPUS;
730 pc->ops->applytranspose = PCApplyTranspose_H2OPUS;
731 pc->ops->reset = PCReset_H2OPUS;
732 pc->ops->setfromoptions = PCSetFromOptions_H2OPUS;
733 pc->ops->view = PCView_H2OPUS;