Lines Matching +full:test +full:- +full:experimental
68 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, nr, nr - i, NULL, B)); in MatDenseOrthogonalRangeOrComplement()
70 PetscCall(PetscArraycpy(data, U + nr * i, (nr - i) * nr)); in MatDenseOrthogonalRangeOrComplement()
155 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCNedelecSupport()
156 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCNedelecSupport()
180 order = pcbddc->nedorder; in PCBDDCNedelecSupport()
181 conforming = pcbddc->conforming; in PCBDDCNedelecSupport()
182 field = pcbddc->nedfield; in PCBDDCNedelecSupport()
183 global = pcbddc->nedglobal; in PCBDDCNedelecSupport()
188 …PetscOptionsBegin(PetscObjectComm((PetscObject)pc), ((PetscObject)pc)->prefix, "BDDC Nedelec optio… in PCBDDCNedelecSupport()
189 …PetscCall(PetscOptionsBool("-pc_bddc_nedelec_field_primal", "All edge dofs set as primals: Toselli… in PCBDDCNedelecSupport()
191 …PetscCall(PetscOptionsInt("-pc_bddc_nedelec_order", "Test variable order code (to be removed)", NU… in PCBDDCNedelecSupport()
192 …PetscCall(PetscOptionsBool("-pc_bddc_nedelec_print", "Print debug info", NULL, print, &print, NULL… in PCBDDCNedelecSupport()
196 PetscCall(MatISGetLocalToGlobalMapping(pc->pmat, &al2g, NULL)); in PCBDDCNedelecSupport()
199 PetscCall(VecGetArrayRead(matis->counter, (const PetscScalar **)&vals)); in PCBDDCNedelecSupport()
207 PetscCall(VecRestoreArrayRead(matis->counter, (const PetscScalar **)&vals)); in PCBDDCNedelecSupport()
212 …->n_ISForDofsLocal || field < pcbddc->n_ISForDofsLocal, comm, PETSC_ERR_USER, "Invalid field for N… in PCBDDCNedelecSupport()
213 if (pcbddc->n_ISForDofsLocal && field >= 0) { in PCBDDCNedelecSupport()
214 PetscCall(PetscObjectReference((PetscObject)pcbddc->ISForDofsLocal[field])); in PCBDDCNedelecSupport()
215 nedfieldlocal = pcbddc->ISForDofsLocal[field]; in PCBDDCNedelecSupport()
217 } else if (!pcbddc->n_ISForDofsLocal && field != PETSC_DECIDE) { in PCBDDCNedelecSupport()
224 PetscCall(PetscArrayzero(matis->sf_leafdata, n)); in PCBDDCNedelecSupport()
225 PetscCall(PetscArrayzero(matis->sf_rootdata, pc->pmat->rmap->n)); in PCBDDCNedelecSupport()
226 PetscCall(MatGetOwnershipRange(pcbddc->discretegradient, &rst, &ren)); in PCBDDCNedelecSupport()
230 PetscCall(MatGetRow(pcbddc->discretegradient, i, &nc, NULL, NULL)); in PCBDDCNedelecSupport()
231 if (nc > 1) matis->sf_rootdata[i - rst] = 1; in PCBDDCNedelecSupport()
232 PetscCall(MatRestoreRow(pcbddc->discretegradient, i, &nc, NULL, NULL)); in PCBDDCNedelecSupport()
234 …PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLA… in PCBDDCNedelecSupport()
235 …PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE… in PCBDDCNedelecSupport()
238 if (matis->sf_leafdata[i]) idx[ne++] = i; in PCBDDCNedelecSupport()
245 …PetscCheck(order || conforming, comm, PETSC_ERR_SUP, "Variable order and non-conforming spaces are… in PCBDDCNedelecSupport()
246 …PetscCheck(!pcbddc->user_ChangeOfBasisMatrix, comm, PETSC_ERR_SUP, "Cannot generate Nedelec suppor… in PCBDDCNedelecSupport()
255 PetscCall(VecGetArrayRead(matis->counter, (const PetscScalar **)&vals)); in PCBDDCNedelecSupport()
267 PetscCall(VecRestoreArrayRead(matis->counter, (const PetscScalar **)&vals)); in PCBDDCNedelecSupport()
306 PetscCall(PetscArrayzero(matis->sf_leafdata, n)); in PCBDDCNedelecSupport()
307 PetscCall(PetscArrayzero(matis->sf_rootdata, pc->pmat->rmap->n)); in PCBDDCNedelecSupport()
310 for (i = 0; i < ne; i++) matis->sf_leafdata[idxs[i]] = 1; in PCBDDCNedelecSupport()
313 for (i = 0; i < ne; i++) matis->sf_leafdata[i] = 1; in PCBDDCNedelecSupport()
315 …PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_SUM)… in PCBDDCNedelecSupport()
316 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_SUM)); in PCBDDCNedelecSupport()
318 …/* There's no way to detect all possible corner candidates in a element-by-element case in a pure … in PCBDDCNedelecSupport()
320 …if (matis->allow_repeated) PetscCall(PetscObjectQuery((PetscObject)pcbddc->discretegradient, "_ele… in PCBDDCNedelecSupport()
323 …PetscCall(MatViewFromOptions(pcbddc->discretegradient, (PetscObject)pc, "-pc_bddc_discrete_gradien… in PCBDDCNedelecSupport()
324 PetscCall(MatDuplicate(pcbddc->discretegradient, MAT_COPY_VALUES, &G)); in PCBDDCNedelecSupport()
330 for (i = 0, cum = 0; i < pc->pmat->rmap->n; i++) { in PCBDDCNedelecSupport()
331 if (matis->sf_rootdata[i] < 2) matis->sf_rootdata[cum++] = i + rst; in PCBDDCNedelecSupport()
334 PetscCall(MatZeroRows(G, cum, matis->sf_rootdata, 0., NULL, NULL)); in PCBDDCNedelecSupport()
339 …PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLA… in PCBDDCNedelecSupport()
340 …PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE… in PCBDDCNedelecSupport()
343 if (matis->sf_leafdata[idxs[i]] == 1) tbz[cum++] = i; in PCBDDCNedelecSupport()
360 if (matis->allow_repeated) { /* multi-element support */ in PCBDDCNedelecSupport()
367 PetscCall(PetscCalloc1(pcbddc->n_local_subs * pcbddc->n_local_subs, &lGn)); in PCBDDCNedelecSupport()
368 PetscCall(PetscMalloc1(pcbddc->n_local_subs, &is_rows)); in PCBDDCNedelecSupport()
369 PetscCall(PetscMalloc1(pcbddc->n_local_subs, &tcols)); in PCBDDCNedelecSupport()
370 for (PetscInt i = 0; i < pcbddc->n_local_subs; i++) { in PCBDDCNedelecSupport()
372 … PetscCall(ISGlobalToLocalMappingApplyIS(fl2g, IS_GTOLM_MASK, pcbddc->local_subs[i], &is_rows[i])); in PCBDDCNedelecSupport()
374 PetscCall(PetscObjectReference((PetscObject)pcbddc->local_subs[i])); in PCBDDCNedelecSupport()
375 is_rows[i] = pcbddc->local_subs[i]; in PCBDDCNedelecSupport()
377 …PetscCall(MatCreateSubMatrix(lG, is_rows[i], NULL, MAT_INITIAL_MATRIX, &lGn[i * (1 + pcbddc->n_loc… in PCBDDCNedelecSupport()
378 PetscCall(MatSeqAIJCompactOutExtraColumns_SeqAIJ(lGn[i * (1 + pcbddc->n_local_subs)], &mapn)); in PCBDDCNedelecSupport()
388 PetscCall(MatSetSizes(B, lGis->rmap->n, lGis->cmap->n, lGis->rmap->N, lGis->cmap->N)); in PCBDDCNedelecSupport()
391 PetscCall(ISConcatenate(PETSC_COMM_SELF, pcbddc->n_local_subs, tcols, &tmap)); in PCBDDCNedelecSupport()
392 PetscCall(ISLocalToGlobalMappingApplyIS(lGis->cmap->mapping, tmap, &nmap)); in PCBDDCNedelecSupport()
401 PetscCall(MatSetLocalToGlobalMapping(B, lGis->rmap->mapping, mapn)); in PCBDDCNedelecSupport()
403 …PetscCall(MatCreateNest(PETSC_COMM_SELF, pcbddc->n_local_subs, is_rows, pcbddc->n_local_subs, NULL… in PCBDDCNedelecSupport()
404 for (PetscInt i = 0; i < pcbddc->n_local_subs; i++) { in PCBDDCNedelecSupport()
405 PetscCall(MatDestroy(&lGn[i * (1 + pcbddc->n_local_subs)])); in PCBDDCNedelecSupport()
419 lGis->assembled = PETSC_TRUE; in PCBDDCNedelecSupport()
421 PetscCall(MatViewFromOptions(lGis, (PetscObject)pc, "-pc_bddc_nedelec_init_G_view")); in PCBDDCNedelecSupport()
430 PetscCall(PetscSFSetGraphLayout(sfv, lGis->cmap, nv, NULL, PETSC_OWN_POINTER, idxs)); in PCBDDCNedelecSupport()
436 Mat_IS *tGis = (Mat_IS *)lGis->data; in PCBDDCNedelecSupport()
439 PetscCall(MatCreateVecs(tGis->A, &local, NULL)); in PCBDDCNedelecSupport()
440 PetscCall(PCBDDCGlobalToLocal(tGis->cctx, global, local, elements_corners, &tmp)); in PCBDDCNedelecSupport()
460 /* Analyze the edge-nodes connections (duplicate lG) */ in PCBDDCNedelecSupport()
469 if (pcbddc->DirichletBoundariesLocal) { in PCBDDCNedelecSupport()
473 …PetscCall(ISGlobalToLocalMappingApplyIS(fl2g, IS_GTOLM_MASK, pcbddc->DirichletBoundariesLocal, &is… in PCBDDCNedelecSupport()
475 is = pcbddc->DirichletBoundariesLocal; in PCBDDCNedelecSupport()
488 if (pcbddc->NeumannBoundariesLocal) { in PCBDDCNedelecSupport()
492 …PetscCall(ISGlobalToLocalMappingApplyIS(fl2g, IS_GTOLM_MASK, pcbddc->NeumannBoundariesLocal, &is)); in PCBDDCNedelecSupport()
494 is = pcbddc->NeumannBoundariesLocal; in PCBDDCNedelecSupport()
531 …if (ii[i + 1] - ii[i] != order + 1) { /* every row of G on the coarse edge should list order+1 nod… in PCBDDCNedelecSupport()
538 - at most 2 endpoints in PCBDDCNedelecSupport()
539 - order-1 interior nodal dofs in PCBDDCNedelecSupport()
540 - no undefined nodal dofs (nconn < order) in PCBDDCNedelecSupport()
545 PetscInt nconn = iit[v + 1] - iit[v]; in PCBDDCNedelecSupport()
547 if (!PetscBTLookup(btee, jjt[k])) nconn--; in PCBDDCNedelecSupport()
552 if (undef || ends > 2 || ints != order - 1) { in PCBDDCNedelecSupport()
559 /* We assume the order on the element edge is ii[i+1]-ii[i]-1 */ in PCBDDCNedelecSupport()
561 PetscScalar val = 1. / (ii[i + 1] - ii[i] - 1); in PCBDDCNedelecSupport()
584 …if (matis->allow_repeated) { /* assign a uniq global id to edge local subsets and communicate it w… in PCBDDCNedelecSupport()
587 PetscInt cum_subs = 0, n_subs = pcbddc->n_local_subs, bs, emnr, emnl, vmnr, vmnl; in PCBDDCNedelecSupport()
599 …PetscCheck(emnr == ne, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of roots in edge multi-lea… in PCBDDCNedelecSupport()
600 …PetscCheck(emnl == j, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of leaves in edge multi-lea… in PCBDDCNedelecSupport()
604 …nv, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of roots in nodal multi-leaves SF %" PetscInt… in PCBDDCNedelecSupport()
605 …j, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of leaves in nodal multi-leaves SF %" PetscInt… in PCBDDCNedelecSupport()
621 PetscCall(ISGetLocalSize(pcbddc->local_subs[i], &ns)); in PCBDDCNedelecSupport()
622 PetscCall(ISGetIndices(pcbddc->local_subs[i], &idxs)); in PCBDDCNedelecSupport()
629 PetscCall(ISRestoreIndices(pcbddc->local_subs[i], &idxs)); in PCBDDCNedelecSupport()
666 PetscInt ord = order, test = ii[i + 1] - ii[i], vc = vcount[i]; in PCBDDCNedelecSupport() local
672 test = PetscFloorReal(vorder + 10. * PETSC_SQRT_MACHINE_EPSILON); in PCBDDCNedelecSupport()
673 …er - test <= PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected value for vo… in PCBDDCNedelecSupport()
696 if (elements_corners) test = 0; in PCBDDCNedelecSupport()
697 if (!sneighs || test >= 3 * ord || bdir) { /* splitpoints */ in PCBDDCNedelecSupport()
698 …OINT %" PetscInt_FMT " (%s %s %s)\n", i, PetscBools[!sneighs], PetscBools[test >= 3 * ord], PetscB… in PCBDDCNedelecSupport()
700 } else if (test == ord) { in PCBDDCNedelecSupport()
701 if (order == 1 || (!order && ii[i + 1] - ii[i] == 1)) { in PCBDDCNedelecSupport()
712 /* a candidate is valid if it is connected to another candidate via a non-primal edge dof */ in PCBDDCNedelecSupport()
755 if (matis->allow_repeated) { in PCBDDCNedelecSupport()
842 …if (!elements_corners) { /* if present, we assume we are in the element-by-element case and the CS… in PCBDDCNedelecSupport()
858 if (pcbddc->mat_graph->nvtxs_csr) { /* the user has passed in a CSR graph */ in PCBDDCNedelecSupport()
859 …->mat_graph->nvtxs_csr == n, PETSC_COMM_SELF, PETSC_ERR_USER, "Invalid size of CSR graph %" PetscI… in PCBDDCNedelecSupport()
860 iiu = pcbddc->mat_graph->xadj; in PCBDDCNedelecSupport()
861 jju = pcbddc->mat_graph->adjncy; in PCBDDCNedelecSupport()
862 } else if (pcbddc->use_local_adj) { in PCBDDCNedelecSupport()
864 …PetscCall(MatGetRowIJ(matis->A, 0, PETSC_TRUE, PETSC_FALSE, &i, (const PetscInt **)&iiu, (const Pe… in PCBDDCNedelecSupport()
871 jju[i] = -1; in PCBDDCNedelecSupport()
877 for (i = 0; i < n; i++) iia[i + 1] = iiu[i + 1] - iiu[i]; in PCBDDCNedelecSupport()
884 iia[idxs[i] + 1] = ii[i + 1] - ii[i]; in PCBDDCNedelecSupport()
894 for (j = 0; j < iiu[i + 1] - iiu[i]; j++) jja[iia[i] + j] = jju[iiu[i] + j]; in PCBDDCNedelecSupport()
901 for (j = 0; j < ii[i + 1] - ii[i]; j++) jja[iia[e] + j] = jj[ii[i] + j]; in PCBDDCNedelecSupport()
906 …if (rest) PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_TRUE, PETSC_FALSE, &i, (const PetscInt **)&… in PCBDDCNedelecSupport()
918 pcbddc->mat_graph->twodim = PETSC_FALSE; in PCBDDCNedelecSupport()
921 …PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, &nee, &alleedges, &allprimals)… in PCBDDCNedelecSupport()
1047 if (cum != size - 1 || found != 2) { in PCBDDCNedelecSupport()
1138 …PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, NULL, NULL, &nee, &alleedges, &allprim… in PCBDDCNedelecSupport()
1144 pcbddc->mat_graph->twodim = PETSC_FALSE; in PCBDDCNedelecSupport()
1145 …PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, &nee, &alleedges, &allprimals)… in PCBDDCNedelecSupport()
1191 if (cum != size - 1) { in PCBDDCNedelecSupport()
1205 if (print) PetscCall(PCBDDCGraphASCIIView(pcbddc->mat_graph, 5, PETSC_VIEWER_STDOUT_SELF)); in PCBDDCNedelecSupport()
1217 cedges[i] = corners[i * 2] = corners[i * 2 + 1] = -1; in PCBDDCNedelecSupport()
1250 cedges[i] = idxs[size - 1]; in PCBDDCNedelecSupport()
1278 …gh the %" PetscInt_FMT " nodal dof at edge dof %" PetscInt_FMT, marks[jj[j]] - 1, eemax, i, jj[j]); in PCBDDCNedelecSupport()
1303 mark--; in PCBDDCNedelecSupport()
1305 size = ii[i + 1] - ii[i]; in PCBDDCNedelecSupport()
1331 PetscCall(MatSetLayouts(T, pc->mat->rmap, pc->mat->cmap)); in PCBDDCNedelecSupport()
1360 for (i = T->rmap->rstart; i < T->rmap->rend; i++) in PCBDDCNedelecSupport()
1361 … if (PetscAbsScalar(wa[i - T->rmap->rstart])) PetscCall(MatSetValue(T, i, i, 1.0, INSERT_VALUES)); in PCBDDCNedelecSupport()
1367 PetscCall(MatDestroy(&pcbddc->nedcG)); in PCBDDCNedelecSupport()
1368 PetscCall(ISDestroy(&pcbddc->nedclocal)); in PCBDDCNedelecSupport()
1369 if (pcbddc->current_level < pcbddc->max_levels) { in PCBDDCNedelecSupport()
1376 PetscCall(ISLocalToGlobalMappingApplyIS(fl2g, wis, &pcbddc->nedclocal)); in PCBDDCNedelecSupport()
1379 pcbddc->nedclocal = wis; in PCBDDCNedelecSupport()
1396 PetscCall(MatCreate(comm, &pcbddc->nedcG)); in PCBDDCNedelecSupport()
1397 PetscCall(MatSetSizes(pcbddc->nedcG, PETSC_DECIDE, PETSC_DECIDE, cne, cnv)); in PCBDDCNedelecSupport()
1398 PetscCall(MatSetType(pcbddc->nedcG, MATAIJ)); in PCBDDCNedelecSupport()
1399 PetscCall(MatSeqAIJSetPreallocation(pcbddc->nedcG, 2, NULL)); in PCBDDCNedelecSupport()
1400 PetscCall(MatMPIAIJSetPreallocation(pcbddc->nedcG, 2, NULL, 2, NULL)); in PCBDDCNedelecSupport()
1401 PetscCall(MatSetLocalToGlobalMapping(pcbddc->nedcG, cel2g, cvl2g)); in PCBDDCNedelecSupport()
1411 PetscCall(MatGetNullSpace(pcbddc->discretegradient, &nnsp)); in PCBDDCNedelecSupport()
1419 PetscCall(MatCreateVecs(pc->pmat, &E, NULL)); in PCBDDCNedelecSupport()
1420 PetscCall(MatCreateVecs(pcbddc->discretegradient, &V, NULL)); in PCBDDCNedelecSupport()
1428 cum += j - 1; in PCBDDCNedelecSupport()
1430 PetscCall(PetscMalloc1(PetscMax(cum, pc->pmat->rmap->n), &eedgesidxs)); in PCBDDCNedelecSupport()
1438 PetscCall(PetscArraycpy(eedgesidxs + cum, idxs, j - 1)); /* last on the edge is primal */ in PCBDDCNedelecSupport()
1440 cum += j - 1; in PCBDDCNedelecSupport()
1450 /* identify dofs we must zero if importing user-defined near nullspace from pmat */ in PCBDDCNedelecSupport()
1456 for (i = 0, cum = 0; i < pc->pmat->rmap->n; i++) in PCBDDCNedelecSupport()
1457 if (evals[i] == 0.0) eedgesidxs[cum++] = i + pc->pmat->rmap->rstart; in PCBDDCNedelecSupport()
1471 lev = pcbddc->current_level; in PCBDDCNedelecSupport()
1480 …if (pcbddc->nedcG) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 2, corners + 2 * i, PETSC_USE_POINTE… in PCBDDCNedelecSupport()
1499 …PetscCheck(ncc == 1 || !pcbddc->nedcG, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot generate the coarse… in PCBDDCNedelecSupport()
1505 if (pcbddc->nedcG) { in PCBDDCNedelecSupport()
1510 PetscCall(MatSetValuesLocal(pcbddc->nedcG, 1, &i, 2, cols, cvals, INSERT_VALUES)); in PCBDDCNedelecSupport()
1531 /* import nearnullspace from preconditioning matrix if user-defined */ in PCBDDCNedelecSupport()
1532 PetscCall(MatGetNearNullSpace(pc->pmat, &onnsp)); in PCBDDCNedelecSupport()
1536 … PetscCall(PetscStrcmp("_internal_BDDC_nedelec_nnsp", ((PetscObject)onnsp)->name, &isinternal)); in PCBDDCNedelecSupport()
1577 PetscCall(MatSetNearNullSpace(pc->pmat, nnsp)); in PCBDDCNedelecSupport()
1587 if (pcbddc->nedcG) PetscCall(MatAssemblyBegin(pcbddc->nedcG, MAT_FINAL_ASSEMBLY)); in PCBDDCNedelecSupport()
1598 PCBDDCGraph graph = pcbddc->mat_graph; in PCBDDCNedelecSupport()
1599 PetscInt *oqueue = graph->queue; in PCBDDCNedelecSupport()
1600 PetscInt *ocptr = graph->cptr; in PCBDDCNedelecSupport()
1604 if (pcbddc->nedclocal) { in PCBDDCNedelecSupport()
1605 PetscCall(ISGetIndices(pcbddc->nedclocal, (const PetscInt **)&idxs)); in PCBDDCNedelecSupport()
1614 PetscCall(PetscMalloc2(graph->nvtxs + 1, &graph->cptr, ocptr[graph->ncc], &graph->queue)); in PCBDDCNedelecSupport()
1615 graph->cptr[0] = 0; in PCBDDCNedelecSupport()
1616 for (i = 0, ncc = 0; i < graph->ncc; i++) { in PCBDDCNedelecSupport()
1617 PetscInt lc = ocptr[i + 1] - ocptr[i]; in PCBDDCNedelecSupport()
1618 if (cum != nee && oqueue[ocptr[i + 1] - 1] == cedges[cum]) { /* this cc has a primal dof */ in PCBDDCNedelecSupport()
1619 graph->cptr[ncc + 1] = graph->cptr[ncc] + 1; in PCBDDCNedelecSupport()
1620 graph->queue[graph->cptr[ncc]] = cedges[cum]; in PCBDDCNedelecSupport()
1622 lc--; in PCBDDCNedelecSupport()
1626 graph->cptr[ncc + 1] = graph->cptr[ncc] + lc; in PCBDDCNedelecSupport()
1627 for (j = 0; j < lc; j++) graph->queue[graph->cptr[ncc] + j] = oqueue[ocptr[i] + j]; in PCBDDCNedelecSupport()
1630 graph->ncc = ncc; in PCBDDCNedelecSupport()
1631 if (pcbddc->nedclocal) PetscCall(ISRestoreIndices(pcbddc->nedclocal, (const PetscInt **)&idxs)); in PCBDDCNedelecSupport()
1635 …PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, NULL, NULL, &nee, &alleedges, &allprim… in PCBDDCNedelecSupport()
1636 PetscCall(PCBDDCGraphResetCSR(pcbddc->mat_graph)); in PCBDDCNedelecSupport()
1649 PetscCall(MatViewFromOptions(T, (PetscObject)pc, "-pc_bddc_nedelec_change_view")); in PCBDDCNedelecSupport()
1650 if (pcbddc->nedcG) { in PCBDDCNedelecSupport()
1651 PetscCall(MatAssemblyEnd(pcbddc->nedcG, MAT_FINAL_ASSEMBLY)); in PCBDDCNedelecSupport()
1652 …PetscCall(MatViewFromOptions(pcbddc->nedcG, (PetscObject)pc, "-pc_bddc_nedelec_coarse_change_view"… in PCBDDCNedelecSupport()
1663 /* the near-null space of BDDC carries information on quadrature weights,
1664 and these can be collinear -> so cheat with MatNullSpaceCreate
1675 …PetscCheck(last - first >= 2 * nvecs || !has_const, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implement… in PCBDDCNullSpaceCreate()
1680 data[i - first] = 1.; in PCBDDCNullSpaceCreate()
1682 data[2 * i - first] = 1. / PetscSqrtReal(2.); in PCBDDCNullSpaceCreate()
1683 data[2 * i - first + 1] = -1. / PetscSqrtReal(2.); in PCBDDCNullSpaceCreate()
1698 data[i - first] = 0.; in PCBDDCNullSpaceCreate()
1700 data[2 * i - first] = 0.; in PCBDDCNullSpaceCreate()
1701 data[2 * i - first + 1] = 0.; in PCBDDCNullSpaceCreate()
1765 PetscCall(PetscSFGetGraph(graph->interface_subset_sf, &nr, NULL, NULL, NULL)); in PCBDDCComputeNoNetFlux()
1766 PetscCall(PetscSFGetMultiSF(graph->interface_subset_sf, &msf)); in PCBDDCComputeNoNetFlux()
1769 PetscCall(PetscSFComputeDegreeBegin(graph->interface_subset_sf, °ree)); in PCBDDCComputeNoNetFlux()
1770 PetscCall(PetscSFComputeDegreeEnd(graph->interface_subset_sf, °ree)); in PCBDDCComputeNoNetFlux()
1775 PetscCall(PetscSFScatterBegin(graph->interface_subset_sf, MPIU_INT, mmask, mask)); in PCBDDCComputeNoNetFlux()
1776 PetscCall(PetscSFScatterEnd(graph->interface_subset_sf, MPIU_INT, mmask, mask)); in PCBDDCComputeNoNetFlux()
1788 PetscCall(VecViewFromOptions(quad_vec, NULL, "-pc_bddc_quad_vec_view")); in PCBDDCComputeNoNetFlux()
1796 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCAddPrimalVerticesLocalIS()
1800 if (pcbddc->user_primal_vertices_local) { in PCBDDCAddPrimalVerticesLocalIS()
1804 list[1] = pcbddc->user_primal_vertices_local; in PCBDDCAddPrimalVerticesLocalIS()
1808 pcbddc->user_primal_vertices_local = newp; in PCBDDCAddPrimalVerticesLocalIS()
1828 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCComputeLocalTopologyInfo()
1829 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCComputeLocalTopologyInfo()
1833 …PetscOptionsBegin(PetscObjectComm((PetscObject)pc), ((PetscObject)pc)->prefix, "BDDC topology opti… in PCBDDCComputeLocalTopologyInfo()
1834 …PetscCall(PetscOptionsBool("-pc_bddc_monolithic", "Discard any information on dofs splitting", NUL… in PCBDDCComputeLocalTopologyInfo()
1837 PetscCall(MatCreateVecs(pc->pmat, &global, NULL)); in PCBDDCComputeLocalTopologyInfo()
1838 PetscCall(MatCreateVecs(matis->A, &local, NULL)); in PCBDDCComputeLocalTopologyInfo()
1842 if (pcbddc->vertex_size == 1) PetscCall(MatGetBlockSize(pc->pmat, &pcbddc->vertex_size)); in PCBDDCComputeLocalTopologyInfo()
1846 if (pcbddc->user_provided_isfordofs) { in PCBDDCComputeLocalTopologyInfo()
1847 if (pcbddc->n_ISForDofs) { in PCBDDCComputeLocalTopologyInfo()
1850 PetscCall(PetscMalloc1(pcbddc->n_ISForDofs, &pcbddc->ISForDofsLocal)); in PCBDDCComputeLocalTopologyInfo()
1851 for (i = 0; i < pcbddc->n_ISForDofs; i++) { in PCBDDCComputeLocalTopologyInfo()
1854 …PetscCall(PCBDDCGlobalToLocal(matis->rctx, global, local, pcbddc->ISForDofs[i], &pcbddc->ISForDofs… in PCBDDCComputeLocalTopologyInfo()
1855 PetscCall(ISGetBlockSize(pcbddc->ISForDofs[i], &bs)); in PCBDDCComputeLocalTopologyInfo()
1856 PetscCall(ISSetBlockSize(pcbddc->ISForDofsLocal[i], bs)); in PCBDDCComputeLocalTopologyInfo()
1857 PetscCall(ISDestroy(&pcbddc->ISForDofs[i])); in PCBDDCComputeLocalTopologyInfo()
1859 pcbddc->n_ISForDofsLocal = pcbddc->n_ISForDofs; in PCBDDCComputeLocalTopologyInfo()
1860 pcbddc->n_ISForDofs = 0; in PCBDDCComputeLocalTopologyInfo()
1861 PetscCall(PetscFree(pcbddc->ISForDofs)); in PCBDDCComputeLocalTopologyInfo()
1864 if (!pcbddc->n_ISForDofsLocal) { /* field split not present */ in PCBDDCComputeLocalTopologyInfo()
1867 PetscCall(MatGetDM(pc->pmat, &dm)); in PCBDDCComputeLocalTopologyInfo()
1874 PetscCall(PetscMalloc1(nf, &pcbddc->ISForDofsLocal)); in PCBDDCComputeLocalTopologyInfo()
1878 … PetscCall(PCBDDCGlobalToLocal(matis->rctx, global, local, fields[i], &pcbddc->ISForDofsLocal[i])); in PCBDDCComputeLocalTopologyInfo()
1880 PetscCall(ISSetBlockSize(pcbddc->ISForDofsLocal[i], bs)); in PCBDDCComputeLocalTopologyInfo()
1884 pcbddc->n_ISForDofsLocal = nf; in PCBDDCComputeLocalTopologyInfo()
1888 … PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "_convert_nest_lfields", (PetscObject *)&c)); in PCBDDCComputeLocalTopologyInfo()
1892 PetscCall(PCBDDCSetDofsSplittingLocal(pc, lf->nr, lf->rf)); in PCBDDCComputeLocalTopologyInfo()
1894 PetscInt i, n = matis->A->rmap->n; in PCBDDCComputeLocalTopologyInfo()
1895 PetscCall(MatGetBlockSize(pc->pmat, &i)); in PCBDDCComputeLocalTopologyInfo()
1897 pcbddc->n_ISForDofsLocal = i; in PCBDDCComputeLocalTopologyInfo()
1898 PetscCall(PetscMalloc1(pcbddc->n_ISForDofsLocal, &pcbddc->ISForDofsLocal)); in PCBDDCComputeLocalTopologyInfo()
1899 …c->n_ISForDofsLocal; i++) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), n / pcbddc->n… in PCBDDCComputeLocalTopologyInfo()
1905 …for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(PCBDDCConsistencyCheckIS(pc, MPI_LAND, &p… in PCBDDCComputeLocalTopologyInfo()
1910 if (!pcbddc->DirichletBoundariesLocal && pcbddc->DirichletBoundaries) { in PCBDDCComputeLocalTopologyInfo()
1911 …PetscCall(PCBDDCGlobalToLocal(matis->rctx, global, local, pcbddc->DirichletBoundaries, &pcbddc->Di… in PCBDDCComputeLocalTopologyInfo()
1912 } else if (pcbddc->DirichletBoundariesLocal) { in PCBDDCComputeLocalTopologyInfo()
1913 PetscCall(PCBDDCConsistencyCheckIS(pc, MPI_LAND, &pcbddc->DirichletBoundariesLocal)); in PCBDDCComputeLocalTopologyInfo()
1915 if (!pcbddc->NeumannBoundariesLocal && pcbddc->NeumannBoundaries) { in PCBDDCComputeLocalTopologyInfo()
1916 …PetscCall(PCBDDCGlobalToLocal(matis->rctx, global, local, pcbddc->NeumannBoundaries, &pcbddc->Neum… in PCBDDCComputeLocalTopologyInfo()
1917 } else if (pcbddc->NeumannBoundariesLocal) { in PCBDDCComputeLocalTopologyInfo()
1918 PetscCall(PCBDDCConsistencyCheckIS(pc, MPI_LOR, &pcbddc->NeumannBoundariesLocal)); in PCBDDCComputeLocalTopologyInfo()
1920 …->user_primal_vertices_local && pcbddc->user_primal_vertices) PetscCall(PCBDDCGlobalToLocal(matis-… in PCBDDCComputeLocalTopologyInfo()
1924 if (pcbddc->detect_disconnected || matis->allow_repeated) { in PCBDDCComputeLocalTopologyInfo()
1927 PetscBool filter = pcbddc->detect_disconnected_filter; in PCBDDCComputeLocalTopologyInfo()
1929 … for (PetscInt i = 0; i < pcbddc->n_local_subs; i++) PetscCall(ISDestroy(&pcbddc->local_subs[i])); in PCBDDCComputeLocalTopologyInfo()
1930 PetscCall(PetscFree(pcbddc->local_subs)); in PCBDDCComputeLocalTopologyInfo()
1931 PetscCall(MatGetVariableBlockSizes(matis->A, &nel, NULL)); in PCBDDCComputeLocalTopologyInfo()
1932 if (matis->allow_repeated && nel) { in PCBDDCComputeLocalTopologyInfo()
1935 pcbddc->n_local_subs = nel; in PCBDDCComputeLocalTopologyInfo()
1936 PetscCall(MatGetVariableBlockSizes(matis->A, NULL, &elsizes)); in PCBDDCComputeLocalTopologyInfo()
1937 PetscCall(PetscMalloc1(nel, &pcbddc->local_subs)); in PCBDDCComputeLocalTopologyInfo()
1939 PetscCall(ISCreateStride(PETSC_COMM_SELF, elsizes[i], c, 1, &pcbddc->local_subs[i])); in PCBDDCComputeLocalTopologyInfo()
1943 …PetscCall(PCBDDCDetectDisconnectedComponents(pc, filter, &pcbddc->n_local_subs, &pcbddc->local_sub… in PCBDDCComputeLocalTopologyInfo()
1952 PetscCall(MatGetDM(pc->pmat, &dm)); in PCBDDCComputeLocalTopologyInfo()
1975 PetscCall(PetscFree(pcbddc->mat_graph->coords)); in PCBDDCComputeLocalTopologyInfo()
1976 PetscCall(PetscMalloc1(dof * n * cdim, &pcbddc->mat_graph->coords)); in PCBDDCComputeLocalTopologyInfo()
1980 PetscCall(PetscArraycpy(pcbddc->mat_graph->coords, coords, cdim * n * dof)); in PCBDDCComputeLocalTopologyInfo()
1982 PetscReal *bcoords = pcbddc->mat_graph->coords; in PCBDDCComputeLocalTopologyInfo()
1992 pcbddc->mat_graph->cdim = cdim; in PCBDDCComputeLocalTopologyInfo()
1993 pcbddc->mat_graph->cnloc = dof * n; in PCBDDCComputeLocalTopologyInfo()
1994 pcbddc->mat_graph->cloc = PETSC_FALSE; in PCBDDCComputeLocalTopologyInfo()
1997 PetscCall(MatISGetLocalMat(pc->pmat, &lA)); in PCBDDCComputeLocalTopologyInfo()
1999 PetscCall(MatISRestoreLocalMat(pc->pmat, &lA)); in PCBDDCComputeLocalTopologyInfo()
2013 } else { /* the original DMDA local-to-local map have been modified */ in PCBDDCComputeLocalTopologyInfo()
2029 pcbddc->corner_selected = PETSC_TRUE; in PCBDDCComputeLocalTopologyInfo()
2030 pcbddc->corner_selection = PETSC_TRUE; in PCBDDCComputeLocalTopologyInfo()
2036 if (pcbddc->corner_selection && !pcbddc->mat_graph->cdim) { in PCBDDCComputeLocalTopologyInfo()
2039 PetscCall(MatGetDM(pc->pmat, &dm)); in PCBDDCComputeLocalTopologyInfo()
2062 for (d = 1; d < nf; d++) ctxs[d] = ctxs[d - 1] + 1; in PCBDDCComputeLocalTopologyInfo()
2066 …ewer(PetscObjectComm((PetscObject)vcoords), ((PetscObject)vcoords)->options, prefix, "-pc_bddc_coo… in PCBDDCComputeLocalTopologyInfo()
2098 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCConsistencyCheckIS()
2101 PetscInt i, nd, n = matis->A->rmap->n, *nidxs, nnd; in PCBDDCConsistencyCheckIS()
2107 for (i = 0; i < pc->pmat->rmap->n; i++) matis->sf_rootdata[i] = 1; in PCBDDCConsistencyCheckIS()
2109 PetscCall(PetscArrayzero(matis->sf_rootdata, pc->pmat->rmap->n)); in PCBDDCConsistencyCheckIS()
2111 PetscCall(PetscArrayzero(matis->sf_leafdata, n)); in PCBDDCConsistencyCheckIS()
2115 if (-1 < idxs[i] && idxs[i] < n) matis->sf_leafdata[idxs[i]] = 1; in PCBDDCConsistencyCheckIS()
2117 PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, mop)); in PCBDDCConsistencyCheckIS()
2118 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, mop)); in PCBDDCConsistencyCheckIS()
2119 …PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLA… in PCBDDCConsistencyCheckIS()
2120 …PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE… in PCBDDCConsistencyCheckIS()
2127 if (matis->sf_leafdata[i]) nidxs[nnd++] = i; in PCBDDCConsistencyCheckIS()
2136 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCBenignRemoveInterior()
2137 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCBenignRemoveInterior()
2140 if (!pcbddc->benign_have_null) PetscFunctionReturn(PETSC_SUCCESS); in PCBDDCBenignRemoveInterior()
2141 if (pcbddc->ChangeOfBasisMatrix) { in PCBDDCBenignRemoveInterior()
2144 PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change)); in PCBDDCBenignRemoveInterior()
2145 swap = pcbddc->work_change; in PCBDDCBenignRemoveInterior()
2146 pcbddc->work_change = r; in PCBDDCBenignRemoveInterior()
2149 PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCBenignRemoveInterior()
2150 PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCBenignRemoveInterior()
2151 PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0)); in PCBDDCBenignRemoveInterior()
2152 PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D)); in PCBDDCBenignRemoveInterior()
2153 PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0)); in PCBDDCBenignRemoveInterior()
2154 PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D)); in PCBDDCBenignRemoveInterior()
2156 PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCBenignRemoveInterior()
2157 PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCBenignRemoveInterior()
2158 if (pcbddc->ChangeOfBasisMatrix) { in PCBDDCBenignRemoveInterior()
2159 pcbddc->work_change = r; in PCBDDCBenignRemoveInterior()
2160 PetscCall(VecCopy(z, pcbddc->work_change)); in PCBDDCBenignRemoveInterior()
2161 PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z)); in PCBDDCBenignRemoveInterior()
2174 apply_right = ctx->apply_left; in PCBDDCBenignMatMult_Private_Private()
2175 apply_left = ctx->apply_right; in PCBDDCBenignMatMult_Private_Private()
2177 apply_right = ctx->apply_right; in PCBDDCBenignMatMult_Private_Private()
2178 apply_left = ctx->apply_left; in PCBDDCBenignMatMult_Private_Private()
2187 PetscCall(PetscArraycpy(ctx->work, ax, nl)); in PCBDDCBenignMatMult_Private_Private()
2189 for (i = 0; i < ctx->benign_n; i++) { in PCBDDCBenignMatMult_Private_Private()
2193 PetscCall(ISGetLocalSize(ctx->benign_zerodiag_subs[i], &nz)); in PCBDDCBenignMatMult_Private_Private()
2194 PetscCall(ISGetIndices(ctx->benign_zerodiag_subs[i], &idxs)); in PCBDDCBenignMatMult_Private_Private()
2196 if (ctx->apply_p0) { in PCBDDCBenignMatMult_Private_Private()
2197 val = ctx->work[idxs[nz - 1]]; in PCBDDCBenignMatMult_Private_Private()
2198 for (j = 0; j < nz - 1; j++) { in PCBDDCBenignMatMult_Private_Private()
2199 sum += ctx->work[idxs[j]]; in PCBDDCBenignMatMult_Private_Private()
2200 ctx->work[idxs[j]] += val; in PCBDDCBenignMatMult_Private_Private()
2203 for (j = 0; j < nz - 1; j++) sum += ctx->work[idxs[j]]; in PCBDDCBenignMatMult_Private_Private()
2205 ctx->work[idxs[nz - 1]] -= sum; in PCBDDCBenignMatMult_Private_Private()
2206 PetscCall(ISRestoreIndices(ctx->benign_zerodiag_subs[i], &idxs)); in PCBDDCBenignMatMult_Private_Private()
2208 PetscCall(VecPlaceArray(x, ctx->work)); in PCBDDCBenignMatMult_Private_Private()
2212 PetscCall(MatMultTranspose(ctx->A, x, y)); in PCBDDCBenignMatMult_Private_Private()
2214 PetscCall(MatMult(ctx->A, x, y)); in PCBDDCBenignMatMult_Private_Private()
2222 for (i = 0; i < ctx->benign_n; i++) { in PCBDDCBenignMatMult_Private_Private()
2226 PetscCall(ISGetLocalSize(ctx->benign_zerodiag_subs[i], &nz)); in PCBDDCBenignMatMult_Private_Private()
2227 PetscCall(ISGetIndices(ctx->benign_zerodiag_subs[i], &idxs)); in PCBDDCBenignMatMult_Private_Private()
2228 val = -ay[idxs[nz - 1]]; in PCBDDCBenignMatMult_Private_Private()
2229 if (ctx->apply_p0) { in PCBDDCBenignMatMult_Private_Private()
2231 for (j = 0; j < nz - 1; j++) { in PCBDDCBenignMatMult_Private_Private()
2235 ay[idxs[nz - 1]] += sum; in PCBDDCBenignMatMult_Private_Private()
2237 for (j = 0; j < nz - 1; j++) ay[idxs[j]] += val; in PCBDDCBenignMatMult_Private_Private()
2238 ay[idxs[nz - 1]] = 0.; in PCBDDCBenignMatMult_Private_Private()
2240 PetscCall(ISRestoreIndices(ctx->benign_zerodiag_subs[i], &idxs)); in PCBDDCBenignMatMult_Private_Private()
2263 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCBenignShellMat()
2264 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCBenignShellMat()
2271 PCBDDCReuseSolvers reuse = pcbddc->sub_schurs ? pcbddc->sub_schurs->reuse_solver : NULL; in PCBDDCBenignShellMat()
2273 …PetscCheck(!pcbddc->benign_original_mat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Benign original mat has… in PCBDDCBenignShellMat()
2274 …if (!pcbddc->benign_change || !pcbddc->benign_n || pcbddc->benign_change_explicit) PetscFunctionRe… in PCBDDCBenignShellMat()
2275 PetscCall(PetscMalloc1(pcis->n, &work)); in PCBDDCBenignShellMat()
2277 PetscCall(MatSetSizes(A_IB, pcis->n - pcis->n_B, pcis->n_B, PETSC_DECIDE, PETSC_DECIDE)); in PCBDDCBenignShellMat()
2283 ctx->apply_left = PETSC_TRUE; in PCBDDCBenignShellMat()
2284 ctx->apply_right = PETSC_FALSE; in PCBDDCBenignShellMat()
2285 ctx->apply_p0 = PETSC_FALSE; in PCBDDCBenignShellMat()
2286 ctx->benign_n = pcbddc->benign_n; in PCBDDCBenignShellMat()
2288 ctx->benign_zerodiag_subs = reuse->benign_zerodiag_subs; in PCBDDCBenignShellMat()
2289 ctx->free = PETSC_FALSE; in PCBDDCBenignShellMat()
2294 PetscCall(ISLocalToGlobalMappingCreateIS(pcis->is_I_local, &N_to_D)); in PCBDDCBenignShellMat()
2295 PetscCall(PetscMalloc1(pcbddc->benign_n, &ctx->benign_zerodiag_subs)); in PCBDDCBenignShellMat()
2296 …i = 0; i < pcbddc->benign_n; i++) PetscCall(ISGlobalToLocalMappingApplyIS(N_to_D, IS_GTOLM_DROP, p… in PCBDDCBenignShellMat()
2298 ctx->free = PETSC_TRUE; in PCBDDCBenignShellMat()
2300 ctx->A = pcis->A_IB; in PCBDDCBenignShellMat()
2301 ctx->work = work; in PCBDDCBenignShellMat()
2305 pcis->A_IB = A_IB; in PCBDDCBenignShellMat()
2309 pcbddc->benign_original_mat = pcis->A_BI; in PCBDDCBenignShellMat()
2310 pcis->A_BI = A_BI; in PCBDDCBenignShellMat()
2312 if (!pcbddc->benign_original_mat) PetscFunctionReturn(PETSC_SUCCESS); in PCBDDCBenignShellMat()
2313 PetscCall(MatShellGetContext(pcis->A_IB, &ctx)); in PCBDDCBenignShellMat()
2314 PetscCall(MatDestroy(&pcis->A_IB)); in PCBDDCBenignShellMat()
2315 pcis->A_IB = ctx->A; in PCBDDCBenignShellMat()
2316 ctx->A = NULL; in PCBDDCBenignShellMat()
2317 PetscCall(MatDestroy(&pcis->A_BI)); in PCBDDCBenignShellMat()
2318 pcis->A_BI = pcbddc->benign_original_mat; in PCBDDCBenignShellMat()
2319 pcbddc->benign_original_mat = NULL; in PCBDDCBenignShellMat()
2320 if (ctx->free) { in PCBDDCBenignShellMat()
2322 for (i = 0; i < ctx->benign_n; i++) PetscCall(ISDestroy(&ctx->benign_zerodiag_subs[i])); in PCBDDCBenignShellMat()
2323 PetscCall(PetscFree(ctx->benign_zerodiag_subs)); in PCBDDCBenignShellMat()
2325 PetscCall(PetscFree(ctx->work)); in PCBDDCBenignShellMat()
2334 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCBenignProject()
2335 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCBenignProject()
2339 PetscCall(MatPtAP(matis->A, pcbddc->benign_change, MAT_INITIAL_MATRIX, 2.0, &An)); in PCBDDCBenignProject()
2340 PetscCall(MatZeroRowsColumns(An, pcbddc->benign_n, pcbddc->benign_p0_lidx, 1.0, NULL, NULL)); in PCBDDCBenignProject()
2388 Mat_SeqAIJ *b = (Mat_SeqAIJ *)Bt->data; in MatSeqAIJCompress()
2389 b->free_a = PETSC_TRUE; in MatSeqAIJCompress()
2390 b->free_ij = PETSC_TRUE; in MatSeqAIJCompress()
2414 PetscCall(MatGetDM(pc->pmat, &dm)); in PCBDDCDetectDisconnectedComponents()
2443 /* Skip non-owned cells in parallel (ParMetis expects no overlap) */ in PCBDDCDetectDisconnectedComponents()
2477 graph->xadj = xadj; in PCBDDCDetectDisconnectedComponents()
2478 graph->adjncy = adjncy; in PCBDDCDetectDisconnectedComponents()
2483 PetscCall(MatISGetLocalMat(pc->pmat, &A)); in PCBDDCDetectDisconnectedComponents()
2484 if (!A->rmap->N || !A->cmap->N) { in PCBDDCDetectDisconnectedComponents()
2497 PetscReal chop = 1.e-6; in PCBDDCDetectDisconnectedComponents()
2539 graph->xadj = xadj_filtered; in PCBDDCDetectDisconnectedComponents()
2540 graph->adjncy = adjncy_filtered; in PCBDDCDetectDisconnectedComponents()
2542 graph->xadj = xadj; in PCBDDCDetectDisconnectedComponents()
2543 graph->adjncy = adjncy; in PCBDDCDetectDisconnectedComponents()
2547 …graph->seq_graph = PETSC_TRUE; /* analyze local connected components (i.e. disconnected subdomains… in PCBDDCDetectDisconnectedComponents()
2570 if (ncc) *ncc = graph->ncc; in PCBDDCDetectDisconnectedComponents()
2584 PetscCall(MatISGetLocalMat(pc->pmat, &A)); in PCBDDCDetectDisconnectedComponents()
2585 PetscCall(PetscMalloc3(A->rmap->n, &ids, graph->ncc + 1, &cids, A->rmap->n, &pids)); in PCBDDCDetectDisconnectedComponents()
2586 PetscCall(PetscBTCreate(A->rmap->n, &btv)); in PCBDDCDetectDisconnectedComponents()
2587 PetscCall(PetscBTCreate(A->rmap->n, &btvt)); in PCBDDCDetectDisconnectedComponents()
2588 PetscCall(PetscBTCreate(pEnd - pStart, &btvc)); in PCBDDCDetectDisconnectedComponents()
2634 if (cnt >= mcnt) PetscCall(PetscBTSet(btvc, v - pStart)); in PCBDDCDetectDisconnectedComponents()
2645 for (i = 0, cump = 0, cum = 0; i < graph->ncc; i++) { in PCBDDCDetectDisconnectedComponents()
2648 PetscCall(PetscBTMemzero(A->rmap->n, btvt)); in PCBDDCDetectDisconnectedComponents()
2649 for (j = graph->cptr[i]; j < graph->cptr[i + 1]; j++) { in PCBDDCDetectDisconnectedComponents()
2650 PetscInt k, size, *closure = NULL, cell = graph->queue[j]; in PCBDDCDetectDisconnectedComponents()
2659 for (s = 0; s < dof - cdof; s++) { in PCBDDCDetectDisconnectedComponents()
2661 if (PetscBTLookup(btvc, p - pStart)) pids[cump++] = off + s; /* subdomain corner */ in PCBDDCDetectDisconnectedComponents()
2663 else pids[cump++] = off + s; /* cross-vertex */ in PCBDDCDetectDisconnectedComponents()
2670 for (s = 0; s < dof - cdof; s++) { in PCBDDCDetectDisconnectedComponents()
2672 if (PetscBTLookup(btvc, pp - pStart)) pids[cump++] = off + s; /* subdomain corner */ in PCBDDCDetectDisconnectedComponents()
2674 else pids[cump++] = off + s; /* cross-vertex */ in PCBDDCDetectDisconnectedComponents()
2685 PetscCall(PetscMalloc1(graph->ncc, &cc_n)); in PCBDDCDetectDisconnectedComponents()
2686 …for (i = 0; i < graph->ncc; i++) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, cids[i + 1] - cids[i],… in PCBDDCDetectDisconnectedComponents()
2697 if (ncc) *ncc = graph->ncc; in PCBDDCDetectDisconnectedComponents()
2699 PetscCall(PetscMalloc1(graph->ncc, &cc_n)); in PCBDDCDetectDisconnectedComponents()
2700 …or (i = 0; i < graph->ncc; i++) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, graph->cptr[i + 1] - gr… in PCBDDCDetectDisconnectedComponents()
2705 graph->xadj = NULL; in PCBDDCDetectDisconnectedComponents()
2706 graph->adjncy = NULL; in PCBDDCDetectDisconnectedComponents()
2713 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCBenignCheck()
2714 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCBenignCheck()
2719 PetscCall(PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirIS)); in PCBDDCBenignCheck()
2728 PetscCall(VecSet(pcis->vec1_N, 0.)); in PCBDDCBenignCheck()
2729 PetscCall(PetscMalloc1(pcis->n, &vals)); in PCBDDCBenignCheck()
2733 PetscCall(VecSetValues(pcis->vec1_N, nz, idxs, vals, INSERT_VALUES)); in PCBDDCBenignCheck()
2734 PetscCall(VecAssemblyBegin(pcis->vec1_N)); in PCBDDCBenignCheck()
2735 PetscCall(VecAssemblyEnd(pcis->vec1_N)); in PCBDDCBenignCheck()
2737 PetscCall(VecSetRandom(pcis->vec2_N, NULL)); in PCBDDCBenignCheck()
2739 PetscCall(VecSetValues(pcis->vec2_N, nz, idxs, vals, INSERT_VALUES)); in PCBDDCBenignCheck()
2741 PetscCall(ISGetIndices(pcis->is_B_local, &idxs)); in PCBDDCBenignCheck()
2742 for (i = 0; i < pcis->n_B; i++) vals[i] = 0.; in PCBDDCBenignCheck()
2743 PetscCall(VecSetValues(pcis->vec2_N, pcis->n_B, idxs, vals, INSERT_VALUES)); in PCBDDCBenignCheck()
2744 PetscCall(ISRestoreIndices(pcis->is_B_local, &idxs)); in PCBDDCBenignCheck()
2751 PetscCall(VecSetValues(pcis->vec2_N, n, idxs, vals, INSERT_VALUES)); in PCBDDCBenignCheck()
2754 PetscCall(VecAssemblyBegin(pcis->vec2_N)); in PCBDDCBenignCheck()
2755 PetscCall(VecAssemblyEnd(pcis->vec2_N)); in PCBDDCBenignCheck()
2756 PetscCall(VecDuplicate(pcis->vec1_N, &vec3_N)); in PCBDDCBenignCheck()
2758 PetscCall(MatISGetLocalMat(pc->pmat, &A)); in PCBDDCBenignCheck()
2759 PetscCall(MatMult(A, pcis->vec1_N, vec3_N)); in PCBDDCBenignCheck()
2760 PetscCall(VecDot(vec3_N, pcis->vec2_N, &vals[0])); in PCBDDCBenignCheck()
2761 …PetscCheck(PetscAbsScalar(vals[0]) <= 1.e-1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Benign trick can not… in PCBDDCBenignCheck()
2766 PetscCall(PetscCalloc1(pcis->n, &count)); in PCBDDCBenignCheck()
2767 PetscCall(ISGetIndices(pcis->is_B_local, &idxs)); in PCBDDCBenignCheck()
2768 for (i = 0; i < pcis->n_B; i++) count[idxs[i]]++; in PCBDDCBenignCheck()
2769 PetscCall(ISRestoreIndices(pcis->is_B_local, &idxs)); in PCBDDCBenignCheck()
2778 PetscCall(VecSetRandom(pcis->vec1_global, NULL)); in PCBDDCBenignCheck()
2779 for (i = 0; i < pcbddc->benign_n; i++) pcbddc->benign_p0[i] = -PetscGlobalRank - i; in PCBDDCBenignCheck()
2780 PetscCall(PCBDDCBenignGetOrSetP0(pc, pcis->vec1_global, PETSC_FALSE)); in PCBDDCBenignCheck()
2781 for (i = 0; i < pcbddc->benign_n; i++) pcbddc->benign_p0[i] = 1; in PCBDDCBenignCheck()
2782 PetscCall(PCBDDCBenignGetOrSetP0(pc, pcis->vec1_global, PETSC_TRUE)); in PCBDDCBenignCheck()
2783 for (i = 0; i < pcbddc->benign_n; i++) { in PCBDDCBenignCheck()
2784 PetscInt val = PetscRealPart(pcbddc->benign_p0[i]); in PCBDDCBenignCheck()
2785 …-PetscGlobalRank - i, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error testing PCBDDCBenignGetOrSetP0! Foun… in PCBDDCBenignCheck()
2792 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCBenignDetectSaddlePoint()
2793 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCBenignDetectSaddlePoint()
2801 PetscCall(PetscSFDestroy(&pcbddc->benign_sf)); in PCBDDCBenignDetectSaddlePoint()
2802 PetscCall(MatDestroy(&pcbddc->benign_B0)); in PCBDDCBenignDetectSaddlePoint()
2803 for (n = 0; n < pcbddc->benign_n; n++) PetscCall(ISDestroy(&pcbddc->benign_zerodiag_subs[n])); in PCBDDCBenignDetectSaddlePoint()
2804 PetscCall(PetscFree(pcbddc->benign_zerodiag_subs)); in PCBDDCBenignDetectSaddlePoint()
2813 if (pcbddc->n_ISForDofsLocal) { in PCBDDCBenignDetectSaddlePoint()
2818 PetscCall(PetscMalloc1(pcbddc->n_ISForDofsLocal, &pp)); in PCBDDCBenignDetectSaddlePoint()
2819 n = pcbddc->n_ISForDofsLocal; in PCBDDCBenignDetectSaddlePoint()
2820 …PetscOptionsBegin(PetscObjectComm((PetscObject)pc), ((PetscObject)pc)->prefix, "BDDC benign option… in PCBDDCBenignDetectSaddlePoint()
2821 …PetscCall(PetscOptionsIntArray("-pc_bddc_pressure_field", "Field id for pressures", NULL, pp, &n, … in PCBDDCBenignDetectSaddlePoint()
2822 …PetscCall(PetscOptionsBool("-pc_bddc_pressure_blocked", "Use blocked pressure fields", NULL, block… in PCBDDCBenignDetectSaddlePoint()
2826 pp[0] = pcbddc->n_ISForDofsLocal - 1; in PCBDDCBenignDetectSaddlePoint()
2833 …PetscCheck(pp[p] >= 0 && pp[p] < pcbddc->n_ISForDofsLocal, PetscObjectComm((PetscObject)pc), PETSC… in PCBDDCBenignDetectSaddlePoint()
2834 if (blocked) PetscCall(ISGetBlockSize(pcbddc->ISForDofsLocal[pp[p]], &bs)); in PCBDDCBenignDetectSaddlePoint()
2843 if (blocked) PetscCall(ISGetBlockSize(pcbddc->ISForDofsLocal[pp[p]], &bs)); in PCBDDCBenignDetectSaddlePoint()
2844 PetscCall(ISGetLocalSize(pcbddc->ISForDofsLocal[pp[p]], &npl)); in PCBDDCBenignDetectSaddlePoint()
2845 PetscCall(ISGetIndices(pcbddc->ISForDofsLocal[pp[p]], &idxs)); in PCBDDCBenignDetectSaddlePoint()
2855 PetscCall(ISRestoreIndices(pcbddc->ISForDofsLocal[pp[p]], &idxs)); in PCBDDCBenignDetectSaddlePoint()
2859 /* remove zeroed out pressures if we are setting up a BDDC solver for a saddle-point FETI-DP */ in PCBDDCBenignDetectSaddlePoint()
2874 PetscCall(MatGetLocalSize(pcbddc->local_mat, &n, NULL)); in PCBDDCBenignDetectSaddlePoint()
2875 if (!n) pcbddc->benign_change_explicit = PETSC_TRUE; in PCBDDCBenignDetectSaddlePoint()
2876 PetscCall(MatFindZeroDiagonals(pcbddc->local_mat, &zerodiag)); in PCBDDCBenignDetectSaddlePoint()
2895 …if (pcbddc->NeumannBoundariesLocal) PetscCall(ISGetLocalSize(pcbddc->NeumannBoundariesLocal, &nneu… in PCBDDCBenignDetectSaddlePoint()
2896 checkb = (PetscBool)(!pcbddc->NeumannBoundariesLocal || pcbddc->current_level); in PCBDDCBenignDetectSaddlePoint()
2902 PetscCall(MatISGetLocalToGlobalMapping(pc->pmat, &mapping, NULL)); in PCBDDCBenignDetectSaddlePoint()
2916 subs = pcbddc->local_subs; in PCBDDCBenignDetectSaddlePoint()
2917 nsubs = pcbddc->n_local_subs; in PCBDDCBenignDetectSaddlePoint()
2920 PetscCall(VecDuplicateVecs(matis->y, 2, &work)); in PCBDDCBenignDetectSaddlePoint()
2955 PetscCall(MatGetLocalSize(pcbddc->local_mat, &nl, NULL)); in PCBDDCBenignDetectSaddlePoint()
2966 PetscCall(VecSet(matis->x, 0)); in PCBDDCBenignDetectSaddlePoint()
2969 PetscCall(VecGetArray(matis->x, &array)); in PCBDDCBenignDetectSaddlePoint()
2971 PetscCall(VecRestoreArray(matis->x, &array)); in PCBDDCBenignDetectSaddlePoint()
2973 PetscCall(VecPointwiseMult(matis->x, work[0], matis->x)); in PCBDDCBenignDetectSaddlePoint()
2974 PetscCall(MatMult(matis->A, matis->x, matis->y)); in PCBDDCBenignDetectSaddlePoint()
2975 PetscCall(VecPointwiseMult(matis->y, work[1], matis->y)); in PCBDDCBenignDetectSaddlePoint()
2976 PetscCall(VecGetArray(matis->y, &array)); in PCBDDCBenignDetectSaddlePoint()
2983 PetscCall(VecRestoreArray(matis->y, &array)); in PCBDDCBenignDetectSaddlePoint()
2989 PetscCall(ISGetIndices(pcbddc->NeumannBoundariesLocal, &idxs)); in PCBDDCBenignDetectSaddlePoint()
2991 PetscCall(ISRestoreIndices(pcbddc->NeumannBoundariesLocal, &idxs)); in PCBDDCBenignDetectSaddlePoint()
3021 PetscCall(MatMult(matis->A, work[0], matis->x)); in PCBDDCBenignDetectSaddlePoint()
3022 PetscCall(VecPointwiseMult(matis->x, work[1], matis->x)); in PCBDDCBenignDetectSaddlePoint()
3023 PetscCall(VecGetArray(matis->x, &array)); in PCBDDCBenignDetectSaddlePoint()
3030 PetscCall(VecRestoreArray(matis->x, &array)); in PCBDDCBenignDetectSaddlePoint()
3048 PetscCall(MatGetLocalSize(pcbddc->local_mat, &n, NULL)); in PCBDDCBenignDetectSaddlePoint()
3086 …PetscCallMPI(MPIU_Allreduce(&have_null, &pcbddc->benign_null, 1, MPI_C_BOOL, MPI_LAND, PetscObject… in PCBDDCBenignDetectSaddlePoint()
3088 /* Prepare matrix to compute no-net-flux */ in PCBDDCBenignDetectSaddlePoint()
3089 if (pcbddc->compute_nonetflux && !pcbddc->divudotp) { in PCBDDCBenignDetectSaddlePoint()
3100 PetscCall(MatISGetLocalToGlobalMapping(pc->pmat, &l2gmap, NULL)); in PCBDDCBenignDetectSaddlePoint()
3101 PetscCall(MatISGetLocalMat(pc->pmat, &A)); in PCBDDCBenignDetectSaddlePoint()
3107 st = st - n_isused; in PCBDDCBenignDetectSaddlePoint()
3122 PetscCall(MatGetSize(pc->pmat, NULL, &N)); in PCBDDCBenignDetectSaddlePoint()
3128 PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pcbddc->divudotp)); in PCBDDCBenignDetectSaddlePoint()
3129 PetscCall(MatSetType(pcbddc->divudotp, MATIS)); in PCBDDCBenignDetectSaddlePoint()
3130 PetscCall(MatSetSizes(pcbddc->divudotp, PETSC_DECIDE, PETSC_DECIDE, M, N)); in PCBDDCBenignDetectSaddlePoint()
3131 PetscCall(MatSetLocalToGlobalMapping(pcbddc->divudotp, rl2g, cl2g)); in PCBDDCBenignDetectSaddlePoint()
3134 PetscCall(MatISSetLocalMat(pcbddc->divudotp, loc_divudotp)); in PCBDDCBenignDetectSaddlePoint()
3136 PetscCall(MatAssemblyBegin(pcbddc->divudotp, MAT_FINAL_ASSEMBLY)); in PCBDDCBenignDetectSaddlePoint()
3137 PetscCall(MatAssemblyEnd(pcbddc->divudotp, MAT_FINAL_ASSEMBLY)); in PCBDDCBenignDetectSaddlePoint()
3147 pcbddc->benign_n = benign_n; in PCBDDCBenignDetectSaddlePoint()
3148 pcbddc->benign_zerodiag_subs = zerodiag_subs; in PCBDDCBenignDetectSaddlePoint()
3151 have_null = (PetscBool)(!!pcbddc->benign_n); in PCBDDCBenignDetectSaddlePoint()
3152 …PetscCallMPI(MPIU_Allreduce(&have_null, &pcbddc->benign_have_null, 1, MPI_C_BOOL, MPI_LOR, PetscOb… in PCBDDCBenignDetectSaddlePoint()
3155 PetscCall(MatGetLocalSize(pcbddc->local_mat, &n, NULL)); in PCBDDCBenignDetectSaddlePoint()
3157 if (pcbddc->benign_n) { in PCBDDCBenignDetectSaddlePoint()
3161 PetscCall(MatDestroy(&pcbddc->benign_change)); in PCBDDCBenignDetectSaddlePoint()
3162 PetscCall(MatCreate(PetscObjectComm((PetscObject)pcbddc->local_mat), &pcbddc->benign_change)); in PCBDDCBenignDetectSaddlePoint()
3163 PetscCall(MatSetType(pcbddc->benign_change, MATAIJ)); in PCBDDCBenignDetectSaddlePoint()
3164 PetscCall(MatSetSizes(pcbddc->benign_change, n, n, PETSC_DECIDE, PETSC_DECIDE)); in PCBDDCBenignDetectSaddlePoint()
3167 for (i = 0; i < pcbddc->benign_n; i++) { in PCBDDCBenignDetectSaddlePoint()
3171 PetscCall(ISGetLocalSize(pcbddc->benign_zerodiag_subs[i], &nzs)); in PCBDDCBenignDetectSaddlePoint()
3172 PetscCall(ISGetIndices(pcbddc->benign_zerodiag_subs[i], &idxs)); in PCBDDCBenignDetectSaddlePoint()
3173 for (j = 0; j < nzs - 1; j++) nnz[idxs[j]] = 2; /* change on pressures */ in PCBDDCBenignDetectSaddlePoint()
3174 nnz[idxs[nzs - 1]] = nzs; /* last local pressure dof in subdomain */ in PCBDDCBenignDetectSaddlePoint()
3175 PetscCall(ISRestoreIndices(pcbddc->benign_zerodiag_subs[i], &idxs)); in PCBDDCBenignDetectSaddlePoint()
3177 PetscCall(MatSeqAIJSetPreallocation(pcbddc->benign_change, 0, nnz)); in PCBDDCBenignDetectSaddlePoint()
3178 PetscCall(MatSetOption(pcbddc->benign_change, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); in PCBDDCBenignDetectSaddlePoint()
3181 for (i = 0; i < n; i++) PetscCall(MatSetValue(pcbddc->benign_change, i, i, 1., INSERT_VALUES)); in PCBDDCBenignDetectSaddlePoint()
3182 PetscCall(PetscFree3(pcbddc->benign_p0_lidx, pcbddc->benign_p0_gidx, pcbddc->benign_p0)); in PCBDDCBenignDetectSaddlePoint()
3183 …all(PetscMalloc3(pcbddc->benign_n, &pcbddc->benign_p0_lidx, pcbddc->benign_n, &pcbddc->benign_p0_g… in PCBDDCBenignDetectSaddlePoint()
3185 for (s = 0; s < pcbddc->benign_n; s++) { in PCBDDCBenignDetectSaddlePoint()
3190 PetscCall(ISGetLocalSize(pcbddc->benign_zerodiag_subs[s], &nzs)); in PCBDDCBenignDetectSaddlePoint()
3191 PetscCall(ISGetIndices(pcbddc->benign_zerodiag_subs[s], &idxs)); in PCBDDCBenignDetectSaddlePoint()
3192 for (i = 0; i < nzs - 1; i++) { in PCBDDCBenignDetectSaddlePoint()
3197 cols[1] = idxs[nzs - 1]; in PCBDDCBenignDetectSaddlePoint()
3200 PetscCall(MatSetValues(pcbddc->benign_change, 1, cols, 2, cols, vals, INSERT_VALUES)); in PCBDDCBenignDetectSaddlePoint()
3203 for (i = 0; i < nzs - 1; i++) array[i] = -1.; in PCBDDCBenignDetectSaddlePoint()
3204 array[nzs - 1] = 1.; in PCBDDCBenignDetectSaddlePoint()
3205 …PetscCall(MatSetValues(pcbddc->benign_change, 1, idxs + nzs - 1, nzs, idxs, array, INSERT_VALUES)); in PCBDDCBenignDetectSaddlePoint()
3207 pcbddc->benign_p0_lidx[s] = idxs[nzs - 1]; in PCBDDCBenignDetectSaddlePoint()
3208 PetscCall(ISRestoreIndices(pcbddc->benign_zerodiag_subs[s], &idxs)); in PCBDDCBenignDetectSaddlePoint()
3211 PetscCall(MatAssemblyBegin(pcbddc->benign_change, MAT_FINAL_ASSEMBLY)); in PCBDDCBenignDetectSaddlePoint()
3212 PetscCall(MatAssemblyEnd(pcbddc->benign_change, MAT_FINAL_ASSEMBLY)); in PCBDDCBenignDetectSaddlePoint()
3215 if (pcbddc->benign_change_explicit) { in PCBDDCBenignDetectSaddlePoint()
3218 PetscCall(MatPtAP(pcbddc->local_mat, pcbddc->benign_change, MAT_INITIAL_MATRIX, 2.0, &M)); in PCBDDCBenignDetectSaddlePoint()
3219 PetscCall(MatDestroy(&pcbddc->local_mat)); in PCBDDCBenignDetectSaddlePoint()
3220 PetscCall(MatSeqAIJCompress(M, &pcbddc->local_mat)); in PCBDDCBenignDetectSaddlePoint()
3224 …PetscCall(ISLocalToGlobalMappingApply(matis->rmapping, pcbddc->benign_n, pcbddc->benign_p0_lidx, p… in PCBDDCBenignDetectSaddlePoint()
3232 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCBenignGetOrSetP0()
3236 if (!pcbddc->benign_sf) { in PCBDDCBenignGetOrSetP0()
3237 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)pc), &pcbddc->benign_sf)); in PCBDDCBenignGetOrSetP0()
3238 …etscCall(PetscSFSetGraphLayout(pcbddc->benign_sf, pc->pmat->rmap, pcbddc->benign_n, NULL, PETSC_OW… in PCBDDCBenignGetOrSetP0()
3242 …PetscCall(PetscSFBcastBegin(pcbddc->benign_sf, MPIU_SCALAR, array, pcbddc->benign_p0, MPI_REPLACE)… in PCBDDCBenignGetOrSetP0()
3243 … PetscCall(PetscSFBcastEnd(pcbddc->benign_sf, MPIU_SCALAR, array, pcbddc->benign_p0, MPI_REPLACE)); in PCBDDCBenignGetOrSetP0()
3247 …PetscCall(PetscSFReduceBegin(pcbddc->benign_sf, MPIU_SCALAR, pcbddc->benign_p0, array, MPI_REPLACE… in PCBDDCBenignGetOrSetP0()
3248 …PetscCall(PetscSFReduceEnd(pcbddc->benign_sf, MPIU_SCALAR, pcbddc->benign_p0, array, MPI_REPLACE)); in PCBDDCBenignGetOrSetP0()
3256 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCBenignPopOrPushB0()
3260 - avoid nested pop (or push) calls. in PCBDDCBenignPopOrPushB0()
3261 - cannot push before pop. in PCBDDCBenignPopOrPushB0()
3262 - cannot call this if pcbddc->local_mat is NULL in PCBDDCBenignPopOrPushB0()
3264 if (!pcbddc->benign_n) PetscFunctionReturn(PETSC_SUCCESS); in PCBDDCBenignPopOrPushB0()
3266 if (pcbddc->benign_change_explicit) { in PCBDDCBenignPopOrPushB0()
3272 if (pcbddc->benign_B0) reuse = MAT_REUSE_MATRIX; in PCBDDCBenignPopOrPushB0()
3273 …PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pcbddc->benign_n, pcbddc->benign_p0_lidx, PETSC_COPY_VA… in PCBDDCBenignPopOrPushB0()
3274 PetscCall(MatCreateSubMatrix(pcbddc->local_mat, is_p0, NULL, reuse, &pcbddc->benign_B0)); in PCBDDCBenignPopOrPushB0()
3276 PetscCall(MatSetOption(pcbddc->local_mat, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); in PCBDDCBenignPopOrPushB0()
3277 PetscCall(MatSetOption(pcbddc->local_mat, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)); in PCBDDCBenignPopOrPushB0()
3278 PetscCall(MatZeroRowsColumnsIS(pcbddc->local_mat, is_p0, 1.0, NULL, NULL)); in PCBDDCBenignPopOrPushB0()
3281 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCBenignPopOrPushB0()
3285 PetscCall(VecGetLocalSize(matis->y, &n)); in PCBDDCBenignPopOrPushB0()
3287 if (!pcbddc->benign_B0) { in PCBDDCBenignPopOrPushB0()
3289 PetscCall(MatCreate(PetscObjectComm((PetscObject)pcbddc->local_mat), &pcbddc->benign_B0)); in PCBDDCBenignPopOrPushB0()
3290 PetscCall(MatSetType(pcbddc->benign_B0, MATAIJ)); in PCBDDCBenignPopOrPushB0()
3291 PetscCall(MatSetSizes(pcbddc->benign_B0, pcbddc->benign_n, n, PETSC_DECIDE, PETSC_DECIDE)); in PCBDDCBenignPopOrPushB0()
3292 PetscCall(PetscMalloc1(pcbddc->benign_n, &nnz)); in PCBDDCBenignPopOrPushB0()
3293 for (i = 0; i < pcbddc->benign_n; i++) { in PCBDDCBenignPopOrPushB0()
3294 PetscCall(ISGetLocalSize(pcbddc->benign_zerodiag_subs[i], &nnz[i])); in PCBDDCBenignPopOrPushB0()
3295 nnz[i] = n - nnz[i]; in PCBDDCBenignPopOrPushB0()
3297 PetscCall(MatSeqAIJSetPreallocation(pcbddc->benign_B0, 0, nnz)); in PCBDDCBenignPopOrPushB0()
3298 PetscCall(MatSetOption(pcbddc->benign_B0, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); in PCBDDCBenignPopOrPushB0()
3302 for (i = 0; i < pcbddc->benign_n; i++) { in PCBDDCBenignPopOrPushB0()
3306 PetscCall(VecSet(matis->x, 0.)); in PCBDDCBenignPopOrPushB0()
3307 PetscCall(ISGetLocalSize(pcbddc->benign_zerodiag_subs[i], &nz)); in PCBDDCBenignPopOrPushB0()
3308 PetscCall(ISGetIndices(pcbddc->benign_zerodiag_subs[i], (const PetscInt **)&idxs)); in PCBDDCBenignPopOrPushB0()
3310 PetscCall(VecSetValues(matis->x, nz, idxs, vals, INSERT_VALUES)); in PCBDDCBenignPopOrPushB0()
3311 PetscCall(VecAssemblyBegin(matis->x)); in PCBDDCBenignPopOrPushB0()
3312 PetscCall(VecAssemblyEnd(matis->x)); in PCBDDCBenignPopOrPushB0()
3313 PetscCall(VecSet(matis->y, 0.)); in PCBDDCBenignPopOrPushB0()
3314 PetscCall(MatMult(matis->A, matis->x, matis->y)); in PCBDDCBenignPopOrPushB0()
3315 PetscCall(VecGetArray(matis->y, &array)); in PCBDDCBenignPopOrPushB0()
3324 PetscCall(MatSetValues(pcbddc->benign_B0, 1, &i, cum, idxs_ins, vals, INSERT_VALUES)); in PCBDDCBenignPopOrPushB0()
3325 PetscCall(VecRestoreArray(matis->y, &array)); in PCBDDCBenignPopOrPushB0()
3326 PetscCall(ISRestoreIndices(pcbddc->benign_zerodiag_subs[i], (const PetscInt **)&idxs)); in PCBDDCBenignPopOrPushB0()
3328 PetscCall(MatAssemblyBegin(pcbddc->benign_B0, MAT_FINAL_ASSEMBLY)); in PCBDDCBenignPopOrPushB0()
3329 PetscCall(MatAssemblyEnd(pcbddc->benign_B0, MAT_FINAL_ASSEMBLY)); in PCBDDCBenignPopOrPushB0()
3334 PetscCheck(pcbddc->benign_change_explicit, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot push B0!"); in PCBDDCBenignPopOrPushB0()
3335 for (PetscInt i = 0; i < pcbddc->benign_n; i++) { in PCBDDCBenignPopOrPushB0()
3339 …PetscCall(MatGetRow(pcbddc->benign_B0, i, &B0_ncol, (const PetscInt **)&B0_cols, (const PetscScala… in PCBDDCBenignPopOrPushB0()
3340 …PetscCall(MatSetValues(pcbddc->local_mat, 1, pcbddc->benign_p0_lidx + i, B0_ncol, B0_cols, B0_vals… in PCBDDCBenignPopOrPushB0()
3341 …PetscCall(MatSetValues(pcbddc->local_mat, B0_ncol, B0_cols, 1, pcbddc->benign_p0_lidx + i, B0_vals… in PCBDDCBenignPopOrPushB0()
3342 …PetscCall(MatSetValue(pcbddc->local_mat, pcbddc->benign_p0_lidx[i], pcbddc->benign_p0_lidx[i], 0.0… in PCBDDCBenignPopOrPushB0()
3343 …PetscCall(MatRestoreRow(pcbddc->benign_B0, i, &B0_ncol, (const PetscInt **)&B0_cols, (const PetscS… in PCBDDCBenignPopOrPushB0()
3345 PetscCall(MatAssemblyBegin(pcbddc->local_mat, MAT_FINAL_ASSEMBLY)); in PCBDDCBenignPopOrPushB0()
3346 PetscCall(MatAssemblyEnd(pcbddc->local_mat, MAT_FINAL_ASSEMBLY)); in PCBDDCBenignPopOrPushB0()
3353 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCAdaptiveSelection()
3354 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCAdaptiveSelection()
3368 if (!pcbddc->adaptive_selection) PetscFunctionReturn(PETSC_SUCCESS); in PCBDDCAdaptiveSelection()
3370 …PetscCheck(sub_schurs->schur_explicit || !sub_schurs->n_subs, PetscObjectComm((PetscObject)pc), PE… in PCBDDCAdaptiveSelection()
3371 …->n_subs || sub_schurs->is_symmetric, PETSC_COMM_SELF, PETSC_ERR_SUP, "Adaptive selection not yet … in PCBDDCAdaptiveSelection()
3372 sub_schurs->is_posdef); in PCBDDCAdaptiveSelection()
3373 PetscCall(PetscLogEventBegin(PC_BDDC_AdaptiveSetUp[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCAdaptiveSelection()
3375 if (pcbddc->dbg_flag) { in PCBDDCAdaptiveSelection()
3376 …if (!pcbddc->dbg_viewer) pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc… in PCBDDCAdaptiveSelection()
3377 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCAdaptiveSelection()
3378 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "--------------------------------------------… in PCBDDCAdaptiveSelection()
3379 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "Check adaptive selection of constraints\n")); in PCBDDCAdaptiveSelection()
3380 PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer)); in PCBDDCAdaptiveSelection()
3383 …->dbg_flag) PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d cc %"… in PCBDDCAdaptiveSelection()
3387 for (i = 0; i < sub_schurs->n_subs; i++) { in PCBDDCAdaptiveSelection()
3390 PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size)); in PCBDDCAdaptiveSelection()
3395 nmax = pcbddc->adaptive_nmax > 0 ? pcbddc->adaptive_nmax : mss; in PCBDDCAdaptiveSelection()
3396 nmin = pcbddc->adaptive_nmin > 0 ? pcbddc->adaptive_nmin : 0; in PCBDDCAdaptiveSelection()
3399 if (nmin || !sub_schurs->is_posdef) { /* XXX */ in PCBDDCAdaptiveSelection()
3406 for (i = 0; i < sub_schurs->n_subs; i++) { in PCBDDCAdaptiveSelection()
3409 PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size)); in PCBDDCAdaptiveSelection()
3423 PetscCheck(sub_schurs->is_symmetric, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented"); in PCBDDCAdaptiveSelection()
3425 B_lwork = -1; in PCBDDCAdaptiveSelection()
3448 …if (sub_schurs->is_vertices && pcbddc->use_vertices) { /* complement set of active subsets, each e… in PCBDDCAdaptiveSelection()
3449 PetscCall(ISGetLocalSize(sub_schurs->is_vertices, &nv)); in PCBDDCAdaptiveSelection()
3457 …->n_subs, &pcbddc->adaptive_constraints_n, nv + sub_schurs->n_subs + 1, &pcbddc->adaptive_constrai… in PCBDDCAdaptiveSelection()
3458 &pcbddc->adaptive_constraints_data)); in PCBDDCAdaptiveSelection()
3459 PetscCall(PetscArrayzero(pcbddc->adaptive_constraints_n, nv + sub_schurs->n_subs)); in PCBDDCAdaptiveSelection()
3463 pcbddc->adaptive_constraints_idxs_ptr[0] = 0; in PCBDDCAdaptiveSelection()
3464 pcbddc->adaptive_constraints_data_ptr[0] = 0; in PCBDDCAdaptiveSelection()
3465 if (sub_schurs->is_vertices && pcbddc->use_vertices) { in PCBDDCAdaptiveSelection()
3468 PetscCall(ISGetIndices(sub_schurs->is_vertices, &idxs)); in PCBDDCAdaptiveSelection()
3470 pcbddc->adaptive_constraints_n[cum] = 1; in PCBDDCAdaptiveSelection()
3471 pcbddc->adaptive_constraints_idxs[cum] = idxs[cum]; in PCBDDCAdaptiveSelection()
3472 pcbddc->adaptive_constraints_data[cum] = 1.0; in PCBDDCAdaptiveSelection()
3473 … pcbddc->adaptive_constraints_idxs_ptr[cum + 1] = pcbddc->adaptive_constraints_idxs_ptr[cum] + 1; in PCBDDCAdaptiveSelection()
3474 … pcbddc->adaptive_constraints_data_ptr[cum + 1] = pcbddc->adaptive_constraints_data_ptr[cum] + 1; in PCBDDCAdaptiveSelection()
3476 PetscCall(ISRestoreIndices(sub_schurs->is_vertices, &idxs)); in PCBDDCAdaptiveSelection()
3480 if (sub_schurs->gdsw) { in PCBDDCAdaptiveSelection()
3481 PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_all, &Sarray)); in PCBDDCAdaptiveSelection()
3482 PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all, &Starray)); in PCBDDCAdaptiveSelection()
3484 PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_inv_all, &Sarray)); in PCBDDCAdaptiveSelection()
3485 PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all, &Starray)); in PCBDDCAdaptiveSelection()
3489 lthresh = pcbddc->adaptive_threshold[0]; in PCBDDCAdaptiveSelection()
3490 uthresh = pcbddc->adaptive_threshold[1]; in PCBDDCAdaptiveSelection()
3491 upart = pcbddc->use_deluxe_scaling; in PCBDDCAdaptiveSelection()
3492 for (i = 0; i < sub_schurs->n_subs; i++) { in PCBDDCAdaptiveSelection()
3504 if (sub_schurs->gdsw) { in PCBDDCAdaptiveSelection()
3508 …PetscCheck(sub_schurs->is_posdef, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented without del… in PCBDDCAdaptiveSelection()
3513 PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size)); in PCBDDCAdaptiveSelection()
3514 PetscCall(ISGetIndices(sub_schurs->is_subs[i], &idxs)); in PCBDDCAdaptiveSelection()
3516 /* this is experimental: we assume the dofs have been properly grouped to have in PCBDDCAdaptiveSelection()
3518 if (!sub_schurs->is_posdef) { in PCBDDCAdaptiveSelection()
3524 PetscCall(MatScale(T, -1.0)); in PCBDDCAdaptiveSelection()
3527 PetscCall(MatScale(T, -1.0)); in PCBDDCAdaptiveSelection()
3529 if (sub_schurs->change_primal_sub) { in PCBDDCAdaptiveSelection()
3533 PetscCall(ISGetLocalSize(sub_schurs->change_primal_sub[i], &nz)); in PCBDDCAdaptiveSelection()
3534 PetscCall(ISGetIndices(sub_schurs->change_primal_sub[i], &idxs)); in PCBDDCAdaptiveSelection()
3536 *(Sarray + cumarray + idxs[k] * (subset_size + 1)) *= -1.0; in PCBDDCAdaptiveSelection()
3539 PetscCall(ISRestoreIndices(sub_schurs->change_primal_sub[i], &idxs)); in PCBDDCAdaptiveSelection()
3548 if (sub_schurs->is_symmetric) { in PCBDDCAdaptiveSelection()
3550 if (sub_schurs->n_subs == 1) { /* zeroing memory to use PetscArraycmp() later */ in PCBDDCAdaptiveSelection()
3569 …if (sub_schurs->n_subs == 1 && pcbddc->use_deluxe_scaling) PetscCall(PetscArraycmp(S, St, subset_s… in PCBDDCAdaptiveSelection()
3571 if (same_data && !sub_schurs->change) { /* there's no need of constraints here */ in PCBDDCAdaptiveSelection()
3575 PetscReal eps = -1.0; /* dlamch? */ in PCBDDCAdaptiveSelection()
3579 PetscCheck(sub_schurs->is_symmetric, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented"); in PCBDDCAdaptiveSelection()
3584 if (pcbddc->dbg_flag) { in PCBDDCAdaptiveSelection()
3585 …PetscInt nc = 0, c = pcbddc->mat_graph->nodes[idxs[0]].count, w = pcbddc->mat_graph->nodes[idxs[0]… in PCBDDCAdaptiveSelection()
3587 …if (sub_schurs->change_primal_sub) PetscCall(ISGetLocalSize(sub_schurs->change_primal_sub[i], &nc)… in PCBDDCAdaptiveSelection()
3589 …->dbg_viewer, "Computing for sub %" PetscInt_FMT "/%" PetscInt_FMT " size %" PetscInt_FMT " count … in PCBDDCAdaptiveSelection()
3595 if (sub_schurs->is_posdef) { in PCBDDCAdaptiveSelection()
3606 …PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)pc)->prefix, "-pc_bddc_adaptive_recipe", &recipe,… in PCBDDCAdaptiveSelection()
3692 …PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)pc)->prefix, "-pc_bddc_adaptive_recipe3_min_scal"… in PCBDDCAdaptiveSelection()
3694 …PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)pc)->prefix, "-pc_bddc_adaptive_recipe3_min", &re… in PCBDDCAdaptiveSelection()
3706 if (recipe_m > 0 && B_N - B_neigs > 0) { in PCBDDCAdaptiveSelection()
3709 PetscCall(PetscBLASIntCast(PetscMin(recipe_m, B_N - B_neigs), &B_IU)); in PCBDDCAdaptiveSelection()
3780 …PetscCheck(sub_schurs->change_primal_sub, PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen… in PCBDDCAdaptiveSelection()
3781 PetscCall(ISGetLocalSize(sub_schurs->change_primal_sub[i], &nmax)); in PCBDDCAdaptiveSelection()
3792 …_ERR_LIB, "Error in SYGVX Lapack routine: illegal value for argument %" PetscBLASInt_FMT, -B_ierr); in PCBDDCAdaptiveSelection()
3794 …routine: leading minor of order %" PetscBLASInt_FMT " is not positive definite", B_ierr - B_N - 1); in PCBDDCAdaptiveSelection()
3798 …if (pcbddc->dbg_flag) PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " found %… in PCBDDCAdaptiveSelection()
3799 if (upart) eigs_start = scal ? 0 : B_neigs - nmax; in PCBDDCAdaptiveSelection()
3812 PetscCall(PetscBLASIntCast(B_N - nmin_s + 1, &B_IL)); in PCBDDCAdaptiveSelection()
3813 B_IU = B_N - B_neigs; in PCBDDCAdaptiveSelection()
3819 if (pcbddc->dbg_flag) { in PCBDDCAdaptiveSelection()
3820 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " found %" PetscBLASInt_FMT " e… in PCBDDCAdaptiveSelection()
3822 if (sub_schurs->is_symmetric) { in PCBDDCAdaptiveSelection()
3845 …_ERR_LIB, "Error in SYGVX Lapack routine: illegal value for argument %" PetscBLASInt_FMT, -B_ierr); in PCBDDCAdaptiveSelection()
3847 …routine: leading minor of order %" PetscBLASInt_FMT " is not positive definite", B_ierr - B_N - 1); in PCBDDCAdaptiveSelection()
3849 if (pcbddc->dbg_flag) { in PCBDDCAdaptiveSelection()
3850 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " -> Got %" PetscBLASInt_FMT " … in PCBDDCAdaptiveSelection()
3852 if (!sub_schurs->gdsw) { in PCBDDCAdaptiveSelection()
3854 PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " Inf\n")); in PCBDDCAdaptiveSelection()
3857 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " %1.6e\n", (double)eigs[j + … in PCBDDCAdaptiveSelection()
3859 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " %1.6e\n", (double)(1 / eigs… in PCBDDCAdaptiveSelection()
3865 PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " %1.6e\n", pg)); in PCBDDCAdaptiveSelection()
3871 if (sub_schurs->change) { in PCBDDCAdaptiveSelection()
3874 if (pcbddc->dbg_flag > 2) { in PCBDDCAdaptiveSelection()
3877 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " -> Eigenvector (old basis) %"… in PCBDDCAdaptiveSelection()
3882 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " %1.4e + %1.4e i\n", (doub… in PCBDDCAdaptiveSelection()
3884 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " %1.4e\n", (double)(eigv[(… in PCBDDCAdaptiveSelection()
3889 PetscCall(KSPGetOperators(sub_schurs->change[i], &change, NULL)); in PCBDDCAdaptiveSelection()
3897 pcbddc->adaptive_constraints_n[i + nv] = B_neigs; in PCBDDCAdaptiveSelection()
3899 …PetscCall(PetscArraycpy(pcbddc->adaptive_constraints_data + pcbddc->adaptive_constraints_data_ptr[… in PCBDDCAdaptiveSelection()
3901 if (pcbddc->dbg_flag > 1) { in PCBDDCAdaptiveSelection()
3904 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " -> Eigenvector %" PetscInt_FM… in PCBDDCAdaptiveSelection()
3907 …PetscReal r = PetscRealPart(pcbddc->adaptive_constraints_data[ii * subset_size + j + pcbddc->adapt… in PCBDDCAdaptiveSelection()
3908 …PetscReal c = PetscImaginaryPart(pcbddc->adaptive_constraints_data[ii * subset_size + j + pcbddc->… in PCBDDCAdaptiveSelection()
3909 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, " %1.4e + %1.4e i\n", (doub… in PCBDDCAdaptiveSelection()
3911 …izedPrintf(pcbddc->dbg_viewer, " %1.4e\n", (double)PetscRealPart(pcbddc->adaptive_constraint… in PCBDDCAdaptiveSelection()
3916 …PetscCall(PetscArraycpy(pcbddc->adaptive_constraints_idxs + pcbddc->adaptive_constraints_idxs_ptr[… in PCBDDCAdaptiveSelection()
3917 …pcbddc->adaptive_constraints_idxs_ptr[cum + 1] = pcbddc->adaptive_constraints_idxs_ptr[cum] + subs… in PCBDDCAdaptiveSelection()
3918 …pcbddc->adaptive_constraints_data_ptr[cum + 1] = pcbddc->adaptive_constraints_data_ptr[cum] + subs… in PCBDDCAdaptiveSelection()
3921 PetscCall(ISRestoreIndices(sub_schurs->is_subs[i], &idxs)); in PCBDDCAdaptiveSelection()
3925 if (pcbddc->dbg_flag) PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCAdaptiveSelection()
3928 if (sub_schurs->gdsw) { in PCBDDCAdaptiveSelection()
3929 PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_all, &Sarray)); in PCBDDCAdaptiveSelection()
3930 PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all, &Starray)); in PCBDDCAdaptiveSelection()
3932 PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_inv_all, &Sarray)); in PCBDDCAdaptiveSelection()
3933 PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all, &Starray)); in PCBDDCAdaptiveSelection()
3935 PetscCall(MatDestroy(&sub_schurs->sum_S_Ej_inv_all)); in PCBDDCAdaptiveSelection()
3936 PetscCall(MatDestroy(&sub_schurs->sum_S_Ej_tilda_all)); in PCBDDCAdaptiveSelection()
3944 if (pcbddc->dbg_flag) { in PCBDDCAdaptiveSelection()
3947 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "Maximum number of constraints per cc %" Pets… in PCBDDCAdaptiveSelection()
3949 PetscCall(PetscLogEventEnd(PC_BDDC_AdaptiveSetUp[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCAdaptiveSelection()
3980 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCResetCustomization()
3983 PetscCall(ISDestroy(&pcbddc->user_primal_vertices)); in PCBDDCResetCustomization()
3984 PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local)); in PCBDDCResetCustomization()
3985 PetscCall(ISDestroy(&pcbddc->NeumannBoundaries)); in PCBDDCResetCustomization()
3986 PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal)); in PCBDDCResetCustomization()
3987 PetscCall(ISDestroy(&pcbddc->DirichletBoundaries)); in PCBDDCResetCustomization()
3988 PetscCall(MatNullSpaceDestroy(&pcbddc->onearnullspace)); in PCBDDCResetCustomization()
3989 PetscCall(PetscFree(pcbddc->onearnullvecs_state)); in PCBDDCResetCustomization()
3990 PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal)); in PCBDDCResetCustomization()
3998 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCResetTopography()
4002 PetscCall(MatDestroy(&pcbddc->nedcG)); in PCBDDCResetTopography()
4003 PetscCall(ISDestroy(&pcbddc->nedclocal)); in PCBDDCResetTopography()
4004 PetscCall(MatDestroy(&pcbddc->discretegradient)); in PCBDDCResetTopography()
4005 PetscCall(MatDestroy(&pcbddc->user_ChangeOfBasisMatrix)); in PCBDDCResetTopography()
4006 PetscCall(MatDestroy(&pcbddc->ChangeOfBasisMatrix)); in PCBDDCResetTopography()
4007 PetscCall(MatDestroy(&pcbddc->switch_static_change)); in PCBDDCResetTopography()
4008 PetscCall(VecDestroy(&pcbddc->work_change)); in PCBDDCResetTopography()
4009 PetscCall(MatDestroy(&pcbddc->ConstraintMatrix)); in PCBDDCResetTopography()
4010 PetscCall(MatDestroy(&pcbddc->divudotp)); in PCBDDCResetTopography()
4011 PetscCall(ISDestroy(&pcbddc->divudotp_vl2l)); in PCBDDCResetTopography()
4012 PetscCall(PCBDDCGraphDestroy(&pcbddc->mat_graph)); in PCBDDCResetTopography()
4013 for (i = 0; i < pcbddc->n_local_subs; i++) PetscCall(ISDestroy(&pcbddc->local_subs[i])); in PCBDDCResetTopography()
4014 pcbddc->n_local_subs = 0; in PCBDDCResetTopography()
4015 PetscCall(PetscFree(pcbddc->local_subs)); in PCBDDCResetTopography()
4016 PetscCall(PCBDDCSubSchursDestroy(&pcbddc->sub_schurs)); in PCBDDCResetTopography()
4017 pcbddc->graphanalyzed = PETSC_FALSE; in PCBDDCResetTopography()
4018 pcbddc->recompute_topography = PETSC_TRUE; in PCBDDCResetTopography()
4019 pcbddc->corner_selected = PETSC_FALSE; in PCBDDCResetTopography()
4025 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCResetSolvers()
4028 PetscCall(VecDestroy(&pcbddc->coarse_vec)); in PCBDDCResetSolvers()
4029 PetscCall(MatDestroy(&pcbddc->coarse_phi_B)); in PCBDDCResetSolvers()
4030 PetscCall(MatDestroy(&pcbddc->coarse_phi_D)); in PCBDDCResetSolvers()
4031 PetscCall(MatDestroy(&pcbddc->coarse_psi_B)); in PCBDDCResetSolvers()
4032 PetscCall(MatDestroy(&pcbddc->coarse_psi_D)); in PCBDDCResetSolvers()
4033 PetscCall(VecDestroy(&pcbddc->vec1_P)); in PCBDDCResetSolvers()
4034 PetscCall(VecDestroy(&pcbddc->vec1_C)); in PCBDDCResetSolvers()
4035 PetscCall(MatDestroy(&pcbddc->local_auxmat2)); in PCBDDCResetSolvers()
4036 PetscCall(MatDestroy(&pcbddc->local_auxmat1)); in PCBDDCResetSolvers()
4037 PetscCall(VecDestroy(&pcbddc->vec1_R)); in PCBDDCResetSolvers()
4038 PetscCall(VecDestroy(&pcbddc->vec2_R)); in PCBDDCResetSolvers()
4039 PetscCall(ISDestroy(&pcbddc->is_R_local)); in PCBDDCResetSolvers()
4040 PetscCall(VecScatterDestroy(&pcbddc->R_to_B)); in PCBDDCResetSolvers()
4041 PetscCall(VecScatterDestroy(&pcbddc->R_to_D)); in PCBDDCResetSolvers()
4042 PetscCall(VecScatterDestroy(&pcbddc->coarse_loc_to_glob)); in PCBDDCResetSolvers()
4043 PetscCall(KSPReset(pcbddc->ksp_D)); in PCBDDCResetSolvers()
4044 PetscCall(KSPReset(pcbddc->ksp_R)); in PCBDDCResetSolvers()
4045 PetscCall(KSPReset(pcbddc->coarse_ksp)); in PCBDDCResetSolvers()
4046 PetscCall(MatDestroy(&pcbddc->local_mat)); in PCBDDCResetSolvers()
4047 PetscCall(PetscFree(pcbddc->primal_indices_local_idxs)); in PCBDDCResetSolvers()
4048 PetscCall(PetscFree2(pcbddc->local_primal_ref_node, pcbddc->local_primal_ref_mult)); in PCBDDCResetSolvers()
4049 PetscCall(PetscFree(pcbddc->global_primal_indices)); in PCBDDCResetSolvers()
4050 PetscCall(ISDestroy(&pcbddc->coarse_subassembling)); in PCBDDCResetSolvers()
4051 PetscCall(MatDestroy(&pcbddc->benign_change)); in PCBDDCResetSolvers()
4052 PetscCall(VecDestroy(&pcbddc->benign_vec)); in PCBDDCResetSolvers()
4054 PetscCall(MatDestroy(&pcbddc->benign_B0)); in PCBDDCResetSolvers()
4055 PetscCall(PetscSFDestroy(&pcbddc->benign_sf)); in PCBDDCResetSolvers()
4056 if (pcbddc->benign_zerodiag_subs) { in PCBDDCResetSolvers()
4058 for (i = 0; i < pcbddc->benign_n; i++) PetscCall(ISDestroy(&pcbddc->benign_zerodiag_subs[i])); in PCBDDCResetSolvers()
4059 PetscCall(PetscFree(pcbddc->benign_zerodiag_subs)); in PCBDDCResetSolvers()
4061 PetscCall(PetscFree3(pcbddc->benign_p0_lidx, pcbddc->benign_p0_gidx, pcbddc->benign_p0)); in PCBDDCResetSolvers()
4067 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCSetUpLocalWorkVectors()
4068 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCSetUpLocalWorkVectors()
4073 n_constraints = pcbddc->local_primal_size - pcbddc->benign_n - pcbddc->n_vertices; in PCBDDCSetUpLocalWorkVectors()
4074 n_R = pcis->n - pcbddc->n_vertices; in PCBDDCSetUpLocalWorkVectors()
4075 PetscCall(VecGetType(pcis->vec1_N, &impVecType)); in PCBDDCSetUpLocalWorkVectors()
4078 old_size = -1; in PCBDDCSetUpLocalWorkVectors()
4079 if (pcbddc->vec1_R) PetscCall(VecGetSize(pcbddc->vec1_R, &old_size)); in PCBDDCSetUpLocalWorkVectors()
4081 PetscCall(VecDestroy(&pcbddc->vec1_R)); in PCBDDCSetUpLocalWorkVectors()
4082 PetscCall(VecDestroy(&pcbddc->vec2_R)); in PCBDDCSetUpLocalWorkVectors()
4083 PetscCall(VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N), &pcbddc->vec1_R)); in PCBDDCSetUpLocalWorkVectors()
4084 PetscCall(VecSetSizes(pcbddc->vec1_R, PETSC_DECIDE, n_R)); in PCBDDCSetUpLocalWorkVectors()
4085 PetscCall(VecSetType(pcbddc->vec1_R, impVecType)); in PCBDDCSetUpLocalWorkVectors()
4086 PetscCall(VecDuplicate(pcbddc->vec1_R, &pcbddc->vec2_R)); in PCBDDCSetUpLocalWorkVectors()
4089 old_size = -1; in PCBDDCSetUpLocalWorkVectors()
4090 if (pcbddc->vec1_P) PetscCall(VecGetSize(pcbddc->vec1_P, &old_size)); in PCBDDCSetUpLocalWorkVectors()
4091 if (pcbddc->local_primal_size != old_size) { in PCBDDCSetUpLocalWorkVectors()
4092 PetscCall(VecDestroy(&pcbddc->vec1_P)); in PCBDDCSetUpLocalWorkVectors()
4093 PetscCall(VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N), &pcbddc->vec1_P)); in PCBDDCSetUpLocalWorkVectors()
4094 PetscCall(VecSetSizes(pcbddc->vec1_P, PETSC_DECIDE, pcbddc->local_primal_size)); in PCBDDCSetUpLocalWorkVectors()
4095 PetscCall(VecSetType(pcbddc->vec1_P, impVecType)); in PCBDDCSetUpLocalWorkVectors()
4098 old_size = -1; in PCBDDCSetUpLocalWorkVectors()
4099 if (pcbddc->vec1_C) PetscCall(VecGetSize(pcbddc->vec1_C, &old_size)); in PCBDDCSetUpLocalWorkVectors()
4101 PetscCall(VecDestroy(&pcbddc->vec1_C)); in PCBDDCSetUpLocalWorkVectors()
4102 PetscCall(VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N), &pcbddc->vec1_C)); in PCBDDCSetUpLocalWorkVectors()
4103 PetscCall(VecSetSizes(pcbddc->vec1_C, PETSC_DECIDE, n_constraints)); in PCBDDCSetUpLocalWorkVectors()
4104 PetscCall(VecSetType(pcbddc->vec1_C, impVecType)); in PCBDDCSetUpLocalWorkVectors()
4143 const PetscInt nci = ii[i + 1] - ii[i]; in MatSetValuesSubMat()
4190 aij = (Mat_SeqAIJ *)(*mat)->data; in MatCreateSeqAIJFromDenseExpand()
4191 aij->free_a = PETSC_TRUE; in MatCreateSeqAIJFromDenseExpand()
4192 aij->free_ij = PETSC_TRUE; in MatCreateSeqAIJFromDenseExpand()
4199 PetscInt n = A->rmap->n, ncnt = 0, ncnt2 = 0, bsizemax = 0, *v_pivots = NULL; in MatSeqAIJInvertVariableBlockDiagonalMat()
4208 …PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Not for rectangular matr… in MatSeqAIJInvertVariableBlockDiagonalMat()
4269 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)(*B)->data; in MatSeqAIJInvertVariableBlockDiagonalMat()
4270 aij->free_a = PETSC_TRUE; in MatSeqAIJInvertVariableBlockDiagonalMat()
4271 aij->free_ij = PETSC_TRUE; in MatSeqAIJInvertVariableBlockDiagonalMat()
4278 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCSetUpCorrection()
4279 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCSetUpCorrection()
4280 PCBDDCGraph graph = pcbddc->mat_graph; in PCBDDCSetUpCorrection()
4281 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCSetUpCorrection()
4300 PetscScalar one = 1.0, m_one = -1.0; in PCBDDCSetUpCorrection()
4302 /* Multi-element support */ in PCBDDCSetUpCorrection()
4303 PetscBool multi_element = graph->multi_element; in PCBDDCSetUpCorrection()
4311 …PetscCheck(pcbddc->symmetric_primal || !pcbddc->benign_n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-sym… in PCBDDCSetUpCorrection()
4312 PetscCall(PetscLogEventBegin(PC_BDDC_CorrectionSetUp[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpCorrection()
4314 /* Set Non-overlapping dimensions */ in PCBDDCSetUpCorrection()
4315 n_vertices = pcbddc->n_vertices; in PCBDDCSetUpCorrection()
4316 n_constraints = pcbddc->local_primal_size - pcbddc->benign_n - n_vertices; in PCBDDCSetUpCorrection()
4317 n_B = pcis->n_B; in PCBDDCSetUpCorrection()
4318 n_D = pcis->n - n_B; in PCBDDCSetUpCorrection()
4319 n_R = pcis->n - n_vertices; in PCBDDCSetUpCorrection()
4323 …PetscCall(ISGlobalToLocalMappingApply(pcis->BtoNmap, IS_GTOLM_DROP, n_vertices, pcbddc->local_prim… in PCBDDCSetUpCorrection()
4327 if (pcbddc->benign_saddle_point || !pcbddc->symmetric_primal) multi_element = PETSC_FALSE; in PCBDDCSetUpCorrection()
4329 /* Subdomain contribution (Non-overlapping) to coarse matrix */ in PCBDDCSetUpCorrection()
4331 PetscCheck(!pcbddc->benign_n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented"); in PCBDDCSetUpCorrection()
4334 …(MatSetSizes(*coarse_submat, pcbddc->local_primal_size, pcbddc->local_primal_size, pcbddc->local_p… in PCBDDCSetUpCorrection()
4340 const PetscInt *vidxs = pcbddc->primal_indices_local_idxs; in PCBDDCSetUpCorrection()
4341 const PetscInt *cidxs = pcbddc->primal_indices_local_idxs + n_vertices; in PCBDDCSetUpCorrection()
4343 PetscInt n_el = PetscMax(graph->n_local_subs, 1); in PCBDDCSetUpCorrection()
4349 PetscInt s = 2 * graph->nodes[vidxs[i]].local_sub; in PCBDDCSetUpCorrection()
4355 PetscInt s = 2 * graph->nodes[cidxs[i]].local_sub + 1; in PCBDDCSetUpCorrection()
4364 PetscInt s = 2 * graph->nodes[vidxs[i]].local_sub; in PCBDDCSetUpCorrection()
4369 PetscInt s = 2 * graph->nodes[cidxs[i]].local_sub; in PCBDDCSetUpCorrection()
4389 const PetscInt e = graph->nodes[vidxs[i]].local_sub; in PCBDDCSetUpCorrection()
4396 const PetscInt e = graph->nodes[cidxs[i]].local_sub; in PCBDDCSetUpCorrection()
4407 for (PetscInt i = 0; i < n_R * n_eff_vertices; i++) R_eff_V_J[i] = -1; in PCBDDCSetUpCorrection()
4408 for (PetscInt i = 0; i < n_R * n_eff_constraints; i++) R_eff_C_J[i] = -1; in PCBDDCSetUpCorrection()
4409 for (PetscInt i = 0; i < n_B * n_eff_vertices; i++) B_eff_V_J[i] = -1; in PCBDDCSetUpCorrection()
4410 for (PetscInt i = 0; i < n_B * n_eff_constraints; i++) B_eff_C_J[i] = -1; in PCBDDCSetUpCorrection()
4412 PetscCall(ISGetIndices(pcbddc->is_R_local, &idx)); in PCBDDCSetUpCorrection()
4414 const PetscInt e = graph->nodes[idx[i]].local_sub; in PCBDDCSetUpCorrection()
4421 PetscCall(ISRestoreIndices(pcbddc->is_R_local, &idx)); in PCBDDCSetUpCorrection()
4422 PetscCall(ISGetIndices(pcis->is_B_local, &idx)); in PCBDDCSetUpCorrection()
4424 const PetscInt e = graph->nodes[idx[i]].local_sub; in PCBDDCSetUpCorrection()
4431 PetscCall(ISRestoreIndices(pcis->is_B_local, &idx)); in PCBDDCSetUpCorrection()
4452 …PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pcbddc->local_primal_size, pcbddc->local_primal_size,… in PCBDDCSetUpCorrection()
4458 PetscCall(KSPGetPC(pcbddc->ksp_R, &pc_R)); in PCBDDCSetUpCorrection()
4460 PetscCall(PetscObjectTypeCompare((PetscObject)pcbddc->ksp_R, KSPPREONLY, &isPreonly)); in PCBDDCSetUpCorrection()
4468 } else if (sub_schurs && sub_schurs->reuse_solver) { in PCBDDCSetUpCorrection()
4469 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpCorrection()
4472 F = reuse_solver->F; in PCBDDCSetUpCorrection()
4477 need_benign_correction = (PetscBool)(!!reuse_solver->benign_n); in PCBDDCSetUpCorrection()
4480 /* determine if we can use a sparse right-hand side */ in PCBDDCSetUpCorrection()
4492 PetscCall(VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N), &dummy_vec)); in PCBDDCSetUpCorrection()
4494 PetscCall(VecSetType(dummy_vec, ((PetscObject)pcis->vec1_N)->type_name)); in PCBDDCSetUpCorrection()
4497 PetscCall(MatDestroy(&pcbddc->local_auxmat1)); in PCBDDCSetUpCorrection()
4498 PetscCall(MatDestroy(&pcbddc->local_auxmat2)); in PCBDDCSetUpCorrection()
4511 …PetscCall(MatCreateSubMatrix(pcbddc->ConstraintMatrix, is_C, pcbddc->is_R_local, MAT_INITIAL_MATRI… in PCBDDCSetUpCorrection()
4512 …PetscCall(MatCreateSubMatrix(pcbddc->ConstraintMatrix, is_C, pcis->is_B_local, MAT_INITIAL_MATRIX,… in PCBDDCSetUpCorrection()
4514 /* Assemble local_auxmat2_R = (- A_{RR}^{-1} C^T_{CR}) needed by BDDC setup */ in PCBDDCSetUpCorrection()
4515 …/* Assemble pcbddc->local_auxmat2 = R_to_B (- A_{RR}^{-1} C^T_{CR}) needed by BDDC application */ in PCBDDCSetUpCorrection()
4527 …for (j = 0; j < size_of_constraint; j++) marr[row_cmat_indices[j] + col * lda_rhs] = -row_cmat_val… in PCBDDCSetUpCorrection()
4534 PetscCall(MatScale(C_CR, -1.0)); in PCBDDCSetUpCorrection()
4556 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpCorrection()
4559 PetscCall(PetscArrayzero(reuse_solver->benign_save_vals, pcbddc->benign_n)); in PCBDDCSetUpCorrection()
4564 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpCorrection()
4576 PetscCall(VecPlaceArray(pcbddc->vec1_R, marr + i * lda_rhs)); in PCBDDCSetUpCorrection()
4577 …PetscCall(PCBDDCReuseSolversBenignAdapt(reuse_solver, pcbddc->vec1_R, NULL, PETSC_TRUE, PETSC_TRUE… in PCBDDCSetUpCorrection()
4578 PetscCall(VecResetArray(pcbddc->vec1_R)); in PCBDDCSetUpCorrection()
4590 PetscCall(VecPlaceArray(pcbddc->vec1_R, barr + i * lda_rhs)); in PCBDDCSetUpCorrection()
4591 PetscCall(VecPlaceArray(pcbddc->vec2_R, marr + i * lda_rhs)); in PCBDDCSetUpCorrection()
4592 PetscCall(KSPSolve(pcbddc->ksp_R, pcbddc->vec1_R, pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
4593 PetscCall(KSPCheckSolve(pcbddc->ksp_R, pc, pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
4594 PetscCall(VecResetArray(pcbddc->vec1_R)); in PCBDDCSetUpCorrection()
4595 PetscCall(VecResetArray(pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
4600 if (sparserhs) PetscCall(MatScale(C_CR, -1.0)); in PCBDDCSetUpCorrection()
4602 /* Assemble explicitly S_CC = ( C_{CR} A_{RR}^{-1} C^T_{CR})^{-1} */ in PCBDDCSetUpCorrection()
4603 if (!pcbddc->switch_static) { in PCBDDCSetUpCorrection()
4604 …PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n_B, n_eff_constraints, NULL, &pcbddc->local_auxmat2)… in PCBDDCSetUpCorrection()
4605 …PetscCall(MatDenseScatter_Private(pcbddc->R_to_B, local_auxmat2_R, pcbddc->local_auxmat2, INSERT_V… in PCBDDCSetUpCorrection()
4612 … PetscCall(MatCreateSeqAIJFromDenseExpand(pcbddc->local_auxmat2, n_constraints, B_eff_C_J, &T)); in PCBDDCSetUpCorrection()
4613 PetscCall(MatDestroy(&pcbddc->local_auxmat2)); in PCBDDCSetUpCorrection()
4614 pcbddc->local_auxmat2 = T; in PCBDDCSetUpCorrection()
4616 PetscCall(MatMatMult(C_B, pcbddc->local_auxmat2, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &S_CC)); in PCBDDCSetUpCorrection()
4626 …PetscCall(MatCreateSubMatrix(local_auxmat2_R, is_R, NULL, MAT_INITIAL_MATRIX, &pcbddc->local_auxma… in PCBDDCSetUpCorrection()
4629 pcbddc->local_auxmat2 = local_auxmat2_R; in PCBDDCSetUpCorrection()
4631 … PetscCall(MatMatMult(C_CR, pcbddc->local_auxmat2, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &S_CC)); in PCBDDCSetUpCorrection()
4657 PetscCall(MatMatMult(S_CC, C_B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &pcbddc->local_auxmat1)); in PCBDDCSetUpCorrection()
4669 …if (sub_schurs && sub_schurs->reuse_solver) { /* is_R_local is not sorted, ISComplement doesn't li… in PCBDDCSetUpCorrection()
4672 PetscCall(ISDuplicate(pcbddc->is_R_local, &tis)); in PCBDDCSetUpCorrection()
4674 PetscCall(ISComplement(tis, 0, pcis->n, &is_aux)); in PCBDDCSetUpCorrection()
4677 PetscCall(ISComplement(pcbddc->is_R_local, 0, pcis->n, &is_aux)); in PCBDDCSetUpCorrection()
4680 oldpin = pcbddc->local_mat->boundtocpu; in PCBDDCSetUpCorrection()
4682 PetscCall(MatBindToCPU(pcbddc->local_mat, PETSC_TRUE)); in PCBDDCSetUpCorrection()
4683 …PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcbddc->is_R_local, is_aux, MAT_INITIAL_MATRIX, &A… in PCBDDCSetUpCorrection()
4684 …PetscCall(MatCreateSubMatrix(pcbddc->local_mat, is_aux, pcbddc->is_R_local, MAT_INITIAL_MATRIX, &A… in PCBDDCSetUpCorrection()
4687 PetscCall(MatCreateSubMatrix(pcbddc->local_mat, is_aux, is_aux, MAT_INITIAL_MATRIX, &A_VV)); in PCBDDCSetUpCorrection()
4689 PetscCall(MatBindToCPU(pcbddc->local_mat, oldpin)); in PCBDDCSetUpCorrection()
4697 if (pcbddc->benign_n && (pcbddc->switch_static || pcbddc->dbg_flag)) { in PCBDDCSetUpCorrection()
4700 PetscCall(ISGetIndices(pcis->is_I_local, &idxs)); in PCBDDCSetUpCorrection()
4701 PetscCall(PetscMalloc1(pcbddc->benign_n, &p0_lidx_I)); in PCBDDCSetUpCorrection()
4702 …for (i = 0; i < pcbddc->benign_n; i++) PetscCall(PetscFindInt(pcbddc->benign_p0_lidx[i], pcis->n -… in PCBDDCSetUpCorrection()
4703 PetscCall(ISRestoreIndices(pcis->is_I_local, &idxs)); in PCBDDCSetUpCorrection()
4709 PetscCall(MatDestroy(&pcbddc->coarse_phi_B)); in PCBDDCSetUpCorrection()
4710 PetscCall(MatDestroy(&pcbddc->coarse_psi_B)); in PCBDDCSetUpCorrection()
4711 PetscCall(MatDestroy(&pcbddc->coarse_phi_D)); in PCBDDCSetUpCorrection()
4712 PetscCall(MatDestroy(&pcbddc->coarse_psi_D)); in PCBDDCSetUpCorrection()
4714 …PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n_B, pcbddc->local_primal_size, NULL, &pcbddc->coarse… in PCBDDCSetUpCorrection()
4715 …if (pcbddc->switch_static || pcbddc->dbg_flag) PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n_D, p… in PCBDDCSetUpCorrection()
4718 IS is_rows[2] = {pcbddc->is_R_local, NULL}; in PCBDDCSetUpCorrection()
4721 …PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n_vertices, pcbddc->local_primal_ref_node, PETSC_USE_PO… in PCBDDCSetUpCorrection()
4743 PetscCall(PetscMalloc1(pcbddc->benign_n, &idxs_p0)); in PCBDDCSetUpCorrection()
4744 PetscCall(ISLocalToGlobalMappingCreateIS(pcbddc->is_R_local, &RtoN)); in PCBDDCSetUpCorrection()
4745 …PetscCall(ISGlobalToLocalMappingApply(RtoN, IS_GTOLM_DROP, pcbddc->benign_n, pcbddc->benign_p0_lid… in PCBDDCSetUpCorrection()
4746 …n == pcbddc->benign_n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in R numbering for benign p0! %" Pe… in PCBDDCSetUpCorrection()
4779 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpCorrection()
4787 L = | 0 0 -1 | (P-p0) in PCBDDCSetUpCorrection()
4788 | 0 0 -1 | (p0) in PCBDDCSetUpCorrection()
4791 for (i = 0; i < reuse_solver->benign_n; i++) { in PCBDDCSetUpCorrection()
4796 PetscCall(ISGetLocalSize(reuse_solver->benign_zerodiag_subs[i], &nz)); in PCBDDCSetUpCorrection()
4797 PetscCall(ISGetIndices(reuse_solver->benign_zerodiag_subs[i], &idxs_zero)); in PCBDDCSetUpCorrection()
4802 for (k = 0; k < nz; k++) marr[idxs_zero[k] + lda_rhs * col] -= val; in PCBDDCSetUpCorrection()
4805 PetscCall(ISRestoreIndices(reuse_solver->benign_zerodiag_subs[i], &idxs_zero)); in PCBDDCSetUpCorrection()
4814 if (!pcbddc->symmetric_primal) { in PCBDDCSetUpCorrection()
4815 /* A_RV already scaled by -1 */ in PCBDDCSetUpCorrection()
4819 PetscCall(MatScale(A_VR, -1.0)); in PCBDDCSetUpCorrection()
4845 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpCorrection()
4857 PetscCall(VecPlaceArray(pcbddc->vec1_R, marr + i * lda_rhs)); in PCBDDCSetUpCorrection()
4858 …PetscCall(PCBDDCReuseSolversBenignAdapt(reuse_solver, pcbddc->vec1_R, NULL, PETSC_FALSE, PETSC_TRU… in PCBDDCSetUpCorrection()
4859 PetscCall(VecResetArray(pcbddc->vec1_R)); in PCBDDCSetUpCorrection()
4865 if (restoreavr) PetscCall(MatScale(A_VR, -1.0)); in PCBDDCSetUpCorrection()
4868 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpCorrection()
4880 PetscCall(VecPlaceArray(pcbddc->vec1_R, marr + i * lda_rhs)); in PCBDDCSetUpCorrection()
4881 …PetscCall(PCBDDCReuseSolversBenignAdapt(reuse_solver, pcbddc->vec1_R, NULL, PETSC_TRUE, PETSC_TRUE… in PCBDDCSetUpCorrection()
4882 PetscCall(VecResetArray(pcbddc->vec1_R)); in PCBDDCSetUpCorrection()
4894 PetscCall(VecPlaceArray(pcbddc->vec1_R, barr + i * lda_rhs)); in PCBDDCSetUpCorrection()
4895 PetscCall(VecPlaceArray(pcbddc->vec2_R, marr + i * lda_rhs)); in PCBDDCSetUpCorrection()
4896 PetscCall(KSPSolve(pcbddc->ksp_R, pcbddc->vec1_R, pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
4897 PetscCall(KSPCheckSolve(pcbddc->ksp_R, pc, pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
4898 PetscCall(VecResetArray(pcbddc->vec1_R)); in PCBDDCSetUpCorrection()
4899 PetscCall(VecResetArray(pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
4911 … PetscCall(MatDenseScatter_Private(pcbddc->R_to_B, A_RRmA_RV, B, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
4913 /* S_CV = pcbddc->local_auxmat1 * B */ in PCBDDCSetUpCorrection()
4921 PetscCall(MatProductCreate(pcbddc->local_auxmat1, B, NULL, &S_CV)); in PCBDDCSetUpCorrection()
4966 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpCorrection()
4972 for (i = 0; i < reuse_solver->benign_n; i++) { in PCBDDCSetUpCorrection()
4977 PetscCall(ISGetLocalSize(reuse_solver->benign_zerodiag_subs[i], &nz)); in PCBDDCSetUpCorrection()
4978 PetscCall(ISGetIndices(reuse_solver->benign_zerodiag_subs[i], &idxs_zero)); in PCBDDCSetUpCorrection()
4989 PetscCall(ISRestoreIndices(reuse_solver->benign_zerodiag_subs[i], &idxs_zero)); in PCBDDCSetUpCorrection()
5014 …PetscCall(MatDenseGetSubMatrix(pcbddc->coarse_phi_B, PETSC_DECIDE, PETSC_DECIDE, 0, n_vertices, &B… in PCBDDCSetUpCorrection()
5015 … PetscCall(MatDenseScatter_Private(pcbddc->R_to_B, A_RRmA_RV, B, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
5016 PetscCall(MatDenseRestoreSubMatrix(pcbddc->coarse_phi_B, &B)); in PCBDDCSetUpCorrection()
5017 if (pcbddc->switch_static || pcbddc->dbg_flag) { in PCBDDCSetUpCorrection()
5018 …PetscCall(MatDenseGetSubMatrix(pcbddc->coarse_phi_D, PETSC_DECIDE, PETSC_DECIDE, 0, n_vertices, &B… in PCBDDCSetUpCorrection()
5019 … PetscCall(MatDenseScatter_Private(pcbddc->R_to_D, A_RRmA_RV, B, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
5020 PetscCall(MatDenseRestoreSubMatrix(pcbddc->coarse_phi_D, &B)); in PCBDDCSetUpCorrection()
5021 if (pcbddc->benign_n) { in PCBDDCSetUpCorrection()
5022 …for (i = 0; i < n_vertices; i++) PetscCall(MatSetValues(pcbddc->coarse_phi_D, pcbddc->benign_n, p0… in PCBDDCSetUpCorrection()
5023 PetscCall(MatAssemblyBegin(pcbddc->coarse_phi_D, MAT_FINAL_ASSEMBLY)); in PCBDDCSetUpCorrection()
5024 PetscCall(MatAssemblyEnd(pcbddc->coarse_phi_D, MAT_FINAL_ASSEMBLY)); in PCBDDCSetUpCorrection()
5028 …for (i = 0; i < n_vertices; i++) PetscCall(MatSetValues(pcbddc->coarse_phi_B, 1, &idx_V_B[i], 1, &… in PCBDDCSetUpCorrection()
5029 PetscCall(MatAssemblyBegin(pcbddc->coarse_phi_B, MAT_FINAL_ASSEMBLY)); in PCBDDCSetUpCorrection()
5030 PetscCall(MatAssemblyEnd(pcbddc->coarse_phi_B, MAT_FINAL_ASSEMBLY)); in PCBDDCSetUpCorrection()
5067 …PetscCall(MatDenseGetSubMatrix(pcbddc->coarse_phi_B, PETSC_DECIDE, PETSC_DECIDE, n_vertices, n_ver… in PCBDDCSetUpCorrection()
5068 PetscCall(MatDenseScatter_Private(pcbddc->R_to_B, B, B2, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
5069 PetscCall(MatDenseRestoreSubMatrix(pcbddc->coarse_phi_B, &B2)); in PCBDDCSetUpCorrection()
5070 if (pcbddc->switch_static || pcbddc->dbg_flag) { in PCBDDCSetUpCorrection()
5071 …PetscCall(MatDenseGetSubMatrix(pcbddc->coarse_phi_D, PETSC_DECIDE, PETSC_DECIDE, n_vertices, n_ver… in PCBDDCSetUpCorrection()
5072 PetscCall(MatDenseScatter_Private(pcbddc->R_to_D, B, B2, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
5073 if (pcbddc->benign_n) { in PCBDDCSetUpCorrection()
5074 …for (i = 0; i < n_constraints; i++) PetscCall(MatSetValues(B2, pcbddc->benign_n, p0_lidx_I, 1, &i,… in PCBDDCSetUpCorrection()
5076 PetscCall(MatDenseRestoreSubMatrix(pcbddc->coarse_phi_D, &B2)); in PCBDDCSetUpCorrection()
5088 …PetscCall(MatCreateSubMatrix(T, pcis->is_B_local, NULL, MAT_INITIAL_MATRIX, &pcbddc->coarse_phi_B)… in PCBDDCSetUpCorrection()
5089 …if (pcbddc->switch_static || pcbddc->dbg_flag) PetscCall(MatCreateSubMatrix(T, pcis->is_I_local, N… in PCBDDCSetUpCorrection()
5096 if (pcbddc->benign_n) { in PCBDDCSetUpCorrection()
5102 PetscCall(ISCreateStride(PETSC_COMM_SELF, pcbddc->benign_n, 0, 1, &is_dummy)); in PCBDDCSetUpCorrection()
5103 …PetscCall(MatCreateSubMatrix(pcbddc->benign_B0, is_dummy, pcis->is_B_local, MAT_INITIAL_MATRIX, &B… in PCBDDCSetUpCorrection()
5105 PetscCall(MatMatMult(B0_B, pcbddc->coarse_phi_B, MAT_INITIAL_MATRIX, 1.0, &B0_BPHI)); in PCBDDCSetUpCorrection()
5108 for (j = 0; j < pcbddc->benign_n; j++) { in PCBDDCSetUpCorrection()
5109 PetscInt primal_idx = pcbddc->local_primal_size - pcbddc->benign_n + j; in PCBDDCSetUpCorrection()
5110 for (i = 0; i < pcbddc->local_primal_size; i++) { in PCBDDCSetUpCorrection()
5111 …PetscCall(MatSetValue(*coarse_submat, primal_idx, i, data[i * pcbddc->benign_n + j], INSERT_VALUES… in PCBDDCSetUpCorrection()
5112 …PetscCall(MatSetValue(*coarse_submat, i, primal_idx, data[i * pcbddc->benign_n + j], INSERT_VALUES… in PCBDDCSetUpCorrection()
5120 /* compute other basis functions for non-symmetric problems */ in PCBDDCSetUpCorrection()
5121 if (!pcbddc->symmetric_primal) { in PCBDDCSetUpCorrection()
5153 /* B_V = B_V - A_VR^T */ in PCBDDCSetUpCorrection()
5160 for (j = xadj[i]; j < xadj[i + 1]; j++) marray[i * n_R + adjncy[j]] -= av[j]; in PCBDDCSetUpCorrection()
5168 PetscCall(PetscMalloc1(n_R * pcbddc->local_primal_size, &work)); in PCBDDCSetUpCorrection()
5172 PetscCall(VecPlaceArray(pcbddc->vec1_R, marray + i * n_R)); in PCBDDCSetUpCorrection()
5173 PetscCall(VecPlaceArray(pcbddc->vec2_R, work + i * n_R)); in PCBDDCSetUpCorrection()
5174 PetscCall(KSPSolveTranspose(pcbddc->ksp_R, pcbddc->vec1_R, pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
5175 PetscCall(KSPCheckSolve(pcbddc->ksp_R, pc, pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
5176 PetscCall(VecResetArray(pcbddc->vec1_R)); in PCBDDCSetUpCorrection()
5177 PetscCall(VecResetArray(pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
5184 PetscCall(VecPlaceArray(pcbddc->vec1_R, marray + (i - n_vertices) * n_R)); in PCBDDCSetUpCorrection()
5185 PetscCall(VecPlaceArray(pcbddc->vec2_R, work + i * n_R)); in PCBDDCSetUpCorrection()
5186 PetscCall(KSPSolveTranspose(pcbddc->ksp_R, pcbddc->vec1_R, pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
5187 PetscCall(KSPCheckSolve(pcbddc->ksp_R, pc, pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
5188 PetscCall(VecResetArray(pcbddc->vec1_R)); in PCBDDCSetUpCorrection()
5189 PetscCall(VecResetArray(pcbddc->vec2_R)); in PCBDDCSetUpCorrection()
5194 …PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n_B, pcbddc->local_primal_size, NULL, &pcbddc->coarse… in PCBDDCSetUpCorrection()
5195 …if (pcbddc->switch_static || pcbddc->dbg_flag) PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n_D, p… in PCBDDCSetUpCorrection()
5196 for (i = 0; i < pcbddc->local_primal_size; i++) { in PCBDDCSetUpCorrection()
5199 PetscCall(VecPlaceArray(pcbddc->vec1_R, work + i * n_R)); in PCBDDCSetUpCorrection()
5200 PetscCall(MatDenseGetColumnVec(pcbddc->coarse_psi_B, i, &v)); in PCBDDCSetUpCorrection()
5201 PetscCall(VecScatterBegin(pcbddc->R_to_B, pcbddc->vec1_R, v, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
5202 PetscCall(VecScatterEnd(pcbddc->R_to_B, pcbddc->vec1_R, v, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
5209 PetscCall(MatDenseRestoreColumnVec(pcbddc->coarse_psi_B, i, &v)); in PCBDDCSetUpCorrection()
5211 if (pcbddc->switch_static || pcbddc->dbg_flag) { in PCBDDCSetUpCorrection()
5212 PetscCall(MatDenseGetColumnVec(pcbddc->coarse_psi_D, i, &v)); in PCBDDCSetUpCorrection()
5213 … PetscCall(VecScatterBegin(pcbddc->R_to_D, pcbddc->vec1_R, v, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
5214 PetscCall(VecScatterEnd(pcbddc->R_to_D, pcbddc->vec1_R, v, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSetUpCorrection()
5215 PetscCall(MatDenseRestoreColumnVec(pcbddc->coarse_psi_D, i, &v)); in PCBDDCSetUpCorrection()
5217 PetscCall(VecResetArray(pcbddc->vec1_R)); in PCBDDCSetUpCorrection()
5223 PetscCall(PetscObjectReference((PetscObject)pcbddc->coarse_phi_B)); in PCBDDCSetUpCorrection()
5224 pcbddc->coarse_psi_B = pcbddc->coarse_phi_B; in PCBDDCSetUpCorrection()
5225 PetscCall(PetscObjectReference((PetscObject)pcbddc->coarse_phi_D)); in PCBDDCSetUpCorrection()
5226 pcbddc->coarse_psi_D = pcbddc->coarse_phi_D; in PCBDDCSetUpCorrection()
5249 PetscCall(PetscLogEventEnd(PC_BDDC_CorrectionSetUp[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpCorrection()
5253 /* Non-symmetric case : It should be \Psi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */ in PCBDDCSetUpCorrection()
5254 if (pcbddc->dbg_flag) { in PCBDDCSetUpCorrection()
5265 if (pcbddc->benign_n && !pcbddc->benign_change_explicit) { in PCBDDCSetUpCorrection()
5268 … PetscCall(MatCreateSubMatrix(A, pcis->is_I_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &A_II)); in PCBDDCSetUpCorrection()
5269 … PetscCall(MatCreateSubMatrix(A, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &A_IB)); in PCBDDCSetUpCorrection()
5270 … PetscCall(MatCreateSubMatrix(A, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &A_BI)); in PCBDDCSetUpCorrection()
5271 … PetscCall(MatCreateSubMatrix(A, pcis->is_B_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &A_BB)); in PCBDDCSetUpCorrection()
5274 PetscCall(MatConvert(pcis->A_II, checkmattype, MAT_INITIAL_MATRIX, &A_II)); in PCBDDCSetUpCorrection()
5275 PetscCall(MatConvert(pcis->A_IB, checkmattype, MAT_INITIAL_MATRIX, &A_IB)); in PCBDDCSetUpCorrection()
5276 PetscCall(MatConvert(pcis->A_BI, checkmattype, MAT_INITIAL_MATRIX, &A_BI)); in PCBDDCSetUpCorrection()
5277 PetscCall(MatConvert(pcis->A_BB, checkmattype, MAT_INITIAL_MATRIX, &A_BB)); in PCBDDCSetUpCorrection()
5279 PetscCall(MatConvert(pcbddc->coarse_phi_D, checkmattype, MAT_INITIAL_MATRIX, &coarse_phi_D)); in PCBDDCSetUpCorrection()
5280 PetscCall(MatConvert(pcbddc->coarse_phi_B, checkmattype, MAT_INITIAL_MATRIX, &coarse_phi_B)); in PCBDDCSetUpCorrection()
5281 if (!pcbddc->symmetric_primal) { in PCBDDCSetUpCorrection()
5282 PetscCall(MatConvert(pcbddc->coarse_psi_D, checkmattype, MAT_INITIAL_MATRIX, &coarse_psi_D)); in PCBDDCSetUpCorrection()
5283 PetscCall(MatConvert(pcbddc->coarse_psi_B, checkmattype, MAT_INITIAL_MATRIX, &coarse_psi_B)); in PCBDDCSetUpCorrection()
5285 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "--------------------------------------------… in PCBDDCSetUpCorrection()
5286 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "Check coarse sub mat computation (symmetric … in PCBDDCSetUpCorrection()
5287 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpCorrection()
5288 if (!pcbddc->symmetric_primal) { in PCBDDCSetUpCorrection()
5315 if (pcbddc->benign_n) { in PCBDDCSetUpCorrection()
5321 PetscCall(ISCreateStride(PETSC_COMM_SELF, pcbddc->benign_n, 0, 1, &is_dummy)); in PCBDDCSetUpCorrection()
5322 …PetscCall(MatCreateSubMatrix(pcbddc->benign_B0, is_dummy, pcis->is_B_local, MAT_INITIAL_MATRIX, &B… in PCBDDCSetUpCorrection()
5327 for (j = 0; j < pcbddc->benign_n; j++) { in PCBDDCSetUpCorrection()
5328 PetscInt primal_idx = pcbddc->local_primal_size - pcbddc->benign_n + j; in PCBDDCSetUpCorrection()
5329 for (i = 0; i < pcbddc->local_primal_size; i++) { in PCBDDCSetUpCorrection()
5330 data[primal_idx * pcbddc->local_primal_size + i] += data2[i * pcbddc->benign_n + j]; in PCBDDCSetUpCorrection()
5331 data[i * pcbddc->local_primal_size + primal_idx] += data2[i * pcbddc->benign_n + j]; in PCBDDCSetUpCorrection()
5342 PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer)); in PCBDDCSetUpCorrection()
5343 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d matrix e… in PCBDDCSetUpCorrection()
5346 …PetscCall(ISCreateStride(PETSC_COMM_SELF, pcbddc->local_primal_size - pcbddc->benign_n, 0, 1, &is_… in PCBDDCSetUpCorrection()
5347 …PetscCall(MatCreateSubMatrix(pcbddc->ConstraintMatrix, is_dummy, pcis->is_B_local, MAT_INITIAL_MAT… in PCBDDCSetUpCorrection()
5348 if (!pcbddc->benign_n) { /* TODO: add benign case */ in PCBDDCSetUpCorrection()
5353 PetscCall(MatDenseGetArray(pcbddc->coarse_phi_B, &data)); in PCBDDCSetUpCorrection()
5354 …PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pcis->n_B, pcbddc->local_primal_size - pcbddc->benign… in PCBDDCSetUpCorrection()
5355 PetscCall(MatDenseRestoreArray(pcbddc->coarse_phi_B, &data)); in PCBDDCSetUpCorrection()
5360 PetscCall(VecSet(mones, -1.0)); in PCBDDCSetUpCorrection()
5363 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d phi constraints e… in PCBDDCSetUpCorrection()
5364 if (!pcbddc->symmetric_primal) { in PCBDDCSetUpCorrection()
5366 PetscCall(VecSet(mones, -1.0)); in PCBDDCSetUpCorrection()
5369 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d psi constraints e… in PCBDDCSetUpCorrection()
5375 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpCorrection()
5386 if (!pcbddc->symmetric_primal) { in PCBDDCSetUpCorrection()
5397 …Y_LENGTH(filename), "details_local_coarse_mat%d_level%d.m",PetscGlobalRank,pcbddc->current_level)); in PCBDDCSetUpCorrection()
5402 if (pcbddc->coarse_phi_B) { in PCBDDCSetUpCorrection()
5403 PetscCall(PetscObjectSetName((PetscObject)pcbddc->coarse_phi_B,"phi_B")); in PCBDDCSetUpCorrection()
5404 PetscCall(MatView(pcbddc->coarse_phi_B,viewer)); in PCBDDCSetUpCorrection()
5406 if (pcbddc->coarse_phi_D) { in PCBDDCSetUpCorrection()
5407 PetscCall(PetscObjectSetName((PetscObject)pcbddc->coarse_phi_D,"phi_D")); in PCBDDCSetUpCorrection()
5408 PetscCall(MatView(pcbddc->coarse_phi_D,viewer)); in PCBDDCSetUpCorrection()
5410 if (pcbddc->coarse_psi_B) { in PCBDDCSetUpCorrection()
5411 PetscCall(PetscObjectSetName((PetscObject)pcbddc->coarse_psi_B,"psi_B")); in PCBDDCSetUpCorrection()
5412 PetscCall(MatView(pcbddc->coarse_psi_B,viewer)); in PCBDDCSetUpCorrection()
5414 if (pcbddc->coarse_psi_D) { in PCBDDCSetUpCorrection()
5415 PetscCall(PetscObjectSetName((PetscObject)pcbddc->coarse_psi_D,"psi_D")); in PCBDDCSetUpCorrection()
5416 PetscCall(MatView(pcbddc->coarse_psi_D,viewer)); in PCBDDCSetUpCorrection()
5418 PetscCall(PetscObjectSetName((PetscObject)pcbddc->local_mat,"A")); in PCBDDCSetUpCorrection()
5419 PetscCall(MatView(pcbddc->local_mat,viewer)); in PCBDDCSetUpCorrection()
5420 PetscCall(PetscObjectSetName((PetscObject)pcbddc->ConstraintMatrix,"C")); in PCBDDCSetUpCorrection()
5421 PetscCall(MatView(pcbddc->ConstraintMatrix,viewer)); in PCBDDCSetUpCorrection()
5422 PetscCall(PetscObjectSetName((PetscObject)pcis->is_I_local,"I")); in PCBDDCSetUpCorrection()
5423 PetscCall(ISView(pcis->is_I_local,viewer)); in PCBDDCSetUpCorrection()
5424 PetscCall(PetscObjectSetName((PetscObject)pcis->is_B_local,"B")); in PCBDDCSetUpCorrection()
5425 PetscCall(ISView(pcis->is_B_local,viewer)); in PCBDDCSetUpCorrection()
5426 PetscCall(PetscObjectSetName((PetscObject)pcbddc->is_R_local,"R")); in PCBDDCSetUpCorrection()
5427 PetscCall(ISView(pcbddc->is_R_local,viewer)); in PCBDDCSetUpCorrection()
5437 …PetscCall(PetscObjectTypeCompareAny((PetscObject)pcis->vec1_N, &iscuda, VECCUDA, VECMPICUDA, VECSE… in PCBDDCSetUpCorrection()
5438 …PetscCall(PetscObjectTypeCompareAny((PetscObject)pcis->vec1_N, &iship, VECHIP, VECMPIHIP, VECSEQHI… in PCBDDCSetUpCorrection()
5439 …PetscCall(PetscObjectTypeCompareAny((PetscObject)pcis->vec1_N, &iskokkos, VECKOKKOS, VECMPIKOKKOS,… in PCBDDCSetUpCorrection()
5448 …if (pcbddc->local_auxmat1) PetscCall(MatConvert(pcbddc->local_auxmat1, mtype, MAT_INPLACE_MATRIX, … in PCBDDCSetUpCorrection()
5449 …if (pcbddc->local_auxmat2) PetscCall(MatConvert(pcbddc->local_auxmat2, mtype, MAT_INPLACE_MATRIX, … in PCBDDCSetUpCorrection()
5450 …if (pcbddc->coarse_phi_B) PetscCall(MatConvert(pcbddc->coarse_phi_B, mtype, MAT_INPLACE_MATRIX, &p… in PCBDDCSetUpCorrection()
5451 …if (pcbddc->coarse_phi_D) PetscCall(MatConvert(pcbddc->coarse_phi_D, mtype, MAT_INPLACE_MATRIX, &p… in PCBDDCSetUpCorrection()
5452 …if (pcbddc->coarse_psi_B) PetscCall(MatConvert(pcbddc->coarse_psi_B, mtype, MAT_INPLACE_MATRIX, &p… in PCBDDCSetUpCorrection()
5453 …if (pcbddc->coarse_psi_D) PetscCall(MatConvert(pcbddc->coarse_psi_D, mtype, MAT_INPLACE_MATRIX, &p… in PCBDDCSetUpCorrection()
5571 (*C)->symmetric = A->symmetric; in MatPtAPWithPrefix_Private()
5572 (*C)->spd = A->spd; in MatPtAPWithPrefix_Private()
5578 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCComputeLocalMatrix()
5579 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCComputeLocalMatrix()
5587 PetscCall(MatDestroy(&pcbddc->local_mat)); in PCBDDCComputeLocalMatrix()
5588 PetscCall(MatGetSize(matis->A, &local_size, NULL)); in PCBDDCComputeLocalMatrix()
5589 if (pcbddc->mat_graph->multi_element) { in PCBDDCComputeLocalMatrix()
5592 PetscInt nsubs = pcbddc->n_local_subs; in PCBDDCComputeLocalMatrix()
5597 … = 0; i < nsubs; i++) PetscCall(ISLocalToGlobalMappingApplyIS(matis->rmapping, pcbddc->local_subs[… in PCBDDCComputeLocalMatrix()
5603 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)matis->A), local_size, 0, 1, &is_local)); in PCBDDCComputeLocalMatrix()
5604 PetscCall(ISLocalToGlobalMappingApplyIS(matis->rmapping, is_local, &is_global)); in PCBDDCComputeLocalMatrix()
5609 …PetscCall(MatCreateSubMatrices(tmats[0], nsubs, pcbddc->local_subs, pcbddc->local_subs, MAT_INITIA… in PCBDDCComputeLocalMatrix()
5613 …PetscCall(MatCreateNest(PETSC_COMM_SELF, nsubs, pcbddc->local_subs, nsubs, pcbddc->local_subs, mat… in PCBDDCComputeLocalMatrix()
5618 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)matis->A), local_size, 0, 1, &is_local)); in PCBDDCComputeLocalMatrix()
5619 PetscCall(ISLocalToGlobalMappingApplyIS(matis->rmapping, is_local, &is_global)); in PCBDDCComputeLocalMatrix()
5624 if (pcbddc->dbg_flag) { in PCBDDCComputeLocalMatrix()
5631 PetscCall(VecScatterBegin(matis->cctx, x, matis->x, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCComputeLocalMatrix()
5632 PetscCall(VecScatterEnd(matis->cctx, x, matis->x, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCComputeLocalMatrix()
5633 PetscCall(MatMult(new_mat, matis->x, matis->y)); in PCBDDCComputeLocalMatrix()
5634 if (!pcbddc->change_interior) { in PCBDDCComputeLocalMatrix()
5639 PetscCall(VecGetArrayRead(matis->x, &x)); in PCBDDCComputeLocalMatrix()
5640 PetscCall(VecGetArrayRead(matis->y, &y)); in PCBDDCComputeLocalMatrix()
5641 PetscCall(VecGetArrayRead(matis->counter, &v)); in PCBDDCComputeLocalMatrix()
5643 …f (PetscRealPart(v[i]) < 1.5 && PetscAbsScalar(x[i] - y[i]) > lerror) lerror = PetscAbsScalar(x[i]… in PCBDDCComputeLocalMatrix()
5644 PetscCall(VecRestoreArrayRead(matis->x, &x)); in PCBDDCComputeLocalMatrix()
5645 PetscCall(VecRestoreArrayRead(matis->y, &y)); in PCBDDCComputeLocalMatrix()
5646 PetscCall(VecRestoreArrayRead(matis->counter, &v)); in PCBDDCComputeLocalMatrix()
5649 if (!pcbddc->user_ChangeOfBasisMatrix || pcbddc->current_level) { in PCBDDCComputeLocalMatrix()
5656 PetscCall(VecScatterBegin(matis->rctx, matis->y, x, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCComputeLocalMatrix()
5657 PetscCall(VecScatterEnd(matis->rctx, matis->y, x, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCComputeLocalMatrix()
5658 PetscCall(VecAXPY(x, -1.0, x_change)); in PCBDDCComputeLocalMatrix()
5661 if (!pcbddc->user_ChangeOfBasisMatrix || pcbddc->current_level) { in PCBDDCComputeLocalMatrix()
5671 /* lA is present if we are setting up an inner BDDC for a saddle point FETI-DP */ in PCBDDCComputeLocalMatrix()
5675 …if (((PetscObject)pc)->prefix) PetscCall(PetscSNPrintf(ptapprefix, sizeof(ptapprefix), "%spc_bddc_… in PCBDDCComputeLocalMatrix()
5677 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij)); in PCBDDCComputeLocalMatrix()
5679 PetscCall(MatDestroy(&pcbddc->local_mat)); in PCBDDCComputeLocalMatrix()
5680 …PetscCall(MatPtAPWithPrefix_Private(matis->A, new_mat, PETSC_DEFAULT, ptapprefix, &pcbddc->local_m… in PCBDDCComputeLocalMatrix()
5690 PetscCall(MatDestroy(&pcbddc->local_mat)); in PCBDDCComputeLocalMatrix()
5691 PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &work_mat)); in PCBDDCComputeLocalMatrix()
5692 …PetscCall(MatPtAPWithPrefix_Private(work_mat, new_mat, PETSC_DEFAULT, ptapprefix, &pcbddc->local_m… in PCBDDCComputeLocalMatrix()
5702 PetscCall(MatIsSymmetricKnown(matis->A, &isset, &issym)); in PCBDDCComputeLocalMatrix()
5703 if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SYMMETRIC, issym)); in PCBDDCComputeLocalMatrix()
5710 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCSetUpLocalScatters()
5711 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCSetUpLocalScatters()
5712 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCSetUpLocalScatters()
5721 - primal space is unchanged in PCBDDCSetUpLocalScatters()
5723 …- we actually have locally some primal dofs (could not be true in multilevel or for isolated subdo… in PCBDDCSetUpLocalScatters()
5725 …- we are not in debugging mode (this is needed since there are Synchronized prints at the end of t… in PCBDDCSetUpLocalScatters()
5727 …if (!pcbddc->new_primal_space_local && pcbddc->local_primal_size && !pcbddc->dbg_flag) PetscFuncti… in PCBDDCSetUpLocalScatters()
5729 PetscCall(ISDestroy(&pcbddc->is_R_local)); in PCBDDCSetUpLocalScatters()
5730 PetscCall(VecScatterDestroy(&pcbddc->R_to_B)); in PCBDDCSetUpLocalScatters()
5731 PetscCall(VecScatterDestroy(&pcbddc->R_to_D)); in PCBDDCSetUpLocalScatters()
5732 /* Set Non-overlapping dimensions */ in PCBDDCSetUpLocalScatters()
5733 n_B = pcis->n_B; in PCBDDCSetUpLocalScatters()
5734 n_D = pcis->n - n_B; in PCBDDCSetUpLocalScatters()
5735 n_vertices = pcbddc->n_vertices; in PCBDDCSetUpLocalScatters()
5740 if (!sub_schurs || !sub_schurs->reuse_solver) { in PCBDDCSetUpLocalScatters()
5741 PetscCall(PetscMalloc1(pcis->n - n_vertices, &idx_R_local)); in PCBDDCSetUpLocalScatters()
5742 PetscCall(PetscBTCreate(pcis->n, &bitmask)); in PCBDDCSetUpLocalScatters()
5743 … for (i = 0; i < n_vertices; i++) PetscCall(PetscBTSet(bitmask, pcbddc->local_primal_ref_node[i])); in PCBDDCSetUpLocalScatters()
5745 for (i = 0, n_R = 0; i < pcis->n; i++) { in PCBDDCSetUpLocalScatters()
5749 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpLocalScatters()
5751 PetscCall(ISGetIndices(reuse_solver->is_R, (const PetscInt **)&idx_R_local)); in PCBDDCSetUpLocalScatters()
5752 PetscCall(ISGetLocalSize(reuse_solver->is_R, &n_R)); in PCBDDCSetUpLocalScatters()
5757 PetscCall(MatGetBlockSize(pcbddc->local_mat, &bs)); in PCBDDCSetUpLocalScatters()
5761 if (!sub_schurs || !sub_schurs->reuse_solver) { in PCBDDCSetUpLocalScatters()
5762 PetscCall(PetscMalloc1(pcis->n / bs, &vary)); in PCBDDCSetUpLocalScatters()
5763 PetscCall(PetscArrayzero(vary, pcis->n / bs)); in PCBDDCSetUpLocalScatters()
5766 for (i = 0; i < n_vertices; i++) vary[pcbddc->local_primal_ref_node[i] / bs]++; in PCBDDCSetUpLocalScatters()
5767 for (i = 0; i < pcis->n / bs; i++) { in PCBDDCSetUpLocalScatters()
5779 if (node != idx_R_local[bs * i + j] - j) { in PCBDDCSetUpLocalScatters()
5791 …PetscCall(ISCreateBlock(PETSC_COMM_SELF, vbs, n_R / vbs, idx_R_local, PETSC_COPY_VALUES, &pcbddc->… in PCBDDCSetUpLocalScatters()
5792 if (sub_schurs && sub_schurs->reuse_solver) { in PCBDDCSetUpLocalScatters()
5793 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpLocalScatters()
5795 PetscCall(ISRestoreIndices(reuse_solver->is_R, (const PetscInt **)&idx_R_local)); in PCBDDCSetUpLocalScatters()
5796 PetscCall(ISDestroy(&reuse_solver->is_R)); in PCBDDCSetUpLocalScatters()
5797 PetscCall(PetscObjectReference((PetscObject)pcbddc->is_R_local)); in PCBDDCSetUpLocalScatters()
5798 reuse_solver->is_R = pcbddc->is_R_local; in PCBDDCSetUpLocalScatters()
5804 if (pcbddc->dbg_flag) { in PCBDDCSetUpLocalScatters()
5805 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "--------------------------------------------… in PCBDDCSetUpLocalScatters()
5806 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpLocalScatters()
5807 PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer)); in PCBDDCSetUpLocalScatters()
5808 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d local dimensions\… in PCBDDCSetUpLocalScatters()
5809 …Printf(pcbddc->dbg_viewer, "local_size = %" PetscInt_FMT ", dirichlet_size = %" PetscInt_FMT ", bo… in PCBDDCSetUpLocalScatters()
5810 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "r_size = %" PetscInt_FMT ", v_si… in PCBDDCSetUpLocalScatters()
5811 … pcbddc->local_primal_size - n_vertices - pcbddc->benign_n, pcbddc->local_primal_size)); in PCBDDCSetUpLocalScatters()
5812 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpLocalScatters()
5815 /* VecScatters pcbddc->R_to_B and (optionally) pcbddc->R_to_D */ in PCBDDCSetUpLocalScatters()
5816 if (!sub_schurs || !sub_schurs->reuse_solver) { in PCBDDCSetUpLocalScatters()
5820 PetscCall(ISGetIndices(pcbddc->is_R_local, (const PetscInt **)&idx_R_local)); in PCBDDCSetUpLocalScatters()
5821 PetscCall(PetscMalloc1(pcis->n_B - n_vertices, &aux_array1)); in PCBDDCSetUpLocalScatters()
5822 PetscCall(PetscMalloc1(pcis->n_B - n_vertices, &aux_array2)); in PCBDDCSetUpLocalScatters()
5823 PetscCall(ISGetIndices(pcis->is_I_local, (const PetscInt **)&is_indices)); in PCBDDCSetUpLocalScatters()
5825 PetscCall(ISRestoreIndices(pcis->is_I_local, (const PetscInt **)&is_indices)); in PCBDDCSetUpLocalScatters()
5830 PetscCall(ISGetIndices(pcis->is_B_local, (const PetscInt **)&is_indices)); in PCBDDCSetUpLocalScatters()
5834 PetscCall(ISRestoreIndices(pcis->is_B_local, (const PetscInt **)&is_indices)); in PCBDDCSetUpLocalScatters()
5836 PetscCall(VecScatterCreate(pcbddc->vec1_R, is_aux1, pcis->vec1_B, is_aux2, &pcbddc->R_to_B)); in PCBDDCSetUpLocalScatters()
5840 if (pcbddc->switch_static || pcbddc->dbg_flag) { in PCBDDCSetUpLocalScatters()
5846 PetscCall(VecScatterCreate(pcbddc->vec1_R, is_aux1, pcis->vec1_D, (IS)0, &pcbddc->R_to_D)); in PCBDDCSetUpLocalScatters()
5850 PetscCall(ISRestoreIndices(pcbddc->is_R_local, (const PetscInt **)&idx_R_local)); in PCBDDCSetUpLocalScatters()
5852 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpLocalScatters()
5856 PetscCall(ISGetLocalSize(reuse_solver->is_B, &schur_size)); in PCBDDCSetUpLocalScatters()
5858 …PetscCall(VecScatterCreate(pcbddc->vec1_R, tis, pcis->vec1_B, reuse_solver->is_B, &pcbddc->R_to_B)… in PCBDDCSetUpLocalScatters()
5860 if (pcbddc->switch_static || pcbddc->dbg_flag) { in PCBDDCSetUpLocalScatters()
5862 PetscCall(VecScatterCreate(pcbddc->vec1_R, tis, pcis->vec1_D, (IS)0, &pcbddc->R_to_D)); in PCBDDCSetUpLocalScatters()
5882 Mat_IS *matis = (Mat_IS *)A->data; in MatNullSpacePropagateAny_Private()
5885 sct = matis->cctx; in MatNullSpacePropagateAny_Private()
5935 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCSetUpLocalSolvers()
5936 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCSetUpLocalSolvers()
5941 PetscScalar m_one = -1.0; in PCBDDCSetUpLocalSolvers()
5950 PetscCall(PetscLogEventBegin(PC_BDDC_LocalSolvers[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpLocalSolvers()
5952 if (!pc->setupcalled && (pcbddc->NullSpace_corr[0] || pcbddc->NullSpace_corr[2])) { in PCBDDCSetUpLocalSolvers()
5956 PetscCall(MatGetNearNullSpace(pcbddc->local_mat, &nnsp)); in PCBDDCSetUpLocalSolvers()
5957 PetscCall(MatGetNearNullSpace(pc->pmat, &gnnsp1)); in PCBDDCSetUpLocalSolvers()
5958 PetscCall(MatGetNullSpace(pc->pmat, &gnnsp2)); in PCBDDCSetUpLocalSolvers()
5961 …if (!ghas && (gnnsp1 || gnnsp2)) PetscCall(MatNullSpacePropagateAny_Private(pc->pmat, NULL, NULL)); in PCBDDCSetUpLocalSolvers()
5967 if (!pcbddc->current_level) { in PCBDDCSetUpLocalSolvers()
5968 PetscCall(PetscStrncpy(dir_prefix, ((PetscObject)pc)->prefix, sizeof(dir_prefix))); in PCBDDCSetUpLocalSolvers()
5969 PetscCall(PetscStrncpy(neu_prefix, ((PetscObject)pc)->prefix, sizeof(neu_prefix))); in PCBDDCSetUpLocalSolvers()
5973 …PetscCall(PetscSNPrintf(str_level, sizeof(str_level), "l%" PetscInt_FMT "_", pcbddc->current_level… in PCBDDCSetUpLocalSolvers()
5974 PetscCall(PetscStrlen(((PetscObject)pc)->prefix, &len)); in PCBDDCSetUpLocalSolvers()
5975 len -= 15; /* remove "pc_bddc_coarse_" */ in PCBDDCSetUpLocalSolvers()
5976 if (pcbddc->current_level > 1) len -= 3; /* remove "lX_" with X level number */ in PCBDDCSetUpLocalSolvers()
5977 if (pcbddc->current_level > 10) len -= 1; /* remove another char from level number */ in PCBDDCSetUpLocalSolvers()
5979 PetscCall(PetscStrncpy(dir_prefix, ((PetscObject)pc)->prefix, len + 1)); in PCBDDCSetUpLocalSolvers()
5980 PetscCall(PetscStrncpy(neu_prefix, ((PetscObject)pc)->prefix, len + 1)); in PCBDDCSetUpLocalSolvers()
5989 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCSetUpLocalSolvers()
5990 if (pcbddc->benign_n && !pcbddc->benign_change_explicit) { in PCBDDCSetUpLocalSolvers()
5991 …PetscCheck(sub_schurs && sub_schurs->reuse_solver, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implem… in PCBDDCSetUpLocalSolvers()
5992 if (pcbddc->dbg_flag) { in PCBDDCSetUpLocalSolvers()
5995 PetscCall(PCBDDCBenignProject(pc, pcis->is_I_local, pcis->is_I_local, &A_IIn)); in PCBDDCSetUpLocalSolvers()
5996 PetscCall(MatDestroy(&pcis->A_II)); in PCBDDCSetUpLocalSolvers()
5997 pcis->A_II = A_IIn; in PCBDDCSetUpLocalSolvers()
6000 PetscCall(MatIsSymmetricKnown(pcbddc->local_mat, &isset, &issym)); in PCBDDCSetUpLocalSolvers()
6001 if (isset) PetscCall(MatSetOption(pcis->A_II, MAT_SYMMETRIC, issym)); in PCBDDCSetUpLocalSolvers()
6003 /* Matrix for Dirichlet problem is pcis->A_II */ in PCBDDCSetUpLocalSolvers()
6004 n_D = pcis->n - pcis->n_B; in PCBDDCSetUpLocalSolvers()
6006 if (!pcbddc->ksp_D) { /* create object if not yet build */ in PCBDDCSetUpLocalSolvers()
6008 PetscCall(KSPCreate(PETSC_COMM_SELF, &pcbddc->ksp_D)); in PCBDDCSetUpLocalSolvers()
6009 PetscCall(KSPSetNestLevel(pcbddc->ksp_D, pc->kspnestlevel)); in PCBDDCSetUpLocalSolvers()
6010 PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_D, (PetscObject)pc, 1)); in PCBDDCSetUpLocalSolvers()
6012 PetscCall(KSPSetType(pcbddc->ksp_D, KSPPREONLY)); in PCBDDCSetUpLocalSolvers()
6013 PetscCall(KSPSetOptionsPrefix(pcbddc->ksp_D, dir_prefix)); in PCBDDCSetUpLocalSolvers()
6014 PetscCall(PetscObjectTypeCompare((PetscObject)pcis->pA_II, MATSEQSBAIJ, &issbaij)); in PCBDDCSetUpLocalSolvers()
6015 PetscCall(KSPGetPC(pcbddc->ksp_D, &pc_temp)); in PCBDDCSetUpLocalSolvers()
6021 PetscCall(KSPSetErrorIfNotConverged(pcbddc->ksp_D, pc->erroriffailure)); in PCBDDCSetUpLocalSolvers()
6023 PetscCall(MatSetOptionsPrefix(pcis->pA_II, ((PetscObject)pcbddc->ksp_D)->prefix)); in PCBDDCSetUpLocalSolvers()
6024 PetscCall(MatViewFromOptions(pcis->pA_II, NULL, "-mat_view")); in PCBDDCSetUpLocalSolvers()
6025 PetscCall(KSPSetOperators(pcbddc->ksp_D, pcis->A_II, pcis->pA_II)); in PCBDDCSetUpLocalSolvers()
6027 if (opts) PetscCall(KSPSetFromOptions(pcbddc->ksp_D)); in PCBDDCSetUpLocalSolvers()
6028 PetscCall(MatGetNearNullSpace(pcis->pA_II, &nnsp)); in PCBDDCSetUpLocalSolvers()
6029 if (pcbddc->NullSpace_corr[0] && !nnsp) { /* approximate solver, propagate NearNullSpace */ in PCBDDCSetUpLocalSolvers()
6030 PetscCall(MatNullSpacePropagateAny_Private(pcbddc->local_mat, pcis->is_I_local, pcis->pA_II)); in PCBDDCSetUpLocalSolvers()
6032 PetscCall(MatGetNearNullSpace(pcis->pA_II, &nnsp)); in PCBDDCSetUpLocalSolvers()
6033 PetscCall(KSPGetPC(pcbddc->ksp_D, &pc_temp)); in PCBDDCSetUpLocalSolvers()
6035 if (f && pcbddc->mat_graph->cloc && !nnsp) { in PCBDDCSetUpLocalSolvers()
6036 PetscReal *coords = pcbddc->mat_graph->coords, *scoords; in PCBDDCSetUpLocalSolvers()
6038 PetscInt cdim = pcbddc->mat_graph->cdim, nl, i, d; in PCBDDCSetUpLocalSolvers()
6040 PetscCall(ISGetLocalSize(pcis->is_I_local, &nl)); in PCBDDCSetUpLocalSolvers()
6041 PetscCall(ISGetIndices(pcis->is_I_local, &idxs)); in PCBDDCSetUpLocalSolvers()
6046 PetscCall(ISRestoreIndices(pcis->is_I_local, &idxs)); in PCBDDCSetUpLocalSolvers()
6050 if (sub_schurs && sub_schurs->reuse_solver) { in PCBDDCSetUpLocalSolvers()
6051 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpLocalSolvers()
6053 PetscCall(KSPSetPC(pcbddc->ksp_D, reuse_solver->interior_solver)); in PCBDDCSetUpLocalSolvers()
6058 PetscCall(KSPGetPC(pcbddc->ksp_D, &pc_temp)); in PCBDDCSetUpLocalSolvers()
6061 PetscCall(KSPSetUp(pcbddc->ksp_D)); in PCBDDCSetUpLocalSolvers()
6063 PetscCall(PetscObjectReference((PetscObject)pcbddc->ksp_D)); in PCBDDCSetUpLocalSolvers()
6064 PetscCall(KSPDestroy(&pcis->ksp_D)); in PCBDDCSetUpLocalSolvers()
6065 pcis->ksp_D = pcbddc->ksp_D; in PCBDDCSetUpLocalSolvers()
6071 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCSetUpLocalSolvers()
6074 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCSetUpLocalSolvers()
6077 if (sub_schurs && sub_schurs->reuse_solver) { in PCBDDCSetUpLocalSolvers()
6081 PetscCall(PetscObjectQuery((PetscObject)sub_schurs->A, "__KSPFETIDP_iP", (PetscObject *)&iP)); in PCBDDCSetUpLocalSolvers()
6084 /* Matrix for Neumann problem is A_RR -> we need to create/reuse it at this point */ in PCBDDCSetUpLocalSolvers()
6085 PetscCall(ISGetSize(pcbddc->is_R_local, &n_R)); in PCBDDCSetUpLocalSolvers()
6086 if (pcbddc->ksp_R) { /* already created ksp */ in PCBDDCSetUpLocalSolvers()
6088 PetscCall(KSPGetOperators(pcbddc->ksp_R, NULL, &A_RR)); in PCBDDCSetUpLocalSolvers()
6092 PetscCall(KSPReset(pcbddc->ksp_R)); in PCBDDCSetUpLocalSolvers()
6096 …if (pcbddc->new_primal_space_local) { /* we are not sure the matrix will have the same nonzero pat… in PCBDDCSetUpLocalSolvers()
6104 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { in PCBDDCSetUpLocalSolvers()
6111 /* convert pcbddc->local_mat if needed later in PCBDDCSetUpCorrection in PCBDDCSetUpLocalSolvers()
6113 PetscCall(MatGetBlockSize(pcbddc->local_mat, &mbs)); in PCBDDCSetUpLocalSolvers()
6114 PetscCall(ISGetBlockSize(pcbddc->is_R_local, &ibs)); in PCBDDCSetUpLocalSolvers()
6115 PetscCall(PetscObjectTypeCompare((PetscObject)pcbddc->local_mat, MATSEQSBAIJ, &issbaij)); in PCBDDCSetUpLocalSolvers()
6117 if (matis->A == pcbddc->local_mat) { in PCBDDCSetUpLocalSolvers()
6118 PetscCall(MatDestroy(&pcbddc->local_mat)); in PCBDDCSetUpLocalSolvers()
6119 PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &pcbddc->local_mat)); in PCBDDCSetUpLocalSolvers()
6121 PetscCall(MatConvert(pcbddc->local_mat, MATSEQAIJ, MAT_INPLACE_MATRIX, &pcbddc->local_mat)); in PCBDDCSetUpLocalSolvers()
6123 } else if (issbaij) { /* need to convert to BAIJ to get off-diagonal blocks */ in PCBDDCSetUpLocalSolvers()
6124 if (matis->A == pcbddc->local_mat) { in PCBDDCSetUpLocalSolvers()
6125 PetscCall(MatDestroy(&pcbddc->local_mat)); in PCBDDCSetUpLocalSolvers()
6126 …PetscCall(MatConvert(matis->A, mbs > 1 ? MATSEQBAIJ : MATSEQAIJ, MAT_INITIAL_MATRIX, &pcbddc->loca… in PCBDDCSetUpLocalSolvers()
6128 …PetscCall(MatConvert(pcbddc->local_mat, mbs > 1 ? MATSEQBAIJ : MATSEQAIJ, MAT_INPLACE_MATRIX, &pcb… in PCBDDCSetUpLocalSolvers()
6133 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpLocalSolvers()
6135 if (pcbddc->dbg_flag) { /* we need A_RR to test the solver later */ in PCBDDCSetUpLocalSolvers()
6137 … if (reuse_solver->benign_n) { /* we are not using the explicit change of basis on the pressures */ in PCBDDCSetUpLocalSolvers()
6138 PetscCall(PCBDDCBenignProject(pc, pcbddc->is_R_local, pcbddc->is_R_local, &A_RR)); in PCBDDCSetUpLocalSolvers()
6140 …PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcbddc->is_R_local, pcbddc->is_R_local, MAT_INITIA… in PCBDDCSetUpLocalSolvers()
6144 PetscCall(PCGetOperators(reuse_solver->correction_solver, &A_RR, NULL)); in PCBDDCSetUpLocalSolvers()
6148 …PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcbddc->is_R_local, pcbddc->is_R_local, reuse, &A_… in PCBDDCSetUpLocalSolvers()
6150 PetscCall(MatIsSymmetricKnown(pcbddc->local_mat, &isset, &issym)); in PCBDDCSetUpLocalSolvers()
6153 if (!pcbddc->ksp_R) { /* create object if not present */ in PCBDDCSetUpLocalSolvers()
6155 PetscCall(KSPCreate(PETSC_COMM_SELF, &pcbddc->ksp_R)); in PCBDDCSetUpLocalSolvers()
6156 PetscCall(KSPSetNestLevel(pcbddc->ksp_R, pc->kspnestlevel)); in PCBDDCSetUpLocalSolvers()
6157 PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_R, (PetscObject)pc, 1)); in PCBDDCSetUpLocalSolvers()
6159 PetscCall(KSPSetType(pcbddc->ksp_R, KSPPREONLY)); in PCBDDCSetUpLocalSolvers()
6160 PetscCall(KSPSetOptionsPrefix(pcbddc->ksp_R, neu_prefix)); in PCBDDCSetUpLocalSolvers()
6161 PetscCall(KSPGetPC(pcbddc->ksp_R, &pc_temp)); in PCBDDCSetUpLocalSolvers()
6168 PetscCall(KSPSetErrorIfNotConverged(pcbddc->ksp_R, pc->erroriffailure)); in PCBDDCSetUpLocalSolvers()
6170 PetscCall(MatSetOptionsPrefix(A_RR, ((PetscObject)pcbddc->ksp_R)->prefix)); in PCBDDCSetUpLocalSolvers()
6171 PetscCall(MatViewFromOptions(A_RR, NULL, "-mat_view")); in PCBDDCSetUpLocalSolvers()
6172 PetscCall(KSPSetOperators(pcbddc->ksp_R, A_RR, A_RR)); in PCBDDCSetUpLocalSolvers()
6174 PetscCall(KSPSetFromOptions(pcbddc->ksp_R)); in PCBDDCSetUpLocalSolvers()
6177 if (pcbddc->NullSpace_corr[2] && !nnsp) { /* approximate solver, propagate NearNullSpace */ in PCBDDCSetUpLocalSolvers()
6178 PetscCall(MatNullSpacePropagateAny_Private(pcbddc->local_mat, pcbddc->is_R_local, A_RR)); in PCBDDCSetUpLocalSolvers()
6181 PetscCall(KSPGetPC(pcbddc->ksp_R, &pc_temp)); in PCBDDCSetUpLocalSolvers()
6183 if (f && pcbddc->mat_graph->cloc && !nnsp) { in PCBDDCSetUpLocalSolvers()
6184 PetscReal *coords = pcbddc->mat_graph->coords, *scoords; in PCBDDCSetUpLocalSolvers()
6186 PetscInt cdim = pcbddc->mat_graph->cdim, nl, i, d; in PCBDDCSetUpLocalSolvers()
6188 PetscCall(ISGetLocalSize(pcbddc->is_R_local, &nl)); in PCBDDCSetUpLocalSolvers()
6189 PetscCall(ISGetIndices(pcbddc->is_R_local, &idxs)); in PCBDDCSetUpLocalSolvers()
6194 PetscCall(ISRestoreIndices(pcbddc->is_R_local, &idxs)); in PCBDDCSetUpLocalSolvers()
6201 PetscCall(KSPGetPC(pcbddc->ksp_R, &pc_temp)); in PCBDDCSetUpLocalSolvers()
6206 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSetUpLocalSolvers()
6208 PetscCall(KSPSetPC(pcbddc->ksp_R, reuse_solver->correction_solver)); in PCBDDCSetUpLocalSolvers()
6210 PetscCall(KSPSetUp(pcbddc->ksp_R)); in PCBDDCSetUpLocalSolvers()
6213 if (pcbddc->dbg_flag) { in PCBDDCSetUpLocalSolvers()
6214 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpLocalSolvers()
6215 PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer)); in PCBDDCSetUpLocalSolvers()
6216 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "--------------------------------------------… in PCBDDCSetUpLocalSolvers()
6218 PetscCall(PetscLogEventEnd(PC_BDDC_LocalSolvers[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpLocalSolvers()
6221 if (pcbddc->NullSpace_corr[0]) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE)); in PCBDDCSetUpLocalSolvers()
6222 …dirichlet && pcbddc->NullSpace_corr[0] && !pcbddc->switch_static) PetscCall(PCBDDCNullSpaceAssembl… in PCBDDCSetUpLocalSolvers()
6223 …if (neumann && pcbddc->NullSpace_corr[2]) PetscCall(PCBDDCNullSpaceAssembleCorrection(pc, PETSC_FA… in PCBDDCSetUpLocalSolvers()
6225 if (pcbddc->dbg_flag) { in PCBDDCSetUpLocalSolvers()
6227 PetscCall(VecSetRandom(pcis->vec1_D, NULL)); in PCBDDCSetUpLocalSolvers()
6228 PetscCall(MatMult(pcis->A_II, pcis->vec1_D, pcis->vec2_D)); in PCBDDCSetUpLocalSolvers()
6229 PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec2_D, pcis->vec2_D)); in PCBDDCSetUpLocalSolvers()
6230 PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D)); in PCBDDCSetUpLocalSolvers()
6231 PetscCall(VecAXPY(pcis->vec1_D, m_one, pcis->vec2_D)); in PCBDDCSetUpLocalSolvers()
6232 PetscCall(VecNorm(pcis->vec1_D, NORM_INFINITY, &value)); in PCBDDCSetUpLocalSolvers()
6233 …intf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet solve (%s) = % 1.14e \n", Pe… in PCBDDCSetUpLocalSolvers()
6234 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpLocalSolvers()
6237 PetscCall(VecSetRandom(pcbddc->vec1_R, NULL)); in PCBDDCSetUpLocalSolvers()
6238 PetscCall(MatMult(A_RR, pcbddc->vec1_R, pcbddc->vec2_R)); in PCBDDCSetUpLocalSolvers()
6239 PetscCall(KSPSolve(pcbddc->ksp_R, pcbddc->vec2_R, pcbddc->vec2_R)); in PCBDDCSetUpLocalSolvers()
6240 PetscCall(KSPCheckSolve(pcbddc->ksp_R, pc, pcbddc->vec2_R)); in PCBDDCSetUpLocalSolvers()
6241 PetscCall(VecAXPY(pcbddc->vec1_R, m_one, pcbddc->vec2_R)); in PCBDDCSetUpLocalSolvers()
6242 PetscCall(VecNorm(pcbddc->vec1_R, NORM_INFINITY, &value)); in PCBDDCSetUpLocalSolvers()
6243 …rintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann solve (%s) = % 1.14e\n", Pets… in PCBDDCSetUpLocalSolvers()
6244 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpLocalSolvers()
6254 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCSolveSubstructureCorrection()
6255 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCSolveSubstructureCorrection()
6256 …PetscBool reuse_solver = sub_schurs ? (sub_schurs->reuse_solver ? PETSC_TRUE : PETSC_FALSE) … in PCBDDCSolveSubstructureCorrection()
6259 if (!reuse_solver) PetscCall(VecSet(pcbddc->vec1_R, 0.)); in PCBDDCSolveSubstructureCorrection()
6260 if (!pcbddc->switch_static) { in PCBDDCSolveSubstructureCorrection()
6261 if (applytranspose && pcbddc->local_auxmat1) { in PCBDDCSolveSubstructureCorrection()
6262 PetscCall(MatMultTranspose(pcbddc->local_auxmat2, inout_B, pcbddc->vec1_C)); in PCBDDCSolveSubstructureCorrection()
6263 PetscCall(MatMultTransposeAdd(pcbddc->local_auxmat1, pcbddc->vec1_C, inout_B, inout_B)); in PCBDDCSolveSubstructureCorrection()
6266 …PetscCall(VecScatterBegin(pcbddc->R_to_B, inout_B, pcbddc->vec1_R, INSERT_VALUES, SCATTER_REVERSE)… in PCBDDCSolveSubstructureCorrection()
6267 … PetscCall(VecScatterEnd(pcbddc->R_to_B, inout_B, pcbddc->vec1_R, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCSolveSubstructureCorrection()
6269 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSolveSubstructureCorrection()
6271 …PetscCall(VecScatterBegin(reuse_solver->correction_scatter_B, inout_B, reuse_solver->rhs_B, INSERT… in PCBDDCSolveSubstructureCorrection()
6272 …PetscCall(VecScatterEnd(reuse_solver->correction_scatter_B, inout_B, reuse_solver->rhs_B, INSERT_V… in PCBDDCSolveSubstructureCorrection()
6275 …PetscCall(VecScatterBegin(pcbddc->R_to_B, inout_B, pcbddc->vec1_R, INSERT_VALUES, SCATTER_REVERSE)… in PCBDDCSolveSubstructureCorrection()
6276 … PetscCall(VecScatterEnd(pcbddc->R_to_B, inout_B, pcbddc->vec1_R, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCSolveSubstructureCorrection()
6277 …PetscCall(VecScatterBegin(pcbddc->R_to_D, inout_D, pcbddc->vec1_R, INSERT_VALUES, SCATTER_REVERSE)… in PCBDDCSolveSubstructureCorrection()
6278 … PetscCall(VecScatterEnd(pcbddc->R_to_D, inout_D, pcbddc->vec1_R, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCSolveSubstructureCorrection()
6279 if (applytranspose && pcbddc->local_auxmat1) { in PCBDDCSolveSubstructureCorrection()
6280 PetscCall(MatMultTranspose(pcbddc->local_auxmat2, pcbddc->vec1_R, pcbddc->vec1_C)); in PCBDDCSolveSubstructureCorrection()
6281 PetscCall(MatMultTransposeAdd(pcbddc->local_auxmat1, pcbddc->vec1_C, inout_B, inout_B)); in PCBDDCSolveSubstructureCorrection()
6282 …PetscCall(VecScatterBegin(pcbddc->R_to_B, inout_B, pcbddc->vec1_R, INSERT_VALUES, SCATTER_REVERSE)… in PCBDDCSolveSubstructureCorrection()
6283 … PetscCall(VecScatterEnd(pcbddc->R_to_B, inout_B, pcbddc->vec1_R, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCSolveSubstructureCorrection()
6286 PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][1], pc, 0, 0, 0)); in PCBDDCSolveSubstructureCorrection()
6287 if (!reuse_solver || pcbddc->switch_static) { in PCBDDCSolveSubstructureCorrection()
6289 PetscCall(KSPSolveTranspose(pcbddc->ksp_R, pcbddc->vec1_R, pcbddc->vec1_R)); in PCBDDCSolveSubstructureCorrection()
6291 PetscCall(KSPSolve(pcbddc->ksp_R, pcbddc->vec1_R, pcbddc->vec1_R)); in PCBDDCSolveSubstructureCorrection()
6293 PetscCall(KSPCheckSolve(pcbddc->ksp_R, pc, pcbddc->vec1_R)); in PCBDDCSolveSubstructureCorrection()
6295 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSolveSubstructureCorrection()
6298 …cCall(MatFactorSolveSchurComplementTranspose(reuse_solver->F, reuse_solver->rhs_B, reuse_solver->s… in PCBDDCSolveSubstructureCorrection()
6300 …PetscCall(MatFactorSolveSchurComplement(reuse_solver->F, reuse_solver->rhs_B, reuse_solver->sol_B)… in PCBDDCSolveSubstructureCorrection()
6303 PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][1], pc, 0, 0, 0)); in PCBDDCSolveSubstructureCorrection()
6305 if (!pcbddc->switch_static) { in PCBDDCSolveSubstructureCorrection()
6307 …PetscCall(VecScatterBegin(pcbddc->R_to_B, pcbddc->vec1_R, inout_B, INSERT_VALUES, SCATTER_FORWARD)… in PCBDDCSolveSubstructureCorrection()
6308 … PetscCall(VecScatterEnd(pcbddc->R_to_B, pcbddc->vec1_R, inout_B, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSolveSubstructureCorrection()
6310 PCBDDCReuseSolvers reuse_solver = sub_schurs->reuse_solver; in PCBDDCSolveSubstructureCorrection()
6312 …PetscCall(VecScatterBegin(reuse_solver->correction_scatter_B, reuse_solver->sol_B, inout_B, INSERT… in PCBDDCSolveSubstructureCorrection()
6313 …PetscCall(VecScatterEnd(reuse_solver->correction_scatter_B, reuse_solver->sol_B, inout_B, INSERT_V… in PCBDDCSolveSubstructureCorrection()
6315 if (!applytranspose && pcbddc->local_auxmat1) { in PCBDDCSolveSubstructureCorrection()
6316 PetscCall(MatMult(pcbddc->local_auxmat1, inout_B, pcbddc->vec1_C)); in PCBDDCSolveSubstructureCorrection()
6317 PetscCall(MatMultAdd(pcbddc->local_auxmat2, pcbddc->vec1_C, inout_B, inout_B)); in PCBDDCSolveSubstructureCorrection()
6320 …PetscCall(VecScatterBegin(pcbddc->R_to_B, pcbddc->vec1_R, inout_B, INSERT_VALUES, SCATTER_FORWARD)… in PCBDDCSolveSubstructureCorrection()
6321 … PetscCall(VecScatterEnd(pcbddc->R_to_B, pcbddc->vec1_R, inout_B, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSolveSubstructureCorrection()
6322 …PetscCall(VecScatterBegin(pcbddc->R_to_D, pcbddc->vec1_R, inout_D, INSERT_VALUES, SCATTER_FORWARD)… in PCBDDCSolveSubstructureCorrection()
6323 … PetscCall(VecScatterEnd(pcbddc->R_to_D, pcbddc->vec1_R, inout_D, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSolveSubstructureCorrection()
6324 if (!applytranspose && pcbddc->local_auxmat1) { in PCBDDCSolveSubstructureCorrection()
6325 PetscCall(MatMult(pcbddc->local_auxmat1, inout_B, pcbddc->vec1_C)); in PCBDDCSolveSubstructureCorrection()
6326 PetscCall(MatMultAdd(pcbddc->local_auxmat2, pcbddc->vec1_C, pcbddc->vec1_R, pcbddc->vec1_R)); in PCBDDCSolveSubstructureCorrection()
6328 …PetscCall(VecScatterBegin(pcbddc->R_to_B, pcbddc->vec1_R, inout_B, INSERT_VALUES, SCATTER_FORWARD)… in PCBDDCSolveSubstructureCorrection()
6329 … PetscCall(VecScatterEnd(pcbddc->R_to_B, pcbddc->vec1_R, inout_B, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSolveSubstructureCorrection()
6330 …PetscCall(VecScatterBegin(pcbddc->R_to_D, pcbddc->vec1_R, inout_D, INSERT_VALUES, SCATTER_FORWARD)… in PCBDDCSolveSubstructureCorrection()
6331 … PetscCall(VecScatterEnd(pcbddc->R_to_D, pcbddc->vec1_R, inout_D, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCSolveSubstructureCorrection()
6339 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCApplyInterfacePreconditioner()
6340 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCApplyInterfacePreconditioner()
6345 if (!pcbddc->benign_apply_coarse_only) { in PCBDDCApplyInterfacePreconditioner()
6347 PetscCall(MatMultTranspose(pcbddc->coarse_phi_B, pcis->vec1_B, pcbddc->vec1_P)); in PCBDDCApplyInterfacePreconditioner()
6348 …if (pcbddc->switch_static) PetscCall(MatMultTransposeAdd(pcbddc->coarse_phi_D, pcis->vec1_D, pcbdd… in PCBDDCApplyInterfacePreconditioner()
6350 PetscCall(MatMultTranspose(pcbddc->coarse_psi_B, pcis->vec1_B, pcbddc->vec1_P)); in PCBDDCApplyInterfacePreconditioner()
6351 …if (pcbddc->switch_static) PetscCall(MatMultTransposeAdd(pcbddc->coarse_psi_D, pcis->vec1_D, pcbdd… in PCBDDCApplyInterfacePreconditioner()
6354 PetscCall(VecSet(pcbddc->vec1_P, zero)); in PCBDDCApplyInterfacePreconditioner()
6358 if (pcbddc->benign_n) { in PCBDDCApplyInterfacePreconditioner()
6362 PetscCall(VecGetArray(pcbddc->vec1_P, &array)); in PCBDDCApplyInterfacePreconditioner()
6363 …for (j = 0; j < pcbddc->benign_n; j++) array[pcbddc->local_primal_size - pcbddc->benign_n + j] += … in PCBDDCApplyInterfacePreconditioner()
6364 PetscCall(VecRestoreArray(pcbddc->vec1_P, &array)); in PCBDDCApplyInterfacePreconditioner()
6368 PetscCall(VecSet(pcbddc->coarse_vec, zero)); in PCBDDCApplyInterfacePreconditioner()
6372 /* Coarse solution -> rhs and sol updated inside PCBDDCScattarCoarseDataBegin/End */ in PCBDDCApplyInterfacePreconditioner()
6373 PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][2], pc, 0, 0, 0)); in PCBDDCApplyInterfacePreconditioner()
6374 if (pcbddc->coarse_ksp) { in PCBDDCApplyInterfacePreconditioner()
6380 if (pcbddc->benign_have_null) { in PCBDDCApplyInterfacePreconditioner()
6383 PetscCall(KSPGetPC(pcbddc->coarse_ksp, &coarse_pc)); in PCBDDCApplyInterfacePreconditioner()
6386 if (isbddc && pcbddc->benign_apply_coarse_only && !pcbddc->benign_skip_correction) { in PCBDDCApplyInterfacePreconditioner()
6387 PC_BDDC *coarsepcbddc = (PC_BDDC *)coarse_pc->data; in PCBDDCApplyInterfacePreconditioner()
6388 coarsepcbddc->benign_skip_correction = PETSC_FALSE; in PCBDDCApplyInterfacePreconditioner()
6389 coarsepcbddc->benign_apply_coarse_only = PETSC_TRUE; in PCBDDCApplyInterfacePreconditioner()
6392 PetscCall(KSPGetRhs(pcbddc->coarse_ksp, &rhs)); in PCBDDCApplyInterfacePreconditioner()
6393 PetscCall(KSPGetSolution(pcbddc->coarse_ksp, &sol)); in PCBDDCApplyInterfacePreconditioner()
6394 PetscCall(KSPGetOperators(pcbddc->coarse_ksp, &coarse_mat, NULL)); in PCBDDCApplyInterfacePreconditioner()
6396 …PetscCheck(!pcbddc->benign_apply_coarse_only, PetscObjectComm((PetscObject)pcbddc->coarse_ksp), PE… in PCBDDCApplyInterfacePreconditioner()
6397 PetscCall(KSPSolveTranspose(pcbddc->coarse_ksp, rhs, sol)); in PCBDDCApplyInterfacePreconditioner()
6398 PetscCall(KSPCheckSolve(pcbddc->coarse_ksp, pc, sol)); in PCBDDCApplyInterfacePreconditioner()
6403 …if (pcbddc->benign_apply_coarse_only && isbddc) { /* need just to apply the coarse preconditioner … in PCBDDCApplyInterfacePreconditioner()
6407 PetscCall(KSPGetPC(pcbddc->coarse_ksp, &coarse_pc)); in PCBDDCApplyInterfacePreconditioner()
6408 PetscCall(PCPreSolve(coarse_pc, pcbddc->coarse_ksp)); in PCBDDCApplyInterfacePreconditioner()
6410 PetscCall(PCPostSolve(coarse_pc, pcbddc->coarse_ksp)); in PCBDDCApplyInterfacePreconditioner()
6412 PetscCall(KSPSolve(pcbddc->coarse_ksp, rhs, sol)); in PCBDDCApplyInterfacePreconditioner()
6413 PetscCall(KSPCheckSolve(pcbddc->coarse_ksp, pc, sol)); in PCBDDCApplyInterfacePreconditioner()
6418 if (pcbddc->benign_have_null && isbddc) { in PCBDDCApplyInterfacePreconditioner()
6422 PetscCall(KSPGetPC(pcbddc->coarse_ksp, &coarse_pc)); in PCBDDCApplyInterfacePreconditioner()
6423 coarsepcbddc = (PC_BDDC *)coarse_pc->data; in PCBDDCApplyInterfacePreconditioner()
6424 coarsepcbddc->benign_skip_correction = PETSC_TRUE; in PCBDDCApplyInterfacePreconditioner()
6425 coarsepcbddc->benign_apply_coarse_only = PETSC_FALSE; in PCBDDCApplyInterfacePreconditioner()
6428 PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][2], pc, 0, 0, 0)); in PCBDDCApplyInterfacePreconditioner()
6431 …if (!pcbddc->benign_apply_coarse_only) PetscCall(PCBDDCSolveSubstructureCorrection(pc, pcis->vec1_… in PCBDDCApplyInterfacePreconditioner()
6437 if (!pcbddc->benign_apply_coarse_only) { in PCBDDCApplyInterfacePreconditioner()
6439 PetscCall(MatMultAdd(pcbddc->coarse_psi_B, pcbddc->vec1_P, pcis->vec1_B, pcis->vec1_B)); in PCBDDCApplyInterfacePreconditioner()
6440 …if (pcbddc->switch_static) PetscCall(MatMultAdd(pcbddc->coarse_psi_D, pcbddc->vec1_P, pcis->vec1_D… in PCBDDCApplyInterfacePreconditioner()
6442 PetscCall(MatMultAdd(pcbddc->coarse_phi_B, pcbddc->vec1_P, pcis->vec1_B, pcis->vec1_B)); in PCBDDCApplyInterfacePreconditioner()
6443 …if (pcbddc->switch_static) PetscCall(MatMultAdd(pcbddc->coarse_phi_D, pcbddc->vec1_P, pcis->vec1_D… in PCBDDCApplyInterfacePreconditioner()
6446 if (pcbddc->benign_n) { in PCBDDCApplyInterfacePreconditioner()
6450 PetscCall(VecGetArray(pcbddc->vec1_P, &array)); in PCBDDCApplyInterfacePreconditioner()
6451 …for (j = 0; j < pcbddc->benign_n; j++) pcbddc->benign_p0[j] = array[pcbddc->local_primal_size - pc… in PCBDDCApplyInterfacePreconditioner()
6452 PetscCall(VecRestoreArray(pcbddc->vec1_P, &array)); in PCBDDCApplyInterfacePreconditioner()
6456 PetscCall(MatMult(pcbddc->coarse_psi_B, pcbddc->vec1_P, pcis->vec1_B)); in PCBDDCApplyInterfacePreconditioner()
6458 PetscCall(MatMult(pcbddc->coarse_phi_B, pcbddc->vec1_P, pcis->vec1_B)); in PCBDDCApplyInterfacePreconditioner()
6466 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCScatterCoarseDataBegin()
6471 if (smode == SCATTER_REVERSE) { /* from global to local -> get data from coarse solution */ in PCBDDCScatterCoarseDataBegin()
6472 from = pcbddc->coarse_vec; in PCBDDCScatterCoarseDataBegin()
6473 to = pcbddc->vec1_P; in PCBDDCScatterCoarseDataBegin()
6474 if (pcbddc->coarse_ksp) { /* get array from coarse processes */ in PCBDDCScatterCoarseDataBegin()
6477 PetscCall(KSPGetRhs(pcbddc->coarse_ksp, &tvec)); in PCBDDCScatterCoarseDataBegin()
6479 PetscCall(KSPGetSolution(pcbddc->coarse_ksp, &tvec)); in PCBDDCScatterCoarseDataBegin()
6484 } else { /* from local to global -> put data in coarse right-hand side */ in PCBDDCScatterCoarseDataBegin()
6485 from = pcbddc->vec1_P; in PCBDDCScatterCoarseDataBegin()
6486 to = pcbddc->coarse_vec; in PCBDDCScatterCoarseDataBegin()
6488 PetscCall(VecScatterBegin(pcbddc->coarse_loc_to_glob, from, to, imode, smode)); in PCBDDCScatterCoarseDataBegin()
6494 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCScatterCoarseDataEnd()
6499 if (smode == SCATTER_REVERSE) { /* from global to local -> get data from coarse solution */ in PCBDDCScatterCoarseDataEnd()
6500 from = pcbddc->coarse_vec; in PCBDDCScatterCoarseDataEnd()
6501 to = pcbddc->vec1_P; in PCBDDCScatterCoarseDataEnd()
6502 } else { /* from local to global -> put data in coarse right-hand side */ in PCBDDCScatterCoarseDataEnd()
6503 from = pcbddc->vec1_P; in PCBDDCScatterCoarseDataEnd()
6504 to = pcbddc->coarse_vec; in PCBDDCScatterCoarseDataEnd()
6506 PetscCall(VecScatterEnd(pcbddc->coarse_loc_to_glob, from, to, imode, smode)); in PCBDDCScatterCoarseDataEnd()
6508 if (pcbddc->coarse_ksp) { /* get array from coarse processes */ in PCBDDCScatterCoarseDataEnd()
6511 PetscCall(KSPGetRhs(pcbddc->coarse_ksp, &tvec)); in PCBDDCScatterCoarseDataEnd()
6517 if (pcbddc->coarse_ksp) { /* restore array of pcbddc->coarse_vec */ in PCBDDCScatterCoarseDataEnd()
6526 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCConstraintsSetUp()
6527 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCConstraintsSetUp()
6528 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCConstraintsSetUp()
6559 PetscCall(MatDestroy(&pcbddc->ChangeOfBasisMatrix)); in PCBDDCConstraintsSetUp()
6560 PetscCall(MatDestroy(&pcbddc->ConstraintMatrix)); in PCBDDCConstraintsSetUp()
6561 PetscCall(MatDestroy(&pcbddc->switch_static_change)); in PCBDDCConstraintsSetUp()
6563 olocal_primal_size = pcbddc->local_primal_size; in PCBDDCConstraintsSetUp()
6564 olocal_primal_size_cc = pcbddc->local_primal_size_cc; in PCBDDCConstraintsSetUp()
6566 …PetscCall(PetscArraycpy(olocal_primal_ref_node, pcbddc->local_primal_ref_node, olocal_primal_size_… in PCBDDCConstraintsSetUp()
6567 …PetscCall(PetscArraycpy(olocal_primal_ref_mult, pcbddc->local_primal_ref_mult, olocal_primal_size_… in PCBDDCConstraintsSetUp()
6568 PetscCall(PetscFree2(pcbddc->local_primal_ref_node, pcbddc->local_primal_ref_mult)); in PCBDDCConstraintsSetUp()
6569 PetscCall(PetscFree(pcbddc->primal_indices_local_idxs)); in PCBDDCConstraintsSetUp()
6571 if (!pcbddc->adaptive_selection) { in PCBDDCConstraintsSetUp()
6592 …PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, &n_ISForFaces, &ISForFaces, &n_ISForEdges,… in PCBDDCConstraintsSetUp()
6598 if (pcbddc->dbg_flag && (!pcbddc->sub_schurs || pcbddc->sub_schurs_rebuild)) { in PCBDDCConstraintsSetUp()
6599 …if (!pcbddc->dbg_viewer) pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc… in PCBDDCConstraintsSetUp()
6600 PetscCall(PCBDDCGraphASCIIView(pcbddc->mat_graph, pcbddc->dbg_flag, pcbddc->dbg_viewer)); in PCBDDCConstraintsSetUp()
6601 PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer)); in PCBDDCConstraintsSetUp()
6602 …scViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "---------------------------------------------… in PCBDDCConstraintsSetUp()
6603 …dPrintf(pcbddc->dbg_viewer, "Subdomain %04d got %02" PetscInt_FMT " local candidate vertices (%d)\… in PCBDDCConstraintsSetUp()
6604 …Printf(pcbddc->dbg_viewer, "Subdomain %04d got %02" PetscInt_FMT " local candidate edges (%d)\n… in PCBDDCConstraintsSetUp()
6605 …Printf(pcbddc->dbg_viewer, "Subdomain %04d got %02" PetscInt_FMT " local candidate faces (%d)\n… in PCBDDCConstraintsSetUp()
6606 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCConstraintsSetUp()
6607 PetscCall(PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer)); in PCBDDCConstraintsSetUp()
6610 if (!pcbddc->use_vertices) n_vertices = 0; in PCBDDCConstraintsSetUp()
6611 if (!pcbddc->use_edges) n_ISForEdges = 0; in PCBDDCConstraintsSetUp()
6612 if (!pcbddc->use_faces) n_ISForFaces = 0; in PCBDDCConstraintsSetUp()
6615 if (pcbddc->use_nnsp) PetscCall(MatGetNearNullSpace(pc->pmat, &nearnullsp)); in PCBDDCConstraintsSetUp()
6621 PetscCall(MatNullSpaceDestroy(&pcbddc->onearnullspace)); in PCBDDCConstraintsSetUp()
6622 PetscCall(PetscFree(pcbddc->onearnullvecs_state)); in PCBDDCConstraintsSetUp()
6625 pcbddc->onearnullspace = nearnullsp; in PCBDDCConstraintsSetUp()
6626 PetscCall(PetscMalloc1(nnsp_size, &pcbddc->onearnullvecs_state)); in PCBDDCConstraintsSetUp()
6627 …ze; i++) PetscCall(PetscObjectStateGet((PetscObject)nearnullvecs[i], &pcbddc->onearnullvecs_state[… in PCBDDCConstraintsSetUp()
6638 - Indices for connected component i stored at "constraints_idxs + constraints_idxs_ptr[i]" in PCBDDCConstraintsSetUp()
6639 …- Values for constraints on connected component i stored at "constraints_data + constraints_data_p… in PCBDDCConstraintsSetUp()
6657 used_is = ISForFaces[i - n_ISForEdges]; in PCBDDCConstraintsSetUp()
6668 PetscCall(VecDuplicate(pcis->vec1_N, &localnearnullsp[k])); in PCBDDCConstraintsSetUp()
6669 …PetscCall(VecScatterBegin(matis->rctx, nearnullvecs[k], localnearnullsp[k], INSERT_VALUES, SCATTER… in PCBDDCConstraintsSetUp()
6670 …PetscCall(VecScatterEnd(matis->rctx, nearnullvecs[k], localnearnullsp[k], INSERT_VALUES, SCATTER_F… in PCBDDCConstraintsSetUp()
6675 …if (n_ISForFaces + n_ISForEdges && max_constraints > 1 && !pcbddc->use_nnsp_true) skip_lapack = PE… in PCBDDCConstraintsSetUp()
6689 /* now we evaluate the optimal workspace using query with lwork=-1 */ in PCBDDCConstraintsSetUp()
6692 lwork = -1; in PCBDDCConstraintsSetUp()
6715 /* now we evaluate the optimal workspace using query with lwork=-1 */ in PCBDDCConstraintsSetUp()
6716 lwork = -1; in PCBDDCConstraintsSetUp()
6762 boolforchange = pcbddc->use_change_of_basis; /* change or not the basis on the edge */ in PCBDDCConstraintsSetUp()
6764 used_is = ISForFaces[ncc - n_ISForEdges]; in PCBDDCConstraintsSetUp()
6765 …boolforchange = (PetscBool)(pcbddc->use_change_of_basis && pcbddc->use_change_on_faces); /* change… in PCBDDCConstraintsSetUp()
6778 if (!pcbddc->use_nnsp_true) { in PCBDDCConstraintsSetUp()
6809 if (!pcbddc->use_nnsp_true && temp_constraints) { in PCBDDCConstraintsSetUp()
6822 /* SVD: Y = U*S*V^H -> U (eigenvectors of Y*Y^H) = Y*V*(S)^\dag in PCBDDCConstraintsSetUp()
6823 POD: Y^H*Y = V*D*V^H, D = S^H*S -> U = Y*V*D^(-1/2) in PCBDDCConstraintsSetUp()
6824 -> When PETSC_USE_COMPLEX and PETSC_MISSING_LAPACK_GESVD are defined in PCBDDCConstraintsSetUp()
6827 -> This is due to a different computation of eigenvectors in LAPACKheev in PCBDDCConstraintsSetUp()
6828 -> The quality of the POD-computed basis will be the same */ in PCBDDCConstraintsSetUp()
6848 … while (j < temp_constraints && singular_vals[j] / singular_vals[temp_constraints - 1] < tol) j++; in PCBDDCConstraintsSetUp()
6849 total_counts = total_counts - j; in PCBDDCConstraintsSetUp()
6850 valid_constraints = temp_constraints - j; in PCBDDCConstraintsSetUp()
6864 for (k = 0; k < temp_constraints - j; k++) { in PCBDDCConstraintsSetUp()
6865 …_of_constraint + ii] = singular_vals[temp_constraints - 1 - k] * temp_basis[(temp_constraints - 1 … in PCBDDCConstraintsSetUp()
6885 while (j < k && singular_vals[k - j - 1] / singular_vals[0] < tol) j++; in PCBDDCConstraintsSetUp()
6886 valid_constraints = k - j; in PCBDDCConstraintsSetUp()
6887 total_counts = total_counts - temp_constraints + valid_constraints; in PCBDDCConstraintsSetUp()
6917 …PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, &o_nf, &ISForFaces, &o_ne, &ISForEdges… in PCBDDCConstraintsSetUp()
6919 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCConstraintsSetUp()
6923 …if (sub_schurs->is_vertices && pcbddc->use_vertices) PetscCall(ISGetLocalSize(sub_schurs->is_verti… in PCBDDCConstraintsSetUp()
6926 for (i = 0; i < sub_schurs->n_subs + n_vertices; i++) { in PCBDDCConstraintsSetUp()
6927 total_counts += pcbddc->adaptive_constraints_n[i]; in PCBDDCConstraintsSetUp()
6928 if (pcbddc->adaptive_constraints_n[i]) total_counts_cc++; in PCBDDCConstraintsSetUp()
6929 max_constraints = PetscMax(max_constraints, pcbddc->adaptive_constraints_n[i]); in PCBDDCConstraintsSetUp()
6931 constraints_idxs_ptr = pcbddc->adaptive_constraints_idxs_ptr; in PCBDDCConstraintsSetUp()
6932 constraints_data_ptr = pcbddc->adaptive_constraints_data_ptr; in PCBDDCConstraintsSetUp()
6933 constraints_idxs = pcbddc->adaptive_constraints_idxs; in PCBDDCConstraintsSetUp()
6934 constraints_data = pcbddc->adaptive_constraints_data; in PCBDDCConstraintsSetUp()
6935 /* constraints_n differs from pcbddc->adaptive_constraints_n */ in PCBDDCConstraintsSetUp()
6938 for (i = 0; i < sub_schurs->n_subs + n_vertices; i++) { in PCBDDCConstraintsSetUp()
6939 …if (pcbddc->adaptive_constraints_n[i]) constraints_n[total_counts_cc++] = pcbddc->adaptive_constra… in PCBDDCConstraintsSetUp()
6943 …f_constraint = PetscMax(max_size_of_constraint, constraints_idxs_ptr[i + 1] - constraints_idxs_ptr… in PCBDDCConstraintsSetUp()
6947 if (pcbddc->use_change_of_basis) { in PCBDDCConstraintsSetUp()
6948 for (i = 0; i < sub_schurs->n_subs; i++) { in PCBDDCConstraintsSetUp()
6949 …if (PetscBTLookup(sub_schurs->is_edge, i) || pcbddc->use_change_on_faces) PetscCall(PetscBTSet(cha… in PCBDDCConstraintsSetUp()
6953 pcbddc->local_primal_size = total_counts; in PCBDDCConstraintsSetUp()
6954 …PetscCall(PetscMalloc1(pcbddc->local_primal_size + pcbddc->benign_n, &pcbddc->primal_indices_local… in PCBDDCConstraintsSetUp()
6957 if (pcbddc->use_change_of_basis) { in PCBDDCConstraintsSetUp()
6958 …PetscCall(ISGlobalToLocalMappingApply(pcis->BtoNmap, IS_GTOLM_DROP, constraints_idxs_ptr[total_cou… in PCBDDCConstraintsSetUp()
6963 PetscCall(MatCreate(PETSC_COMM_SELF, &pcbddc->ConstraintMatrix)); in PCBDDCConstraintsSetUp()
6964 PetscCall(MatSetType(pcbddc->ConstraintMatrix, MATAIJ)); in PCBDDCConstraintsSetUp()
6965 …PetscCall(MatSetSizes(pcbddc->ConstraintMatrix, pcbddc->local_primal_size, pcis->n, pcbddc->local_… in PCBDDCConstraintsSetUp()
6969 qr_needed = pcbddc->use_qr_single; in PCBDDCConstraintsSetUp()
6972 pcbddc->local_primal_size_cc = 0; in PCBDDCConstraintsSetUp()
6974 size_of_constraint = constraints_idxs_ptr[i + 1] - constraints_idxs_ptr[i]; in PCBDDCConstraintsSetUp()
6975 if (size_of_constraint == 1 && pcbddc->mat_graph->custom_minimal_size) { in PCBDDCConstraintsSetUp()
6976 …pcbddc->primal_indices_local_idxs[total_primal_vertices++] = constraints_idxs[constraints_idxs_ptr… in PCBDDCConstraintsSetUp()
6977 pcbddc->local_primal_size_cc += 1; in PCBDDCConstraintsSetUp()
6979 …for (k = 0; k < constraints_n[i]; k++) pcbddc->primal_indices_local_idxs[total_primal_vertices++] … in PCBDDCConstraintsSetUp()
6980 pcbddc->local_primal_size_cc += constraints_n[i]; in PCBDDCConstraintsSetUp()
6981 if (constraints_n[i] > 1 || pcbddc->use_qr_single) { in PCBDDCConstraintsSetUp()
6986 pcbddc->local_primal_size_cc += 1; in PCBDDCConstraintsSetUp()
6990 pcbddc->n_vertices = total_primal_vertices; in PCBDDCConstraintsSetUp()
6992 PetscCall(PetscSortInt(total_primal_vertices, pcbddc->primal_indices_local_idxs)); in PCBDDCConstraintsSetUp()
6993 …oc2(pcbddc->local_primal_size_cc + pcbddc->benign_n, &pcbddc->local_primal_ref_node, pcbddc->local… in PCBDDCConstraintsSetUp()
6994 …PetscCall(PetscArraycpy(pcbddc->local_primal_ref_node, pcbddc->primal_indices_local_idxs, total_pr… in PCBDDCConstraintsSetUp()
6995 for (i = 0; i < total_primal_vertices; i++) pcbddc->local_primal_ref_mult[i] = 1; in PCBDDCConstraintsSetUp()
6999 PetscCall(PetscMalloc1(pcbddc->local_primal_size, &nnz)); in PCBDDCConstraintsSetUp()
7007 pcbddc->local_primal_ref_node[cum] = constraints_idxs[constraints_idxs_ptr[i]]; in PCBDDCConstraintsSetUp()
7008 pcbddc->local_primal_ref_mult[cum] = constraints_n[i]; in PCBDDCConstraintsSetUp()
7010 size_of_constraint = constraints_idxs_ptr[i + 1] - constraints_idxs_ptr[i]; in PCBDDCConstraintsSetUp()
7012 … pcbddc->primal_indices_local_idxs[total_counts++] = constraints_idxs[constraints_idxs_ptr[i] + k]; in PCBDDCConstraintsSetUp()
7018 PetscCall(MatSeqAIJSetPreallocation(pcbddc->ConstraintMatrix, 0, nnz)); in PCBDDCConstraintsSetUp()
7019 PetscCall(MatSetOption(pcbddc->ConstraintMatrix, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); in PCBDDCConstraintsSetUp()
7020 PetscCall(MatSetOption(pcbddc->ConstraintMatrix, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); in PCBDDCConstraintsSetUp()
7024 … i < total_primal_vertices; i++) PetscCall(MatSetValue(pcbddc->ConstraintMatrix, i, pcbddc->local_… in PCBDDCConstraintsSetUp()
7030 size_of_constraint = constraints_idxs_ptr[i + 1] - constraints_idxs_ptr[i]; in PCBDDCConstraintsSetUp()
7037 …PetscCall(MatSetValues(pcbddc->ConstraintMatrix, 1, &row, size_of_constraint, cols, vals, INSERT_V… in PCBDDCConstraintsSetUp()
7043 PetscCall(MatAssemblyBegin(pcbddc->ConstraintMatrix, MAT_FINAL_ASSEMBLY)); in PCBDDCConstraintsSetUp()
7044 PetscCall(MatAssemblyEnd(pcbddc->ConstraintMatrix, MAT_FINAL_ASSEMBLY)); in PCBDDCConstraintsSetUp()
7045 …PetscCall(MatViewFromOptions(pcbddc->ConstraintMatrix, (PetscObject)pc, "-pc_bddc_constraint_mat_v… in PCBDDCConstraintsSetUp()
7047 …/* Create matrix for change of basis. We don't need it in case pcbddc->use_change_of_basis is FALS… in PCBDDCConstraintsSetUp()
7048 if (pcbddc->use_change_of_basis) { in PCBDDCConstraintsSetUp()
7075 PetscCall(MatSetSizes(localChangeOfBasisMatrix, pcis->n, pcis->n, pcis->n, pcis->n)); in PCBDDCConstraintsSetUp()
7077 PetscCall(PetscMalloc1(pcis->n, &nnz)); in PCBDDCConstraintsSetUp()
7078 if (!pcbddc->benign_change || pcbddc->fake_change) { in PCBDDCConstraintsSetUp()
7079 for (i = 0; i < pcis->n; i++) nnz[i] = 1; in PCBDDCConstraintsSetUp()
7084 …PetscCall(MatGetRowIJ(pcbddc->benign_change, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, NULL, &flg_row)… in PCBDDCConstraintsSetUp()
7085 for (i = 0; i < n; i++) nnz[i] = ii[i + 1] - ii[i]; in PCBDDCConstraintsSetUp()
7086 …PetscCall(MatRestoreRowIJ(pcbddc->benign_change, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, NULL, &flg_… in PCBDDCConstraintsSetUp()
7090 size_of_constraint = constraints_idxs_ptr[i + 1] - constraints_idxs_ptr[i]; in PCBDDCConstraintsSetUp()
7103 if (!pcbddc->benign_change || pcbddc->fake_change) { in PCBDDCConstraintsSetUp()
7104 …for (i = 0; i < pcis->n; i++) PetscCall(MatSetValue(localChangeOfBasisMatrix, i, i, 1.0, INSERT_VA… in PCBDDCConstraintsSetUp()
7110 …PetscCall(MatGetRowIJ(pcbddc->benign_change, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &flg_row)); in PCBDDCConstraintsSetUp()
7111 PetscCall(MatSeqAIJGetArray(pcbddc->benign_change, &aa)); in PCBDDCConstraintsSetUp()
7112 …for (i = 0; i < n; i++) PetscCall(MatSetValues(localChangeOfBasisMatrix, 1, &i, ii[i + 1] - ii[i],… in PCBDDCConstraintsSetUp()
7113 PetscCall(MatSeqAIJRestoreArray(pcbddc->benign_change, &aa)); in PCBDDCConstraintsSetUp()
7114 …PetscCall(MatRestoreRowIJ(pcbddc->benign_change, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &flg_r… in PCBDDCConstraintsSetUp()
7117 if (pcbddc->dbg_flag) { in PCBDDCConstraintsSetUp()
7118 …scViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "---------------------------------------------… in PCBDDCConstraintsSetUp()
7119 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Checking change of basis computa… in PCBDDCConstraintsSetUp()
7125 Klawonn and Widlund, Dual-primal FETI-DP methods for linear elasticity, (see Sect 6.2.1) in PCBDDCConstraintsSetUp()
7129 …- By using the following block transformation if there is only a primal dof on the cc (and -pc_bdd… in PCBDDCConstraintsSetUp()
7134 | 0 ... 1 s_{n-1}/S | in PCBDDCConstraintsSetUp()
7135 | -s_1/s_n ... -s_{n-1}/s_n s_n/S | in PCBDDCConstraintsSetUp()
7141 - QR decomposition of constraints otherwise in PCBDDCConstraintsSetUp()
7152 lqr_work = -1; in PCBDDCConstraintsSetUp()
7157 lgqr_work = -1; in PCBDDCConstraintsSetUp()
7170 …if (pcbddc->dbg_flag) PetscCall(PetscMalloc1(max_size_of_constraint * (max_constraints + max_size_… in PCBDDCConstraintsSetUp()
7173 PetscCall(PetscBTCreate(pcis->n_B, &is_primal)); in PCBDDCConstraintsSetUp()
7175 …PetscCall(ISGlobalToLocalMappingApply(pcis->BtoNmap, IS_GTOLM_DROP, total_primal_vertices, pcbddc-… in PCBDDCConstraintsSetUp()
7182 … size_of_constraint = constraints_idxs_ptr[total_counts + 1] - constraints_idxs_ptr[total_counts]; in PCBDDCConstraintsSetUp()
7186 dual_dofs = size_of_constraint - primal_dofs; in PCBDDCConstraintsSetUp()
7188 …if (pcbddc->dbg_flag) PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Constraint… in PCBDDCConstraintsSetUp()
7193 …if (pcbddc->dbg_flag) PetscCall(PetscArraycpy(dbg_work, &constraints_data[constraints_data_ptr[tot… in PCBDDCConstraintsSetUp()
7206 /* explicitly compute R^-T */ in PCBDDCConstraintsSetUp()
7228 … /* first primal_dofs columns of Q need to be re-scaled in order to be unitary w.r.t constraints in PCBDDCConstraintsSetUp()
7261 if (pcbddc->dbg_flag) { in PCBDDCConstraintsSetUp()
7275 …ar(dbg_work[size_of_constraint * primal_dofs + jj * primal_dofs + ii]) > 1.e-12) valid_qr = PETSC_… in PCBDDCConstraintsSetUp()
7276 …g_work[size_of_constraint * primal_dofs + jj * primal_dofs + ii] - (PetscReal)1) > 1.e-12) valid_q… in PCBDDCConstraintsSetUp()
7280 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "\t-> wrong change of basis!\n")); in PCBDDCConstraintsSetUp()
7283 …j && PetscAbsScalar(dbg_work[size_of_constraint * primal_dofs + jj * primal_dofs + ii]) > 1.e-12) { in PCBDDCConstraintsSetUp()
7284 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "\tQr basis function %" PetscInt_… in PCBDDCConstraintsSetUp()
7286 …alar(dbg_work[size_of_constraint * primal_dofs + jj * primal_dofs + ii] - (PetscReal)1) > 1.e-12) { in PCBDDCConstraintsSetUp()
7287 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "\tQr basis function %" PetscInt_… in PCBDDCConstraintsSetUp()
7292 …PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "\t-> right change of basis!\n")); in PCBDDCConstraintsSetUp()
7312 …val = -constraints_data[constraints_data_ptr[total_counts] + k] / constraints_data[constraints_dat… in PCBDDCConstraintsSetUp()
7320 …if (pcbddc->dbg_flag) PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "\t-> using… in PCBDDCConstraintsSetUp()
7323 …if (pcbddc->dbg_flag) PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Constraint… in PCBDDCConstraintsSetUp()
7329 if (pcbddc->dbg_flag) PetscCall(PetscFree(dbg_work)); in PCBDDCConstraintsSetUp()
7341 if (!pcbddc->fake_change) { in PCBDDCConstraintsSetUp()
7344 PetscCall(VecGetSize(pcis->vec1_global, &global_size)); in PCBDDCConstraintsSetUp()
7345 PetscCall(VecGetLocalSize(pcis->vec1_global, &local_size)); in PCBDDCConstraintsSetUp()
7346 PetscCall(MatDuplicate(pc->pmat, MAT_DO_NOT_COPY_VALUES, &tmat)); in PCBDDCConstraintsSetUp()
7350 PetscCall(MatConvert(tmat, MATAIJ, MAT_INITIAL_MATRIX, &pcbddc->ChangeOfBasisMatrix)); in PCBDDCConstraintsSetUp()
7352 PetscCall(VecSet(pcis->vec1_global, 0.0)); in PCBDDCConstraintsSetUp()
7353 PetscCall(VecSet(pcis->vec1_N, 1.0)); in PCBDDCConstraintsSetUp()
7354 …PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERS… in PCBDDCConstraintsSetUp()
7355 …PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE)… in PCBDDCConstraintsSetUp()
7356 PetscCall(VecReciprocal(pcis->vec1_global)); in PCBDDCConstraintsSetUp()
7357 PetscCall(MatDiagonalScale(pcbddc->ChangeOfBasisMatrix, pcis->vec1_global, NULL)); in PCBDDCConstraintsSetUp()
7360 if (pcbddc->dbg_flag) { in PCBDDCConstraintsSetUp()
7364 PetscCall(VecDuplicate(pcis->vec1_global, &x)); in PCBDDCConstraintsSetUp()
7365 PetscCall(VecDuplicate(pcis->vec1_global, &x_change)); in PCBDDCConstraintsSetUp()
7367 PetscCall(VecCopy(x, pcis->vec1_global)); in PCBDDCConstraintsSetUp()
7368 PetscCall(VecScatterBegin(matis->rctx, x, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCConstraintsSetUp()
7369 PetscCall(VecScatterEnd(matis->rctx, x, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD)); in PCBDDCConstraintsSetUp()
7370 PetscCall(MatMult(localChangeOfBasisMatrix, pcis->vec1_N, pcis->vec2_N)); in PCBDDCConstraintsSetUp()
7371 PetscCall(VecScatterBegin(matis->rctx, pcis->vec2_N, x, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCConstraintsSetUp()
7372 PetscCall(VecScatterEnd(matis->rctx, pcis->vec2_N, x, INSERT_VALUES, SCATTER_REVERSE)); in PCBDDCConstraintsSetUp()
7373 PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcis->vec1_global, x_change)); in PCBDDCConstraintsSetUp()
7374 PetscCall(VecAXPY(x, -1.0, x_change)); in PCBDDCConstraintsSetUp()
7381 if (pcbddc->use_deluxe_scaling) { in PCBDDCConstraintsSetUp()
7382 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCConstraintsSetUp()
7384 …->use_change_of_basis || !pcbddc->adaptive_userdefined, PetscObjectComm((PetscObject)pc), PETSC_ER… in PCBDDCConstraintsSetUp()
7385 if (sub_schurs && sub_schurs->S_Ej_all) { in PCBDDCConstraintsSetUp()
7389 PetscCall(ISLocalToGlobalMappingApplyIS(pcis->BtoNmap, sub_schurs->is_Ej_all, &is_all_N)); in PCBDDCConstraintsSetUp()
7391 if (pcbddc->deluxe_zerorows) { in PCBDDCConstraintsSetUp()
7394 …PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pcbddc->n_vertices, pcbddc->local_primal_ref_node, PETS… in PCBDDCConstraintsSetUp()
7401 PetscCall(MatPtAP(sub_schurs->S_Ej_all, tmat, MAT_INITIAL_MATRIX, 1.0, &S_new)); in PCBDDCConstraintsSetUp()
7402 PetscCall(MatDestroy(&sub_schurs->S_Ej_all)); in PCBDDCConstraintsSetUp()
7404 if (pcbddc->deluxe_zerorows) { in PCBDDCConstraintsSetUp()
7412 PetscCall(ISGetIndices(sub_schurs->is_Ej_all, &idxs_all)); in PCBDDCConstraintsSetUp()
7413 PetscCall(VecGetArrayRead(pcis->D, &array)); in PCBDDCConstraintsSetUp()
7424 PetscCall(VecRestoreArrayRead(pcis->D, &array)); in PCBDDCConstraintsSetUp()
7425 PetscCall(ISRestoreIndices(sub_schurs->is_Ej_all, &idxs_all)); in PCBDDCConstraintsSetUp()
7428 sub_schurs->S_Ej_all = S_new; in PCBDDCConstraintsSetUp()
7430 if (sub_schurs->sum_S_Ej_all) { in PCBDDCConstraintsSetUp()
7431 PetscCall(MatPtAP(sub_schurs->sum_S_Ej_all, tmat, MAT_INITIAL_MATRIX, 1.0, &S_new)); in PCBDDCConstraintsSetUp()
7432 PetscCall(MatDestroy(&sub_schurs->sum_S_Ej_all)); in PCBDDCConstraintsSetUp()
7434 … if (pcbddc->deluxe_zerorows) PetscCall(MatZeroRowsColumnsIS(S_new, is_V_Sall, 1., NULL, NULL)); in PCBDDCConstraintsSetUp()
7435 sub_schurs->sum_S_Ej_all = S_new; in PCBDDCConstraintsSetUp()
7442 if (sub_schurs && sub_schurs->change) { in PCBDDCConstraintsSetUp()
7445 for (i = 0; i < sub_schurs->n_subs; i++) PetscCall(KSPDestroy(&sub_schurs->change[i])); in PCBDDCConstraintsSetUp()
7446 PetscCall(PetscFree(sub_schurs->change)); in PCBDDCConstraintsSetUp()
7449 if (pcbddc->switch_static) { /* need to save the local change */ in PCBDDCConstraintsSetUp()
7450 pcbddc->switch_static_change = localChangeOfBasisMatrix; in PCBDDCConstraintsSetUp()
7455 pcbddc->change_interior = pcbddc->benign_have_null; in PCBDDCConstraintsSetUp()
7457 PetscCall(MatDestroy(&pcbddc->ConstraintMatrix)); in PCBDDCConstraintsSetUp()
7458 pcbddc->ConstraintMatrix = localChangeOfBasisMatrix; in PCBDDCConstraintsSetUp()
7459 pcbddc->use_qr_single = qr_needed; in PCBDDCConstraintsSetUp()
7461 } else if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->benign_saddle_point) { in PCBDDCConstraintsSetUp()
7462 if (!pcbddc->benign_have_null && pcbddc->user_ChangeOfBasisMatrix) { in PCBDDCConstraintsSetUp()
7463 PetscCall(PetscObjectReference((PetscObject)pcbddc->user_ChangeOfBasisMatrix)); in PCBDDCConstraintsSetUp()
7464 pcbddc->ChangeOfBasisMatrix = pcbddc->user_ChangeOfBasisMatrix; in PCBDDCConstraintsSetUp()
7467 if (pcbddc->benign_have_null) { in PCBDDCConstraintsSetUp()
7470 pcbddc->change_interior = PETSC_TRUE; in PCBDDCConstraintsSetUp()
7471 PetscCall(VecCopy(matis->counter, pcis->vec1_N)); in PCBDDCConstraintsSetUp()
7472 PetscCall(VecReciprocal(pcis->vec1_N)); in PCBDDCConstraintsSetUp()
7473 PetscCall(MatDuplicate(pc->pmat, MAT_DO_NOT_COPY_VALUES, &benign_global)); in PCBDDCConstraintsSetUp()
7474 if (pcbddc->benign_change) { in PCBDDCConstraintsSetUp()
7475 PetscCall(MatDuplicate(pcbddc->benign_change, MAT_COPY_VALUES, &M)); in PCBDDCConstraintsSetUp()
7476 PetscCall(MatDiagonalScale(M, pcis->vec1_N, NULL)); in PCBDDCConstraintsSetUp()
7478 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, pcis->n, pcis->n, 1, NULL, &M)); in PCBDDCConstraintsSetUp()
7479 PetscCall(MatDiagonalSet(M, pcis->vec1_N, INSERT_VALUES)); in PCBDDCConstraintsSetUp()
7486 if (pcbddc->user_ChangeOfBasisMatrix) { in PCBDDCConstraintsSetUp()
7487 …PetscCall(MatMatMult(pcbddc->user_ChangeOfBasisMatrix, benign_global, MAT_INITIAL_MATRIX, PETSC_DE… in PCBDDCConstraintsSetUp()
7489 } else if (pcbddc->benign_have_null) { in PCBDDCConstraintsSetUp()
7490 pcbddc->ChangeOfBasisMatrix = benign_global; in PCBDDCConstraintsSetUp()
7493 if (pcbddc->switch_static && pcbddc->ChangeOfBasisMatrix) { /* need to save the local change */ in PCBDDCConstraintsSetUp()
7497 PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs)); in PCBDDCConstraintsSetUp()
7498 …PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), pcis->n, gidxs, PETSC_COPY_VALUES, &is… in PCBDDCConstraintsSetUp()
7499 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs)); in PCBDDCConstraintsSetUp()
7500 …PetscCall(MatCreateSubMatrixUnsorted(pcbddc->ChangeOfBasisMatrix, is_global, is_global, &pcbddc->s… in PCBDDCConstraintsSetUp()
7504 …if (!pcbddc->fake_change && pcbddc->ChangeOfBasisMatrix && !pcbddc->work_change) PetscCall(VecDupl… in PCBDDCConstraintsSetUp()
7506 if (!pcbddc->fake_change) { in PCBDDCConstraintsSetUp()
7508 for (i = 0; i < pcbddc->benign_n; i++) { in PCBDDCConstraintsSetUp()
7509 pcbddc->local_primal_ref_node[pcbddc->local_primal_size_cc] = pcbddc->benign_p0_lidx[i]; in PCBDDCConstraintsSetUp()
7510 pcbddc->primal_indices_local_idxs[pcbddc->local_primal_size] = pcbddc->benign_p0_lidx[i]; in PCBDDCConstraintsSetUp()
7511 pcbddc->local_primal_ref_mult[pcbddc->local_primal_size_cc] = 1; in PCBDDCConstraintsSetUp()
7512 pcbddc->local_primal_size_cc++; in PCBDDCConstraintsSetUp()
7513 pcbddc->local_primal_size++; in PCBDDCConstraintsSetUp()
7517 pcbddc->new_primal_space_local = PETSC_TRUE; in PCBDDCConstraintsSetUp()
7518 if (olocal_primal_size == pcbddc->local_primal_size) { in PCBDDCConstraintsSetUp()
7519 …PetscCall(PetscArraycmp(pcbddc->local_primal_ref_node, olocal_primal_ref_node, olocal_primal_size_… in PCBDDCConstraintsSetUp()
7520 pcbddc->new_primal_space_local = (PetscBool)(!pcbddc->new_primal_space_local); in PCBDDCConstraintsSetUp()
7521 if (!pcbddc->new_primal_space_local) { in PCBDDCConstraintsSetUp()
7522 …PetscCall(PetscArraycmp(pcbddc->local_primal_ref_mult, olocal_primal_ref_mult, olocal_primal_size_… in PCBDDCConstraintsSetUp()
7523 pcbddc->new_primal_space_local = (PetscBool)(!pcbddc->new_primal_space_local); in PCBDDCConstraintsSetUp()
7527 …PetscCallMPI(MPIU_Allreduce(&pcbddc->new_primal_space_local, &pcbddc->new_primal_space, 1, MPI_C_B… in PCBDDCConstraintsSetUp()
7532 if (pcbddc->dbg_flag) PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCConstraintsSetUp()
7537 if (!pcbddc->adaptive_selection) { in PCBDDCConstraintsSetUp()
7541 …ddc->adaptive_constraints_n, pcbddc->adaptive_constraints_idxs_ptr, pcbddc->adaptive_constraints_d… in PCBDDCConstraintsSetUp()
7551 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCAnalyzeInterface()
7552 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCAnalyzeInterface()
7557 if (pcbddc->recompute_topography) { in PCBDDCAnalyzeInterface()
7558 pcbddc->graphanalyzed = PETSC_FALSE; in PCBDDCAnalyzeInterface()
7560 PetscCall(PCBDDCGraphReset(pcbddc->mat_graph)); in PCBDDCAnalyzeInterface()
7562 PetscCall(MatGetSize(pc->pmat, &N, NULL)); in PCBDDCAnalyzeInterface()
7563 PetscCall(MatISGetLocalToGlobalMapping(pc->pmat, &map, NULL)); in PCBDDCAnalyzeInterface()
7564 PetscCall(PCBDDCGraphInit(pcbddc->mat_graph, map, N, pcbddc->graphmaxcount)); in PCBDDCAnalyzeInterface()
7566 …if (pcbddc->user_primal_vertices_local && !pcbddc->user_primal_vertices) PetscCall(PCBDDCConsisten… in PCBDDCAnalyzeInterface()
7568 …->mat_graph->nvtxs_csr || pcbddc->mat_graph->nvtxs_csr == pcbddc->mat_graph->nvtxs, PETSC_COMM_SEL… in PCBDDCAnalyzeInterface()
7569 pcbddc->mat_graph->nvtxs); in PCBDDCAnalyzeInterface()
7572 if (!pcbddc->mat_graph->xadj && pcbddc->use_local_adj) { in PCBDDCAnalyzeInterface()
7578 PetscCall(PetscObjectReference((PetscObject)matis->A)); in PCBDDCAnalyzeInterface()
7579 A = matis->A; in PCBDDCAnalyzeInterface()
7580 for (PetscInt i = 0; i < pcbddc->local_adj_square; i++) { in PCBDDCAnalyzeInterface()
7590 AtA->assembled = PETSC_TRUE; in PCBDDCAnalyzeInterface()
7597 pcbddc->computed_rowadj = PETSC_TRUE; in PCBDDCAnalyzeInterface()
7603 if (pcbddc->dbg_flag) PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCAnalyzeInterface()
7605 if (pcbddc->mat_graph->cdim && !pcbddc->mat_graph->cloc) { in PCBDDCAnalyzeInterface()
7612 …->mat_graph->cnloc == pc->pmat->rmap->n, PETSC_COMM_SELF, PETSC_ERR_USER, "Invalid number of local… in PCBDDCAnalyzeInterface()
7613 PetscCall(MatGetLocalSize(matis->A, &n, NULL)); in PCBDDCAnalyzeInterface()
7614 PetscCall(PetscMalloc1(pcbddc->mat_graph->cdim * n, &lcoords)); in PCBDDCAnalyzeInterface()
7615 PetscCall(PetscMPIIntCast(pcbddc->mat_graph->cdim, &cdimi)); in PCBDDCAnalyzeInterface()
7618 …PetscCall(PetscSFBcastBegin(matis->sf, dimrealtype, pcbddc->mat_graph->coords, lcoords, MPI_REPLAC… in PCBDDCAnalyzeInterface()
7619 …PetscCall(PetscSFBcastEnd(matis->sf, dimrealtype, pcbddc->mat_graph->coords, lcoords, MPI_REPLACE)… in PCBDDCAnalyzeInterface()
7621 PetscCall(PetscFree(pcbddc->mat_graph->coords)); in PCBDDCAnalyzeInterface()
7623 pcbddc->mat_graph->coords = lcoords; in PCBDDCAnalyzeInterface()
7624 pcbddc->mat_graph->cloc = PETSC_TRUE; in PCBDDCAnalyzeInterface()
7625 pcbddc->mat_graph->cnloc = n; in PCBDDCAnalyzeInterface()
7627 …->mat_graph->cnloc || pcbddc->mat_graph->cnloc == pcbddc->mat_graph->nvtxs, PETSC_COMM_SELF, PETSC… in PCBDDCAnalyzeInterface()
7628 pcbddc->mat_graph->nvtxs); in PCBDDCAnalyzeInterface()
7629 …pcbddc->mat_graph->active_coords = (PetscBool)(pcbddc->corner_selection && pcbddc->mat_graph->cdim… in PCBDDCAnalyzeInterface()
7632 if (pcbddc->n_local_subs) { in PCBDDCAnalyzeInterface()
7635 PetscCall(MatGetLocalSize(matis->A, &n, NULL)); in PCBDDCAnalyzeInterface()
7637 for (i = 0; i < n; i++) local_subs[i] = pcbddc->n_local_subs; in PCBDDCAnalyzeInterface()
7638 for (i = 0; i < pcbddc->n_local_subs; i++) { in PCBDDCAnalyzeInterface()
7642 PetscCall(ISGetLocalSize(pcbddc->local_subs[i], &nl)); in PCBDDCAnalyzeInterface()
7643 PetscCall(ISGetIndices(pcbddc->local_subs[i], &idxs)); in PCBDDCAnalyzeInterface()
7645 PetscCall(ISRestoreIndices(pcbddc->local_subs[i], &idxs)); in PCBDDCAnalyzeInterface()
7648 pcbddc->mat_graph->n_local_subs = totn + 1; in PCBDDCAnalyzeInterface()
7649 pcbddc->mat_graph->local_subs = local_subs; in PCBDDCAnalyzeInterface()
7653 …->mat_graph, pcbddc->vertex_size, pcbddc->NeumannBoundariesLocal, pcbddc->DirichletBoundariesLocal… in PCBDDCAnalyzeInterface()
7656 if (!pcbddc->graphanalyzed) { in PCBDDCAnalyzeInterface()
7658 PetscCall(PCBDDCGraphComputeConnectedComponents(pcbddc->mat_graph)); in PCBDDCAnalyzeInterface()
7659 pcbddc->graphanalyzed = PETSC_TRUE; in PCBDDCAnalyzeInterface()
7660 pcbddc->corner_selected = pcbddc->corner_selection; in PCBDDCAnalyzeInterface()
7662 if (rcsr) pcbddc->mat_graph->nvtxs_csr = 0; in PCBDDCAnalyzeInterface()
7686 for (j = 0; j < i; j++) alphas[j] = PetscConj(-alphas[j]); in PCBDDCOrthonormalizeVecs()
7744 void_procs = size - active_procs; in PCBDDCMatISGetSubassemblingPattern()
7745 /* get ranks of non-active processes in mat communicator */ in PCBDDCMatISGetSubassemblingPattern()
7755 /* force n_subdomains to be not greater that the number of non-active processes */ in PCBDDCMatISGetSubassemblingPattern()
7759 …/* number of subdomains requested greater than active processes or matrix size -> just shift the m… in PCBDDCMatISGetSubassemblingPattern()
7760 number of subdomains requested 1 -> send to rank-0 or first candidate in voids */ in PCBDDCMatISGetSubassemblingPattern()
7770 if (procs_candidates) { /* shift the pattern on non-active candidates (if any) */ in PCBDDCMatISGetSubassemblingPattern()
7787 …PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_partitioning_use_vwgt", &us… in PCBDDCMatISGetSubassemblingPattern()
7788 …PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)A)->prefix, "-mat_is_partitioning_threshold", &th… in PCBDDCMatISGetSubassemblingPattern()
7798 xadj[1] = PetscMax(n_neighs - 1, 0); in PCBDDCMatISGetSubassemblingPattern()
7879 ncols = xadj[1] - xadj[0]; in PCBDDCMatISGetSubassemblingPattern()
7937 if (procs_candidates) { /* shift the pattern on non-active candidates (if any) */ in PCBDDCMatISGetSubassemblingPattern()
7951 PetscCall(PetscMalloc1(rend - rstart, &reqs)); in PCBDDCMatISGetSubassemblingPattern()
7952 … irend; i++) PetscCallMPI(MPIU_Isend(is_indices + i - rstart, 1, MPIU_INT, i, tag, subcomm, &reqs[… in PCBDDCMatISGetSubassemblingPattern()
7954 PetscCallMPI(MPI_Waitall(irend - irstart, reqs, MPI_STATUSES_IGNORE)); in PCBDDCMatISGetSubassemblingPattern()
7956 if (procs_candidates) { /* shift the pattern on non-active candidates (if any) */ in PCBDDCMatISGetSubassemblingPattern()
8089 subcommsize = size - subcommsize; in PCBDDCMatISSubassemble()
8128 - MatType_PRIVATE type in PCBDDCMatISSubassemble()
8129 - PetscInt size_of_l2gmap in PCBDDCMatISSubassemble()
8130 - PetscInt global_row_indices[size_of_l2gmap] in PCBDDCMatISSubassemble()
8131 - PetscInt all_other_info_which_is_needed_to_compute_preallocation_and_set_values in PCBDDCMatISSubassemble()
8224 …PetscCallMPI(MPIU_Irecv(ptr_vecs, olengths_idxs[i] - 2, MPIU_SCALAR, onodes[i], tag_vecs, comm, &r… in PCBDDCMatISSubassemble()
8225 ptr_vecs += olengths_idxs[i] - 2; in PCBDDCMatISSubassemble()
8235 …PetscCallMPI(MPIU_Isend(send_buffer_vecs, ilengths_idxs[source_dest] - 2, MPIU_SCALAR, source_dest… in PCBDDCMatISSubassemble()
8319 /* restore attributes -> type of incoming data and its size */ in PCBDDCMatISSubassemble()
8349 for (i = 0; i < new_local_rows; i++) new_local_nnz[i] = PetscMax(new_local_nnz[i] - i, 0); in PCBDDCMatISSubassemble()
8390 PetscCall(VecScale(lvec,-1.0)); in PCBDDCMatISSubassemble()
8410 count_is[j] += plen; /* increment counting of buffer for j-th IS */ in PCBDDCMatISSubassemble()
8417 …(i = 1; i < nis; i++) temp_idxs[i] = PetscSafePointerPlusOffset(temp_idxs[i - 1], count_is[i - 1]); in PCBDDCMatISSubassemble()
8424 count_is[j] += plen; /* increment starting point of buffer for j-th IS */ in PCBDDCMatISSubassemble()
8470 ptr_vals += olengths_idxs[i] - 2; in PCBDDCMatISSubassemble()
8514 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCSetUpCoarseSolver()
8515 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCSetUpCoarseSolver()
8516 PCBDDCGraph graph = pcbddc->mat_graph; in PCBDDCSetUpCoarseSolver()
8523 PetscInt i, im_active = -1, active_procs = -1; in PCBDDCSetUpCoarseSolver()
8530 PetscBool coarse_reuse, multi_element = graph->multi_element; in PCBDDCSetUpCoarseSolver()
8539 PetscCall(PetscLogEventBegin(PC_BDDC_CoarseSetUp[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpCoarseSolver()
8541 …if (pcbddc->new_primal_space || pcbddc->coarse_size == -1) { /* a new primal space is present or i… in PCBDDCSetUpCoarseSolver()
8545 pcbddc->new_primal_space = PETSC_TRUE; in PCBDDCSetUpCoarseSolver()
8546 ocoarse_size = pcbddc->coarse_size; in PCBDDCSetUpCoarseSolver()
8547 PetscCall(PetscFree(pcbddc->global_primal_indices)); in PCBDDCSetUpCoarseSolver()
8548 … PetscCall(PCBDDCComputePrimalNumbering(pc, &pcbddc->coarse_size, &pcbddc->global_primal_indices)); in PCBDDCSetUpCoarseSolver()
8550 if (pcbddc->coarse_ksp) { /* coarse ksp has already been created */ in PCBDDCSetUpCoarseSolver()
8552 if (ocoarse_size != pcbddc->coarse_size || pcbddc->adaptive_selection) { in PCBDDCSetUpCoarseSolver()
8553 PetscCall(KSPReset(pcbddc->coarse_ksp)); in PCBDDCSetUpCoarseSolver()
8562 …if (!coarse_reuse || pcbddc->recompute_topography) PetscCall(ISDestroy(&pcbddc->coarse_subassembli… in PCBDDCSetUpCoarseSolver()
8566 if (coarse_reuse && pcbddc->coarse_ksp) { in PCBDDCSetUpCoarseSolver()
8567 PetscCall(KSPGetOperators(pcbddc->coarse_ksp, &coarse_mat, NULL)); in PCBDDCSetUpCoarseSolver()
8576 …PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), pcbddc->local_primal_size, pcbddc->glo… in PCBDDCSetUpCoarseSolver()
8582 …PetscCall(MatSetSizes(t_coarse_mat_is, PETSC_DECIDE, PETSC_DECIDE, pcbddc->coarse_size, pcbddc->co… in PCBDDCSetUpCoarseSolver()
8588 PetscCall(MatViewFromOptions(t_coarse_mat_is, (PetscObject)pc, "-pc_bddc_coarse_mat_is_view")); in PCBDDCSetUpCoarseSolver()
8591 im_active = !!pcis->n; in PCBDDCSetUpCoarseSolver()
8601 coarse_eqs_per_proc = PetscMin(PetscMax(pcbddc->coarse_size, 1), pcbddc->coarse_eqs_per_proc); in PCBDDCSetUpCoarseSolver()
8602 if (coarse_eqs_per_proc < 0 || size == 1) coarse_eqs_per_proc = PetscMax(pcbddc->coarse_size, 1); in PCBDDCSetUpCoarseSolver()
8603 if (pcbddc->current_level < pcbddc->max_levels) multilevel_requested = PETSC_TRUE; in PCBDDCSetUpCoarseSolver()
8604 if (pcbddc->coarse_size <= pcbddc->coarse_eqs_limit) multilevel_requested = PETSC_FALSE; in PCBDDCSetUpCoarseSolver()
8605 coarsening_ratio = multi_element ? 1 : pcbddc->coarsening_ratio; in PCBDDCSetUpCoarseSolver()
8611 …ncoarse = pcbddc->coarse_size / coarse_eqs_per_proc + !!(pcbddc->coarse_size % coarse_eqs_per_p… in PCBDDCSetUpCoarseSolver()
8615 …if (!pcbddc->coarse_size || (size == 1 && !multi_element)) multilevel_allowed = multilevel_request… in PCBDDCSetUpCoarseSolver()
8617 if (!pcbddc->coarse_subassembling) { in PCBDDCSetUpCoarseSolver()
8620 …PetscCall(PCBDDCMatISGetSubassemblingPattern(pc->pmat, &ncoarse, pcbddc->coarse_adj_red, &pcbddc->… in PCBDDCSetUpCoarseSolver()
8622 …CMatISGetSubassemblingPattern(t_coarse_mat_is, &ncoarse, pcbddc->coarse_adj_red, &pcbddc->coarse_s… in PCBDDCSetUpCoarseSolver()
8629 …PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), 1, rank, 1, &pcbddc->coarse_subassembli… in PCBDDCSetUpCoarseSolver()
8630 … PetscCall(PetscObjectSetName((PetscObject)pcbddc->coarse_subassembling, "default subassembling")); in PCBDDCSetUpCoarseSolver()
8634 if (pcbddc->coarse_ksp) psum = 1; in PCBDDCSetUpCoarseSolver()
8647 …if (pcbddc->dbg_flag && multilevel_allowed) PetscCall(ISView(pcbddc->coarse_subassembling, pcbddc-… in PCBDDCSetUpCoarseSolver()
8649 nedcfield = -1; in PCBDDCSetUpCoarseSolver()
8651 …llowed && !coarse_reuse && (pcbddc->n_ISForDofsLocal || pcbddc->NeumannBoundariesLocal || pcbddc->… in PCBDDCSetUpCoarseSolver()
8657 …PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1, pcbddc->local_primal_size, pcbddc->prim… in PCBDDCSetUpCoarseSolver()
8659 PetscCall(PetscMalloc1(pcbddc->local_primal_size, &tidxs)); in PCBDDCSetUpCoarseSolver()
8660 PetscCall(PetscMalloc1(pcbddc->local_primal_size, &tidxs2)); in PCBDDCSetUpCoarseSolver()
8662 nisdofs = pcbddc->n_ISForDofsLocal; in PCBDDCSetUpCoarseSolver()
8663 if (pcbddc->nedclocal) { in PCBDDCSetUpCoarseSolver()
8664 if (pcbddc->nedfield > -1) { in PCBDDCSetUpCoarseSolver()
8665 nedcfield = pcbddc->nedfield; in PCBDDCSetUpCoarseSolver()
8672 nisneu = !!pcbddc->NeumannBoundariesLocal; in PCBDDCSetUpCoarseSolver()
8678 /* PetscCall(ISView(pcbddc->ISForDofsLocal[i],0)); */ in PCBDDCSetUpCoarseSolver()
8680 PetscCall(ISGetLocalSize(pcbddc->ISForDofsLocal[i], &tsize)); in PCBDDCSetUpCoarseSolver()
8681 PetscCall(ISGetIndices(pcbddc->ISForDofsLocal[i], &idxs)); in PCBDDCSetUpCoarseSolver()
8683 PetscCall(ISRestoreIndices(pcbddc->ISForDofsLocal[i], &idxs)); in PCBDDCSetUpCoarseSolver()
8685 PetscCall(ISGetLocalSize(pcbddc->nedclocal, &tsize)); in PCBDDCSetUpCoarseSolver()
8686 PetscCall(ISGetIndices(pcbddc->nedclocal, &idxs)); in PCBDDCSetUpCoarseSolver()
8689 PetscCall(ISRestoreIndices(pcbddc->nedclocal, &idxs)); in PCBDDCSetUpCoarseSolver()
8696 if (pcbddc->NeumannBoundariesLocal) { in PCBDDCSetUpCoarseSolver()
8697 /* PetscCall(ISView(pcbddc->NeumannBoundariesLocal,0)); */ in PCBDDCSetUpCoarseSolver()
8698 PetscCall(ISGetLocalSize(pcbddc->NeumannBoundariesLocal, &tsize)); in PCBDDCSetUpCoarseSolver()
8699 PetscCall(ISGetIndices(pcbddc->NeumannBoundariesLocal, &idxs)); in PCBDDCSetUpCoarseSolver()
8701 PetscCall(ISRestoreIndices(pcbddc->NeumannBoundariesLocal, &idxs)); in PCBDDCSetUpCoarseSolver()
8707 if (pcbddc->corner_selected) { in PCBDDCSetUpCoarseSolver()
8708 PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &corners)); in PCBDDCSetUpCoarseSolver()
8714 … PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &corners)); in PCBDDCSetUpCoarseSolver()
8739 if (pcbddc->benign_have_null) { /* propagate no-net-flux quadrature to coarser level */ in PCBDDCSetUpCoarseSolver()
8741 PetscCall(VecSetSizes(vp[0], pcbddc->local_primal_size, PETSC_DECIDE)); in PCBDDCSetUpCoarseSolver()
8745 if (pcbddc->divudotp) { in PCBDDCSetUpCoarseSolver()
8751 PetscCall(MatISGetLocalMat(pcbddc->divudotp, &loc_divudotp)); in PCBDDCSetUpCoarseSolver()
8754 … PetscCall(MatCreateSubMatrix(loc_divudotp, dummy, pcis->is_B_local, MAT_INITIAL_MATRIX, &B)); in PCBDDCSetUpCoarseSolver()
8761 PetscCall(VecPlaceArray(pcbddc->vec1_P, array)); in PCBDDCSetUpCoarseSolver()
8762 PetscCall(MatMultTranspose(pcbddc->coarse_phi_B, v, pcbddc->vec1_P)); in PCBDDCSetUpCoarseSolver()
8763 PetscCall(VecResetArray(pcbddc->vec1_P)); in PCBDDCSetUpCoarseSolver()
8777 …PetscCall(PCBDDCMatISSubassemble(t_coarse_mat_is, pcbddc->coarse_subassembling, 0, restr, full_res… in PCBDDCSetUpCoarseSolver()
8779 …PetscCall(PCBDDCMatISSubassemble(t_coarse_mat_is, pcbddc->coarse_subassembling, 0, restr, full_res… in PCBDDCSetUpCoarseSolver()
8800 …PetscCall(PetscStrcmp(((PetscObject)pcbddc->coarse_subassembling)->name, "default subassembling", … in PCBDDCSetUpCoarseSolver()
8801 …if (!default_sub) PetscCall(PCBDDCMatISSubassemble(t_coarse_mat_is, pcbddc->coarse_subassembling, … in PCBDDCSetUpCoarseSolver()
8825 PetscCall(VecDestroy(&pcbddc->coarse_vec)); in PCBDDCSetUpCoarseSolver()
8831 PetscCall(VecCreate(PetscObjectComm((PetscObject)pc), &pcbddc->coarse_vec)); in PCBDDCSetUpCoarseSolver()
8832 PetscCall(VecSetSizes(pcbddc->coarse_vec, lrows, PETSC_DECIDE)); in PCBDDCSetUpCoarseSolver()
8833 … PetscCall(VecSetType(pcbddc->coarse_vec, coarse_mat ? coarse_mat->defaultvectype : VECSTANDARD)); in PCBDDCSetUpCoarseSolver()
8834 PetscCall(VecScatterDestroy(&pcbddc->coarse_loc_to_glob)); in PCBDDCSetUpCoarseSolver()
8835 …PetscCall(VecScatterCreate(pcbddc->vec1_P, NULL, pcbddc->coarse_vec, coarse_is, &pcbddc->coarse_lo… in PCBDDCSetUpCoarseSolver()
8849 if (pcbddc->dbg_flag) { in PCBDDCSetUpCoarseSolver()
8851 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "--------------------------------------------… in PCBDDCSetUpCoarseSolver()
8853 …->dbg_viewer, "Not enough active processes on level %" PetscInt_FMT " (active processes %" PetscIn… in PCBDDCSetUpCoarseSolver()
8854 } else if (pcbddc->max_levels) { in PCBDDCSetUpCoarseSolver()
8855 …ll(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "Maximum number of requested levels reached (%" Pets… in PCBDDCSetUpCoarseSolver()
8857 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpCoarseSolver()
8863 if (pcbddc->nedcG && multilevel_allowed) { in PCBDDCSetUpCoarseSolver()
8870 PetscCall(MatMPIAIJRestrict(pcbddc->nedcG, ccomm, &coarseG)); in PCBDDCSetUpCoarseSolver()
8879 if (pcbddc->dbg_flag) { in PCBDDCSetUpCoarseSolver()
8881 PetscCall(PetscViewerASCIIAddTab(dbg_viewer, 2 * pcbddc->current_level)); in PCBDDCSetUpCoarseSolver()
8883 if (!pcbddc->coarse_ksp) { in PCBDDCSetUpCoarseSolver()
8887 PetscCall(KSPCreate(PetscObjectComm((PetscObject)coarse_mat), &pcbddc->coarse_ksp)); in PCBDDCSetUpCoarseSolver()
8888 PetscCall(KSPSetNestLevel(pcbddc->coarse_ksp, pc->kspnestlevel)); in PCBDDCSetUpCoarseSolver()
8889 PetscCall(KSPSetErrorIfNotConverged(pcbddc->coarse_ksp, pc->erroriffailure)); in PCBDDCSetUpCoarseSolver()
8890 PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcbddc->coarse_ksp, (PetscObject)pc, 1)); in PCBDDCSetUpCoarseSolver()
8891 … PetscCall(KSPSetTolerances(pcbddc->coarse_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 1)); in PCBDDCSetUpCoarseSolver()
8892 PetscCall(KSPSetOperators(pcbddc->coarse_ksp, coarse_mat, coarse_mat)); in PCBDDCSetUpCoarseSolver()
8893 PetscCall(KSPSetType(pcbddc->coarse_ksp, coarse_ksp_type)); in PCBDDCSetUpCoarseSolver()
8894 PetscCall(KSPSetNormType(pcbddc->coarse_ksp, KSP_NORM_NONE)); in PCBDDCSetUpCoarseSolver()
8895 PetscCall(KSPGetPC(pcbddc->coarse_ksp, &pc_temp)); in PCBDDCSetUpCoarseSolver()
8901 if (!pcbddc->current_level) { in PCBDDCSetUpCoarseSolver()
8902 PetscCall(PetscStrncpy(prefix, ((PetscObject)pc)->prefix, sizeof(prefix))); in PCBDDCSetUpCoarseSolver()
8905 PetscCall(PetscStrlen(((PetscObject)pc)->prefix, &len)); in PCBDDCSetUpCoarseSolver()
8906 if (pcbddc->current_level > 1) len -= 3; /* remove "lX_" with X level number */ in PCBDDCSetUpCoarseSolver()
8907 if (pcbddc->current_level > 10) len -= 1; /* remove another char from level number */ in PCBDDCSetUpCoarseSolver()
8909 PetscCall(PetscStrncpy(prefix, ((PetscObject)pc)->prefix, len + 1)); in PCBDDCSetUpCoarseSolver()
8910 …PetscCall(PetscSNPrintf(str_level, sizeof(str_level), "l%" PetscInt_FMT "_", pcbddc->current_level… in PCBDDCSetUpCoarseSolver()
8913 PetscCall(KSPSetOptionsPrefix(pcbddc->coarse_ksp, prefix)); in PCBDDCSetUpCoarseSolver()
8915 PetscCall(PCBDDCSetLevel(pc_temp, pcbddc->current_level + 1)); in PCBDDCSetUpCoarseSolver()
8916 PetscCall(PCBDDCSetCoarseningRatio(pc_temp, pcbddc->coarsening_ratio)); in PCBDDCSetUpCoarseSolver()
8917 PetscCall(PCBDDCSetLevels(pc_temp, pcbddc->max_levels)); in PCBDDCSetUpCoarseSolver()
8919 PetscCall(KSPSetFromOptions(pcbddc->coarse_ksp)); in PCBDDCSetUpCoarseSolver()
8921 PetscCall(KSPGetPC(pcbddc->coarse_ksp, &pc_temp)); in PCBDDCSetUpCoarseSolver()
8924 …PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)pc_temp)->prefix, "-pc_type_forced", &force, NUL… in PCBDDCSetUpCoarseSolver()
8930 PetscCall(PCBDDCSetLevel(pc_temp, pcbddc->current_level + 1)); in PCBDDCSetUpCoarseSolver()
8931 PetscCall(PCBDDCSetCoarseningRatio(pc_temp, pcbddc->coarsening_ratio)); in PCBDDCSetUpCoarseSolver()
8932 PetscCall(PCBDDCSetLevels(pc_temp, pcbddc->max_levels)); in PCBDDCSetUpCoarseSolver()
8933 if (pc_temp->ops->setfromoptions) { /* need to setfromoptions again, skipping the pc_type */ in PCBDDCSetUpCoarseSolver()
8935 PetscCall((*pc_temp->ops->setfromoptions)(pc_temp, PetscOptionsObject)); in PCBDDCSetUpCoarseSolver()
8938 pc_temp->setfromoptionscalled++; in PCBDDCSetUpCoarseSolver()
8943 PetscCall(KSPGetPC(pcbddc->coarse_ksp, &pc_temp)); in PCBDDCSetUpCoarseSolver()
8953 PetscCall(PCBDDCSetPrimalVerticesIS(pc_temp, isarray[nis - 1])); in PCBDDCSetUpCoarseSolver()
8954 PetscCall(ISDestroy(&isarray[nis - 1])); in PCBDDCSetUpCoarseSolver()
8961 /* multilevel can only be requested via -pc_bddc_levels or PCBDDCSetLevels */ in PCBDDCSetUpCoarseSolver()
8965 …PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)pc_temp)->prefix, "-pc_type_forced", &force, NUL… in PCBDDCSetUpCoarseSolver()
8980 PC_BDDC *pcbddc_coarse = (PC_BDDC *)pc_temp->data; in PCBDDCSetUpCoarseSolver()
8982 pcbddc_coarse->detect_disconnected = PETSC_TRUE; in PCBDDCSetUpCoarseSolver()
8983 pcbddc_coarse->coarse_eqs_per_proc = pcbddc->coarse_eqs_per_proc; in PCBDDCSetUpCoarseSolver()
8984 pcbddc_coarse->coarse_eqs_limit = pcbddc->coarse_eqs_limit; in PCBDDCSetUpCoarseSolver()
8985 pcbddc_coarse->benign_saddle_point = pcbddc->benign_have_null; in PCBDDCSetUpCoarseSolver()
8986 if (pcbddc_coarse->benign_saddle_point) { in PCBDDCSetUpCoarseSolver()
8995 st = st - n; in PCBDDCSetUpCoarseSolver()
9018 pcbddc_coarse->adaptive_userdefined = PETSC_TRUE; in PCBDDCSetUpCoarseSolver()
9019 if (pcbddc->adaptive_threshold[0] == 0.0) pcbddc_coarse->deluxe_zerorows = PETSC_TRUE; in PCBDDCSetUpCoarseSolver()
9025 PetscCall(MatIsSymmetricKnown(pc->pmat, &isset, &issym)); in PCBDDCSetUpCoarseSolver()
9027 PetscCall(MatIsHermitianKnown(pc->pmat, &isset, &isher)); in PCBDDCSetUpCoarseSolver()
9029 PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd)); in PCBDDCSetUpCoarseSolver()
9032 …if (pcbddc->benign_saddle_point && !pcbddc->benign_have_null) PetscCall(MatSetOption(coarse_mat, M… in PCBDDCSetUpCoarseSolver()
9034 PetscCall(MatViewFromOptions(coarse_mat, (PetscObject)pc, "-pc_bddc_coarse_mat_view")); in PCBDDCSetUpCoarseSolver()
9035 PetscCall(MatSetOptionsPrefix(coarse_mat, ((PetscObject)pcbddc->coarse_ksp)->prefix)); in PCBDDCSetUpCoarseSolver()
9036 PetscCall(KSPSetOperators(pcbddc->coarse_ksp, coarse_mat, coarse_mat)); in PCBDDCSetUpCoarseSolver()
9037 …if (pcbddc->dbg_flag) PetscCall(PetscViewerASCIISubtractTab(dbg_viewer, 2 * pcbddc->current_level)… in PCBDDCSetUpCoarseSolver()
9045 …intf(filename, PETSC_STATIC_ARRAY_LENGTH(filename), "coarse_mat_level%d.m",pcbddc->current_level)); in PCBDDCSetUpCoarseSolver()
9058 PetscInt i, d, N, n, cdim = pcbddc->mat_graph->cdim; in PCBDDCSetUpCoarseSolver()
9061 … PetscCheck(pcbddc->mat_graph->cloc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing local coordinates"); in PCBDDCSetUpCoarseSolver()
9062 PetscCall(VecGetSize(pcbddc->coarse_vec, &N)); in PCBDDCSetUpCoarseSolver()
9063 PetscCall(VecGetLocalSize(pcbddc->coarse_vec, &n)); in PCBDDCSetUpCoarseSolver()
9064 PetscCall(VecCreate(PetscObjectComm((PetscObject)pcbddc->coarse_vec), &gv)); in PCBDDCSetUpCoarseSolver()
9071 PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &is)); in PCBDDCSetUpCoarseSolver()
9076 … for (d = 0; d < cdim; d++) coords[cdim * i + d] = pcbddc->mat_graph->coords[cdim * idxs[i] + d]; in PCBDDCSetUpCoarseSolver()
9079 PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &is)); in PCBDDCSetUpCoarseSolver()
9089 if (pcbddc->coarse_ksp) { in PCBDDCSetUpCoarseSolver()
9093 PetscCall(KSPGetPC(pcbddc->coarse_ksp, &coarse_pc)); in PCBDDCSetUpCoarseSolver()
9116 if (pcbddc->coarse_ksp) { in PCBDDCSetUpCoarseSolver()
9119 PetscCall(KSPGetSolution(pcbddc->coarse_ksp, &csol)); in PCBDDCSetUpCoarseSolver()
9120 PetscCall(KSPGetRhs(pcbddc->coarse_ksp, &crhs)); in PCBDDCSetUpCoarseSolver()
9121 if (!csol) PetscCall(MatCreateVecs(coarse_mat, &pcbddc->coarse_ksp->vec_sol, NULL)); in PCBDDCSetUpCoarseSolver()
9122 if (!crhs) PetscCall(MatCreateVecs(coarse_mat, NULL, &pcbddc->coarse_ksp->vec_rhs)); in PCBDDCSetUpCoarseSolver()
9127 if (pcbddc->benign_null) { in PCBDDCSetUpCoarseSolver()
9128 PetscCall(VecSet(pcbddc->vec1_P, 0.)); in PCBDDCSetUpCoarseSolver()
9129 …for (i = 0; i < pcbddc->benign_n; i++) PetscCall(VecSetValue(pcbddc->vec1_P, pcbddc->local_primal_… in PCBDDCSetUpCoarseSolver()
9130 PetscCall(VecAssemblyBegin(pcbddc->vec1_P)); in PCBDDCSetUpCoarseSolver()
9131 PetscCall(VecAssemblyEnd(pcbddc->vec1_P)); in PCBDDCSetUpCoarseSolver()
9132 …PetscCall(VecScatterBegin(pcbddc->coarse_loc_to_glob, pcbddc->vec1_P, pcbddc->coarse_vec, INSERT_V… in PCBDDCSetUpCoarseSolver()
9133 …PetscCall(VecScatterEnd(pcbddc->coarse_loc_to_glob, pcbddc->vec1_P, pcbddc->coarse_vec, INSERT_VAL… in PCBDDCSetUpCoarseSolver()
9141 PetscCall(VecGetArrayRead(pcbddc->coarse_vec, (const PetscScalar **)&array)); in PCBDDCSetUpCoarseSolver()
9145 PetscCall(VecRestoreArrayRead(pcbddc->coarse_vec, (const PetscScalar **)&array)); in PCBDDCSetUpCoarseSolver()
9151 PetscCall(PetscLogEventEnd(PC_BDDC_CoarseSetUp[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpCoarseSolver()
9153 PetscCall(PetscLogEventBegin(PC_BDDC_CoarseSolver[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpCoarseSolver()
9154 if (pcbddc->coarse_ksp) { in PCBDDCSetUpCoarseSolver()
9165 PetscCall(KSPSetUp(pcbddc->coarse_ksp)); in PCBDDCSetUpCoarseSolver()
9167 PetscCall(PetscObjectTypeCompare((PetscObject)pcbddc->coarse_ksp, KSPPREONLY, &ispreonly)); in PCBDDCSetUpCoarseSolver()
9168 if (pcbddc->dbg_flag || (!ispreonly && pcbddc->use_coarse_estimates)) { in PCBDDCSetUpCoarseSolver()
9181 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pcbddc->coarse_ksp), &check_ksp)); in PCBDDCSetUpCoarseSolver()
9182 PetscCall(KSPSetNestLevel(check_ksp, pc->kspnestlevel)); in PCBDDCSetUpCoarseSolver()
9183 …PetscCall(PetscObjectIncrementTabLevel((PetscObject)check_ksp, (PetscObject)pcbddc->coarse_ksp, 0)… in PCBDDCSetUpCoarseSolver()
9184 PetscCall(KSPSetErrorIfNotConverged(pcbddc->coarse_ksp, PETSC_FALSE)); in PCBDDCSetUpCoarseSolver()
9186 PetscCall(KSPSetTolerances(check_ksp, 1.e-12, 1.e-12, PETSC_CURRENT, pcbddc->coarse_size)); in PCBDDCSetUpCoarseSolver()
9200 PetscCall(KSPGMRESSetRestart(check_ksp, pcbddc->coarse_size + 1)); in PCBDDCSetUpCoarseSolver()
9201 PetscCall(KSPGetOptionsPrefix(pcbddc->coarse_ksp, &prefix)); in PCBDDCSetUpCoarseSolver()
9206 PetscCall(KSPGetPC(pcbddc->coarse_ksp, &check_pc)); in PCBDDCSetUpCoarseSolver()
9217 PetscCall(PetscMalloc1(pcbddc->coarse_size + 1, &eigs_r)); in PCBDDCSetUpCoarseSolver()
9218 PetscCall(PetscMalloc1(pcbddc->coarse_size + 1, &eigs_c)); in PCBDDCSetUpCoarseSolver()
9219 … PetscCall(KSPComputeEigenvalues(check_ksp, pcbddc->coarse_size + 1, eigs_r, eigs_c, &neigs)); in PCBDDCSetUpCoarseSolver()
9221 lambda_max = eigs_r[neigs - 1]; in PCBDDCSetUpCoarseSolver()
9223 if (pcbddc->use_coarse_estimates) { in PCBDDCSetUpCoarseSolver()
9225 … PetscCall(KSPChebyshevSetEigenvalues(pcbddc->coarse_ksp, lambda_max + PETSC_SMALL, lambda_min)); in PCBDDCSetUpCoarseSolver()
9226 PetscCall(KSPRichardsonSetScale(pcbddc->coarse_ksp, 2.0 / (lambda_max + lambda_min))); in PCBDDCSetUpCoarseSolver()
9233 if (pcbddc->dbg_flag) { in PCBDDCSetUpCoarseSolver()
9234 … PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pcbddc->coarse_ksp)); in PCBDDCSetUpCoarseSolver()
9235 PetscCall(PetscViewerASCIIAddTab(dbg_viewer, 2 * (pcbddc->current_level + 1))); in PCBDDCSetUpCoarseSolver()
9236 PetscCall(VecAXPY(check_vec, -1.0, coarse_vec)); in PCBDDCSetUpCoarseSolver()
9240 …SCIIPrintf(dbg_viewer, "Coarse problem details (use estimates %d)\n", pcbddc->use_coarse_estimates… in PCBDDCSetUpCoarseSolver()
9241 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pcbddc->coarse_ksp, dbg_viewer)); in PCBDDCSetUpCoarseSolver()
9257 PetscCall(PetscViewerASCIISubtractTab(dbg_viewer, 2 * (pcbddc->current_level + 1))); in PCBDDCSetUpCoarseSolver()
9270 if (pcbddc->dbg_flag) { in PCBDDCSetUpCoarseSolver()
9273 …cCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "Coarse solver setup completed at level %" PetscI… in PCBDDCSetUpCoarseSolver()
9274 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCSetUpCoarseSolver()
9279 PetscCall(PetscLogEventEnd(PC_BDDC_CoarseSolver[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpCoarseSolver()
9285 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCComputePrimalNumbering()
9286 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCComputePrimalNumbering()
9294 …PetscCheck(!pcbddc->local_primal_size || pcbddc->local_primal_ref_node, PETSC_COMM_SELF, PETSC_ERR… in PCBDDCComputePrimalNumbering()
9295 …cCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc->pmat), pcbddc->local_primal_size_cc, pcbddc… in PCBDDCComputePrimalNumbering()
9296 PetscCall(ISLocalToGlobalMappingApplyIS(pcis->mapping, subset_n, &subset)); in PCBDDCComputePrimalNumbering()
9298 …cCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc->pmat), pcbddc->local_primal_size_cc, pcbddc… in PCBDDCComputePrimalNumbering()
9303 …->local_primal_size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of local primal indices comp… in PCBDDCComputePrimalNumbering()
9310 if (pcbddc->dbg_flag) { in PCBDDCComputePrimalNumbering()
9311 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCComputePrimalNumbering()
9312 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "--------------------------------------------… in PCBDDCComputePrimalNumbering()
9313 …PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "Size of coarse problem is %" PetscInt_FMT "\… in PCBDDCComputePrimalNumbering()
9314 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCComputePrimalNumbering()
9365 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCComputeFakeChange()
9366 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCComputeFakeChange()
9373 PetscCall(PCSetOperators(pcf, pc->mat, pc->pmat)); in PCBDDCComputeFakeChange()
9376 pcisf = (PC_IS *)pcf->data; in PCBDDCComputeFakeChange()
9377 pcbddcf = (PC_BDDC *)pcf->data; in PCBDDCComputeFakeChange()
9379 pcisf->is_B_local = pcis->is_B_local; in PCBDDCComputeFakeChange()
9380 pcisf->vec1_N = pcis->vec1_N; in PCBDDCComputeFakeChange()
9381 pcisf->BtoNmap = pcis->BtoNmap; in PCBDDCComputeFakeChange()
9382 pcisf->n = pcis->n; in PCBDDCComputeFakeChange()
9383 pcisf->n_B = pcis->n_B; in PCBDDCComputeFakeChange()
9385 PetscCall(PetscFree(pcbddcf->mat_graph)); in PCBDDCComputeFakeChange()
9386 PetscCall(PetscFree(pcbddcf->sub_schurs)); in PCBDDCComputeFakeChange()
9387 pcbddcf->mat_graph = graph ? graph : pcbddc->mat_graph; in PCBDDCComputeFakeChange()
9388 pcbddcf->sub_schurs = schurs; in PCBDDCComputeFakeChange()
9389 pcbddcf->adaptive_selection = schurs ? PETSC_TRUE : PETSC_FALSE; in PCBDDCComputeFakeChange()
9390 pcbddcf->adaptive_threshold[0] = pcbddc->adaptive_threshold[0]; in PCBDDCComputeFakeChange()
9391 pcbddcf->adaptive_threshold[1] = pcbddc->adaptive_threshold[1]; in PCBDDCComputeFakeChange()
9392 pcbddcf->adaptive_nmin = pcbddc->adaptive_nmin; in PCBDDCComputeFakeChange()
9393 pcbddcf->adaptive_nmax = pcbddc->adaptive_nmax; in PCBDDCComputeFakeChange()
9394 pcbddcf->use_faces = PETSC_TRUE; in PCBDDCComputeFakeChange()
9395 pcbddcf->use_change_of_basis = (PetscBool)!constraints; in PCBDDCComputeFakeChange()
9396 pcbddcf->use_change_on_faces = (PetscBool)!constraints; in PCBDDCComputeFakeChange()
9397 pcbddcf->use_qr_single = (PetscBool)!constraints; in PCBDDCComputeFakeChange()
9398 pcbddcf->fake_change = PETSC_TRUE; in PCBDDCComputeFakeChange()
9399 pcbddcf->dbg_flag = pcbddc->dbg_flag; in PCBDDCComputeFakeChange()
9404 *change = pcbddcf->ConstraintMatrix; in PCBDDCComputeFakeChange()
9405 …Call(ISCreateGeneral(PetscObjectComm((PetscObject)pc->pmat), pcbddcf->local_primal_size_cc, pcbddc… in PCBDDCComputeFakeChange()
9406 …Call(ISCreateGeneral(PetscObjectComm((PetscObject)pc->pmat), pcbddcf->local_primal_size_cc, pcbddc… in PCBDDCComputeFakeChange()
9407 if (change_with_qr) *change_with_qr = pcbddcf->use_qr_single; in PCBDDCComputeFakeChange()
9409 if (schurs) pcbddcf->sub_schurs = NULL; in PCBDDCComputeFakeChange()
9410 pcbddcf->ConstraintMatrix = NULL; in PCBDDCComputeFakeChange()
9411 pcbddcf->mat_graph = NULL; in PCBDDCComputeFakeChange()
9412 pcisf->is_B_local = NULL; in PCBDDCComputeFakeChange()
9413 pcisf->vec1_N = NULL; in PCBDDCComputeFakeChange()
9414 pcisf->BtoNmap = NULL; in PCBDDCComputeFakeChange()
9421 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCSetUpSubSchurs()
9422 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCSetUpSubSchurs()
9423 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; in PCBDDCSetUpSubSchurs()
9429 PetscCall(PetscLogEventBegin(PC_BDDC_Schurs[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpSubSchurs()
9432 if (pcbddc->sub_schurs_layers == -1) { in PCBDDCSetUpSubSchurs()
9436 if (pcbddc->sub_schurs_use_useradj && pcbddc->mat_graph->xadj) { in PCBDDCSetUpSubSchurs()
9437 used_xadj = pcbddc->mat_graph->xadj; in PCBDDCSetUpSubSchurs()
9438 used_adjncy = pcbddc->mat_graph->adjncy; in PCBDDCSetUpSubSchurs()
9439 } else if (pcbddc->computed_rowadj) { in PCBDDCSetUpSubSchurs()
9440 used_xadj = pcbddc->mat_graph->xadj; in PCBDDCSetUpSubSchurs()
9441 used_adjncy = pcbddc->mat_graph->adjncy; in PCBDDCSetUpSubSchurs()
9447 …PetscCall(MatGetRowIJ(pcbddc->local_mat, 0, PETSC_TRUE, PETSC_FALSE, &nvtxs, &xadj, &adjncy, &flg_… in PCBDDCSetUpSubSchurs()
9454 pcbddc->sub_schurs_layers = -1; in PCBDDCSetUpSubSchurs()
9458 …PetscCall(MatRestoreRowIJ(pcbddc->local_mat, 0, PETSC_TRUE, PETSC_FALSE, &nvtxs, &xadj, &adjncy, &… in PCBDDCSetUpSubSchurs()
9463 …PetscCall(MatCreateSchurComplement(pcis->A_II, pcis->pA_II, pcis->A_IB, pcis->A_BI, pcis->A_BB, &S… in PCBDDCSetUpSubSchurs()
9464 if (!sub_schurs->schur_explicit) { in PCBDDCSetUpSubSchurs()
9465 /* pcbddc->ksp_D up to date only if not using MatFactor with Schur complement support */ in PCBDDCSetUpSubSchurs()
9466 PetscCall(MatSchurComplementSetKSP(S_j, pcbddc->ksp_D)); in PCBDDCSetUpSubSchurs()
9467 …urs, NULL, S_j, PETSC_FALSE, used_xadj, used_adjncy, pcbddc->sub_schurs_layers, NULL, pcbddc->adap… in PCBDDCSetUpSubSchurs()
9473 PetscBool reuse_solvers = (PetscBool)!pcbddc->use_change_of_basis; in PCBDDCSetUpSubSchurs()
9477 if (!pcbddc->use_vertices && reuse_solvers) { in PCBDDCSetUpSubSchurs()
9480 PetscCall(ISGetLocalSize(sub_schurs->is_vertices, &n_vertices)); in PCBDDCSetUpSubSchurs()
9483 if (!pcbddc->benign_change_explicit) { in PCBDDCSetUpSubSchurs()
9484 benign_n = pcbddc->benign_n; in PCBDDCSetUpSubSchurs()
9488 …/* sub_schurs->change is a local object; instead, PCBDDCConstraintsSetUp and the quantities used i… in PCBDDCSetUpSubSchurs()
9490 …We assume that sub_schurs->change is created once, and then reused for different solves, unless th… in PCBDDCSetUpSubSchurs()
9491 if (pcbddc->adaptive_userdefined || (pcbddc->deluxe_zerorows && !pcbddc->use_change_of_basis)) { in PCBDDCSetUpSubSchurs()
9492 PetscBool have_loc_change = (PetscBool)(!!sub_schurs->change); in PCBDDCSetUpSubSchurs()
9498 …PetscCheck(!pcbddc->sub_schurs_rebuild, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot compute change of … in PCBDDCSetUpSubSchurs()
9499 …keChange(pc, PETSC_FALSE, NULL, NULL, &change, &change_primal, NULL, &sub_schurs->change_with_qr)); in PCBDDCSetUpSubSchurs()
9501 if (!pcbddc->use_deluxe_scaling) scaling = pcis->D; in PCBDDCSetUpSubSchurs()
9505 …PetscOptionsBegin(PetscObjectComm((PetscObject)iP), sub_schurs->prefix, "BDDC sub_schurs options",… in PCBDDCSetUpSubSchurs()
9506 …PetscCall(PetscOptionsBool("-sub_schurs_discrete_harmonic", NULL, NULL, discrete_harmonic, &discre… in PCBDDCSetUpSubSchurs()
9511 PetscCall(MatDuplicate(pcbddc->local_mat, MAT_COPY_VALUES, &A)); in PCBDDCSetUpSubSchurs()
9514 …->sub_schurs_exact_schur, used_xadj, used_adjncy, pcbddc->sub_schurs_layers, scaling, pcbddc->adap… in PCBDDCSetUpSubSchurs()
9515 pcbddc->benign_zerodiag_subs, change, change_primal)); in PCBDDCSetUpSubSchurs()
9518 …->local_mat, S_j, pcbddc->sub_schurs_exact_schur, used_xadj, used_adjncy, pcbddc->sub_schurs_layer… in PCBDDCSetUpSubSchurs()
9519 … pcbddc->benign_p0_lidx, pcbddc->benign_zerodiag_subs, change, change_primal)); in PCBDDCSetUpSubSchurs()
9528 PetscCall(PetscLogEventEnd(PC_BDDC_Schurs[pcbddc->current_level], pc, 0, 0, 0)); in PCBDDCSetUpSubSchurs()
9534 PC_IS *pcis = (PC_IS *)pc->data; in PCBDDCInitSubSchurs()
9535 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCInitSubSchurs()
9540 …if (pcbddc->sub_schurs_rebuild) { /* in case rebuild has been requested, it uses a graph generated… in PCBDDCInitSubSchurs()
9544 PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &verticesIS)); in PCBDDCInitSubSchurs()
9549 … PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &verticesIS)); in PCBDDCInitSubSchurs()
9551 …PetscCall(PCBDDCGraphInit(graph, pcbddc->mat_graph->l2gmap, pcbddc->mat_graph->nvtxs_global, pcbdd… in PCBDDCInitSubSchurs()
9552 …PetscCall(PCBDDCGraphSetUp(graph, pcbddc->mat_graph->custom_minimal_size, NULL, pcbddc->DirichletB… in PCBDDCInitSubSchurs()
9556 graph = pcbddc->mat_graph; in PCBDDCInitSubSchurs()
9559 if (pcbddc->dbg_flag && !pcbddc->sub_schurs_rebuild) { in PCBDDCInitSubSchurs()
9562 PetscCall(PCBDDCGraphASCIIView(graph, pcbddc->dbg_flag, pcbddc->dbg_viewer)); in PCBDDCInitSubSchurs()
9565 PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer)); in PCBDDCInitSubSchurs()
9566 …scViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "---------------------------------------------… in PCBDDCInitSubSchurs()
9567 …nizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d got %02" PetscInt_FMT " local candidate vertices (… in PCBDDCInitSubSchurs()
9568 …zedPrintf(pcbddc->dbg_viewer, "Subdomain %04d got %02" PetscInt_FMT " local candidate edges (%d… in PCBDDCInitSubSchurs()
9569 …zedPrintf(pcbddc->dbg_viewer, "Subdomain %04d got %02" PetscInt_FMT " local candidate faces (%d… in PCBDDCInitSubSchurs()
9570 PetscCall(PetscViewerFlush(pcbddc->dbg_viewer)); in PCBDDCInitSubSchurs()
9571 PetscCall(PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer)); in PCBDDCInitSubSchurs()
9576 if (!pcbddc->sub_schurs) PetscCall(PCBDDCSubSchursCreate(&pcbddc->sub_schurs)); in PCBDDCInitSubSchurs()
9577 …DDCSubSchursInit(pcbddc->sub_schurs, ((PetscObject)pc)->prefix, pcis->is_I_local, pcis->is_B_local… in PCBDDCInitSubSchurs()
9580 if (pcbddc->sub_schurs_rebuild) PetscCall(PCBDDCGraphDestroy(&graph)); in PCBDDCInitSubSchurs()
9586 Mat_IS *matis = (Mat_IS *)pc->pmat->data; in PCBDDCViewGlobalIS()
9587 PetscInt n = pc->pmat->rmap->n, ln, ni, st; in PCBDDCViewGlobalIS()
9593 PetscCall(MatGetOwnershipRange(pc->pmat, &st, NULL)); in PCBDDCViewGlobalIS()
9594 PetscCall(MatGetLocalSize(matis->A, NULL, &ln)); in PCBDDCViewGlobalIS()
9595 PetscCall(PetscArrayzero(matis->sf_leafdata, ln)); in PCBDDCViewGlobalIS()
9596 PetscCall(PetscArrayzero(matis->sf_rootdata, n)); in PCBDDCViewGlobalIS()
9601 matis->sf_leafdata[idxs[i]] = 1; in PCBDDCViewGlobalIS()
9604 …PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_SUM)… in PCBDDCViewGlobalIS()
9605 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_SUM)); in PCBDDCViewGlobalIS()
9608 if (matis->sf_rootdata[i]) matis->sf_rootdata[ln++] = i + st; in PCBDDCViewGlobalIS()
9610 …PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), ln, matis->sf_rootdata, PETSC_USE_POIN… in PCBDDCViewGlobalIS()
9619 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; in PCBDDCLoadOrViewCustomization()
9685 header[0] = (PetscInt)!!pcbddc->DirichletBoundariesLocal; in PCBDDCLoadOrViewCustomization()
9686 header[1] = (PetscInt)!!pcbddc->NeumannBoundariesLocal; in PCBDDCLoadOrViewCustomization()
9687 header[2] = pcbddc->n_ISForDofsLocal; in PCBDDCLoadOrViewCustomization()
9688 header[3] = (PetscInt)!!pcbddc->user_primal_vertices_local; in PCBDDCLoadOrViewCustomization()
9689 header[4] = (PetscInt)!!pcbddc->discretegradient; in PCBDDCLoadOrViewCustomization()
9690 header[5] = pcbddc->nedorder; in PCBDDCLoadOrViewCustomization()
9691 header[6] = pcbddc->nedfield; in PCBDDCLoadOrViewCustomization()
9692 header[7] = (PetscInt)pcbddc->nedglobal; in PCBDDCLoadOrViewCustomization()
9693 header[8] = (PetscInt)pcbddc->conforming; in PCBDDCLoadOrViewCustomization()
9694 header[9] = (PetscInt)!!pcbddc->divudotp; in PCBDDCLoadOrViewCustomization()
9695 header[10] = (PetscInt)pcbddc->divudotp_trans; in PCBDDCLoadOrViewCustomization()
9699 PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->DirichletBoundariesLocal, viewer)); in PCBDDCLoadOrViewCustomization()
9700 PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->NeumannBoundariesLocal, viewer)); in PCBDDCLoadOrViewCustomization()
9701 …for (PetscInt i = 0; i < header[2]; i++) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->ISForDofsLocal[i… in PCBDDCLoadOrViewCustomization()
9702 if (header[3]) PetscCall(PCBDDCViewGlobalIS(pc, pcbddc->user_primal_vertices_local, viewer)); in PCBDDCLoadOrViewCustomization()
9703 if (header[4]) PetscCall(MatView(pcbddc->discretegradient, viewer)); in PCBDDCLoadOrViewCustomization()
9704 if (header[9]) PetscCall(MatView(pcbddc->divudotp, viewer)); in PCBDDCLoadOrViewCustomization()
9722 PetscCall(PetscLayoutSetSize(rmap, A->rmap->N)); in MatMPIAIJRestrict()
9727 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), ren - rst, rst, 1, &rows)); in MatMPIAIJRestrict()
9738 PetscCall(MatSetSizes(*B, ren - rst, PETSC_DECIDE, PETSC_DECIDE, At->cmap->N)); in MatMPIAIJRestrict()
9740 PetscCall(PetscLayoutDestroy(&(*B)->rmap)); in MatMPIAIJRestrict()
9741 PetscCall(PetscLayoutSetUp((*B)->cmap)); in MatMPIAIJRestrict()
9742 a = (Mat_MPIAIJ *)At->data; in MatMPIAIJRestrict()
9743 b = (Mat_MPIAIJ *)(*B)->data; in MatMPIAIJRestrict()
9744 PetscCallMPI(MPI_Comm_size(ccomm, &b->size)); in MatMPIAIJRestrict()
9745 PetscCallMPI(MPI_Comm_rank(ccomm, &b->rank)); in MatMPIAIJRestrict()
9746 PetscCall(PetscObjectReference((PetscObject)a->A)); in MatMPIAIJRestrict()
9747 PetscCall(PetscObjectReference((PetscObject)a->B)); in MatMPIAIJRestrict()
9748 b->A = a->A; in MatMPIAIJRestrict()
9749 b->B = a->B; in MatMPIAIJRestrict()
9751 b->donotstash = a->donotstash; in MatMPIAIJRestrict()
9752 b->roworiented = a->roworiented; in MatMPIAIJRestrict()
9753 b->rowindices = NULL; in MatMPIAIJRestrict()
9754 b->rowvalues = NULL; in MatMPIAIJRestrict()
9755 b->getrowactive = PETSC_FALSE; in MatMPIAIJRestrict()
9757 (*B)->rmap = rmap; in MatMPIAIJRestrict()
9758 (*B)->factortype = A->factortype; in MatMPIAIJRestrict()
9759 (*B)->assembled = PETSC_TRUE; in MatMPIAIJRestrict()
9760 (*B)->insertmode = NOT_SET_VALUES; in MatMPIAIJRestrict()
9761 (*B)->preallocated = PETSC_TRUE; in MatMPIAIJRestrict()
9763 if (a->colmap) { in MatMPIAIJRestrict()
9765 PetscCall(PetscHMapIDuplicate(a->colmap, &b->colmap)); in MatMPIAIJRestrict()
9767 PetscCall(PetscMalloc1(At->cmap->N, &b->colmap)); in MatMPIAIJRestrict()
9768 PetscCall(PetscArraycpy(b->colmap, a->colmap, At->cmap->N)); in MatMPIAIJRestrict()
9770 } else b->colmap = NULL; in MatMPIAIJRestrict()
9771 if (a->garray) { in MatMPIAIJRestrict()
9773 len = a->B->cmap->n; in MatMPIAIJRestrict()
9774 PetscCall(PetscMalloc1(len + 1, &b->garray)); in MatMPIAIJRestrict()
9775 if (len) PetscCall(PetscArraycpy(b->garray, a->garray, len)); in MatMPIAIJRestrict()
9776 } else b->garray = NULL; in MatMPIAIJRestrict()
9778 PetscCall(PetscObjectReference((PetscObject)a->lvec)); in MatMPIAIJRestrict()
9779 b->lvec = a->lvec; in MatMPIAIJRestrict()
9782 PetscCall(VecGetLocalSize(b->lvec, &lsize)); in MatMPIAIJRestrict()
9783 PetscCall(ISCreateGeneral(ccomm, lsize, b->garray, PETSC_USE_POINTER, &from)); in MatMPIAIJRestrict()
9786 PetscCall(VecScatterCreate(gvec, from, b->lvec, to, &b->Mvctx)); in MatMPIAIJRestrict()
9812 PetscInt ni, *di, *dj, m = A->rmap->n, c, *ldata, *rdata; in MatAIJExtractRows()
9822 PetscCall(PetscSFSetGraphLayout(sf, A->rmap, ni, NULL, PETSC_USE_POINTER, idxs)); in MatAIJExtractRows()
9826 da = (Mat_SeqAIJ *)A_loc->data; in MatAIJExtractRows()
9829 rdata[2 * i + 0] = da->i[i + 1] - da->i[i]; in MatAIJExtractRows()
9830 rdata[2 * i + 1] = da->i[i]; in MatAIJExtractRows()
9855 …PetscCall(PetscSFSetGraph(sf, da->i[m], di[ni], NULL, PETSC_USE_POINTER, remotes, PETSC_USE_POINTE… in MatAIJExtractRows()
9856 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, da->j, dj, MPI_REPLACE)); in MatAIJExtractRows()
9857 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, da->j, dj, MPI_REPLACE)); in MatAIJExtractRows()
9858 PetscCall(PetscSFBcastBegin(sf, MPIU_SCALAR, da->a, daa, MPI_REPLACE)); in MatAIJExtractRows()
9859 PetscCall(PetscSFBcastEnd(sf, MPIU_SCALAR, da->a, daa, MPI_REPLACE)); in MatAIJExtractRows()
9861 …PetscCall(MatCreateMPIAIJWithArrays(comm, ni, A->cmap->n, PETSC_DECIDE, A->cmap->N, di, dj, daa, s… in MatAIJExtractRows()