Actual source code: bddcfetidp.c

      1: #include <petsc/private/pcbddcimpl.h>
      2: #include <petsc/private/pcbddcprivateimpl.h>
      3: #include <petscblaslapack.h>

      5: static PetscErrorCode MatMult_BDdelta_deluxe_nonred(Mat A, Vec x, Vec y)
      6: {
      7:   BDdelta_DN ctx;

      9:   PetscFunctionBegin;
     10:   PetscCall(MatShellGetContext(A, &ctx));
     11:   PetscCall(MatMultTranspose(ctx->BD, x, ctx->work));
     12:   PetscCall(KSPSolveTranspose(ctx->kBD, ctx->work, y));
     13:   /* No PC so cannot propagate up failure in KSPSolveTranspose() */
     14:   PetscFunctionReturn(PETSC_SUCCESS);
     15: }

     17: static PetscErrorCode MatMultTranspose_BDdelta_deluxe_nonred(Mat A, Vec x, Vec y)
     18: {
     19:   BDdelta_DN ctx;

     21:   PetscFunctionBegin;
     22:   PetscCall(MatShellGetContext(A, &ctx));
     23:   PetscCall(KSPSolve(ctx->kBD, x, ctx->work));
     24:   /* No PC so cannot propagate up failure in KSPSolve() */
     25:   PetscCall(MatMult(ctx->BD, ctx->work, y));
     26:   PetscFunctionReturn(PETSC_SUCCESS);
     27: }

     29: static PetscErrorCode MatDestroy_BDdelta_deluxe_nonred(Mat A)
     30: {
     31:   BDdelta_DN ctx;

     33:   PetscFunctionBegin;
     34:   PetscCall(MatShellGetContext(A, &ctx));
     35:   PetscCall(MatDestroy(&ctx->BD));
     36:   PetscCall(KSPDestroy(&ctx->kBD));
     37:   PetscCall(VecDestroy(&ctx->work));
     38:   PetscCall(PetscFree(ctx));
     39:   PetscFunctionReturn(PETSC_SUCCESS);
     40: }

     42: PetscErrorCode PCBDDCCreateFETIDPMatContext(PC pc, FETIDPMat_ctx *fetidpmat_ctx)
     43: {
     44:   FETIDPMat_ctx newctx;

     46:   PetscFunctionBegin;
     47:   PetscCall(PetscNew(&newctx));
     48:   /* increase the reference count for BDDC preconditioner */
     49:   PetscCall(PetscObjectReference((PetscObject)pc));
     50:   newctx->pc     = pc;
     51:   *fetidpmat_ctx = newctx;
     52:   PetscFunctionReturn(PETSC_SUCCESS);
     53: }

     55: PetscErrorCode PCBDDCCreateFETIDPPCContext(PC pc, FETIDPPC_ctx *fetidppc_ctx)
     56: {
     57:   FETIDPPC_ctx newctx;

     59:   PetscFunctionBegin;
     60:   PetscCall(PetscNew(&newctx));
     61:   /* increase the reference count for BDDC preconditioner */
     62:   PetscCall(PetscObjectReference((PetscObject)pc));
     63:   newctx->pc    = pc;
     64:   *fetidppc_ctx = newctx;
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: PetscErrorCode PCBDDCDestroyFETIDPMat(Mat A)
     69: {
     70:   FETIDPMat_ctx mat_ctx;

     72:   PetscFunctionBegin;
     73:   PetscCall(MatShellGetContext(A, &mat_ctx));
     74:   PetscCall(VecDestroy(&mat_ctx->lambda_local));
     75:   PetscCall(VecDestroy(&mat_ctx->temp_solution_D));
     76:   PetscCall(VecDestroy(&mat_ctx->temp_solution_B));
     77:   PetscCall(MatDestroy(&mat_ctx->B_delta));
     78:   PetscCall(MatDestroy(&mat_ctx->B_Ddelta));
     79:   PetscCall(ISDestroy(&mat_ctx->lP_I));
     80:   PetscCall(ISDestroy(&mat_ctx->lP_B));
     81:   PetscCall(MatDestroy(&mat_ctx->B_BB));
     82:   PetscCall(MatDestroy(&mat_ctx->B_BI));
     83:   PetscCall(MatDestroy(&mat_ctx->Bt_BB));
     84:   PetscCall(MatDestroy(&mat_ctx->Bt_BI));
     85:   PetscCall(MatDestroy(&mat_ctx->C));
     86:   PetscCall(VecDestroy(&mat_ctx->rhs_flip));
     87:   PetscCall(VecDestroy(&mat_ctx->vP));
     88:   PetscCall(VecDestroy(&mat_ctx->xPg));
     89:   PetscCall(VecDestroy(&mat_ctx->yPg));
     90:   PetscCall(VecScatterDestroy(&mat_ctx->l2g_lambda));
     91:   PetscCall(VecScatterDestroy(&mat_ctx->l2g_lambda_only));
     92:   PetscCall(VecScatterDestroy(&mat_ctx->l2g_p));
     93:   PetscCall(VecScatterDestroy(&mat_ctx->g2g_p));
     94:   PetscCall(PCDestroy(&mat_ctx->pc)); /* decrease PCBDDC reference count */
     95:   PetscCall(ISDestroy(&mat_ctx->pressure));
     96:   PetscCall(ISDestroy(&mat_ctx->lagrange));
     97:   PetscCall(PetscFree(mat_ctx));
     98:   PetscFunctionReturn(PETSC_SUCCESS);
     99: }

    101: PetscErrorCode PCBDDCDestroyFETIDPPC(PC pc)
    102: {
    103:   FETIDPPC_ctx pc_ctx;

    105:   PetscFunctionBegin;
    106:   PetscCall(PCShellGetContext(pc, &pc_ctx));
    107:   PetscCall(VecDestroy(&pc_ctx->lambda_local));
    108:   PetscCall(MatDestroy(&pc_ctx->B_Ddelta));
    109:   PetscCall(VecScatterDestroy(&pc_ctx->l2g_lambda));
    110:   PetscCall(MatDestroy(&pc_ctx->S_j));
    111:   PetscCall(PCDestroy(&pc_ctx->pc)); /* decrease PCBDDC reference count */
    112:   PetscCall(VecDestroy(&pc_ctx->xPg));
    113:   PetscCall(VecDestroy(&pc_ctx->yPg));
    114:   PetscCall(PetscFree(pc_ctx));
    115:   PetscFunctionReturn(PETSC_SUCCESS);
    116: }

    118: PetscErrorCode PCBDDCSetupFETIDPMatContext(FETIDPMat_ctx fetidpmat_ctx)
    119: {
    120:   PC_IS          *pcis      = (PC_IS *)fetidpmat_ctx->pc->data;
    121:   PC_BDDC        *pcbddc    = (PC_BDDC *)fetidpmat_ctx->pc->data;
    122:   PCBDDCGraph     mat_graph = pcbddc->mat_graph;
    123:   Mat_IS         *matis     = (Mat_IS *)fetidpmat_ctx->pc->pmat->data;
    124:   MPI_Comm        comm;
    125:   Mat             ScalingMat, BD1, BD2;
    126:   Vec             fetidp_global;
    127:   IS              IS_l2g_lambda;
    128:   IS              subset, subset_mult, subset_n, isvert;
    129:   PetscBool       skip_node, fully_redundant;
    130:   PetscInt        i, j, k, s, n_boundary_dofs, n_global_lambda, n_vertices, partial_sum;
    131:   PetscInt        cum, n_local_lambda, n_lambda_for_dof, dual_size, n_neg_values, n_pos_values, buf_size;
    132:   PetscMPIInt     rank, size, neigh;
    133:   PetscScalar     scalar_value;
    134:   const PetscInt *vertex_indices;
    135:   PetscInt       *dual_dofs_boundary_indices, *aux_local_numbering_1;
    136:   const PetscInt *aux_global_numbering;
    137:   PetscInt       *aux_sums, *cols_B_delta, *l2g_indices;
    138:   PetscScalar    *array, *scaling_factors, *vals_B_delta;
    139:   PetscScalar   **all_factors;
    140:   PetscInt       *aux_local_numbering_2;
    141:   PetscInt       *count, **neighbours_set;
    142:   PetscLayout     llay;

    144:   /* saddlepoint */
    145:   ISLocalToGlobalMapping l2gmap_p;
    146:   PetscLayout            play;
    147:   IS                     gP, pP;
    148:   PetscInt               nPl, nPg, nPgl;

    150:   PetscFunctionBegin;
    151:   PetscCall(PetscObjectGetComm((PetscObject)fetidpmat_ctx->pc, &comm));
    152:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    153:   PetscCallMPI(MPI_Comm_size(comm, &size));

    155:   /* saddlepoint */
    156:   nPl      = 0;
    157:   nPg      = 0;
    158:   nPgl     = 0;
    159:   gP       = NULL;
    160:   pP       = NULL;
    161:   l2gmap_p = NULL;
    162:   play     = NULL;
    163:   PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_pP", (PetscObject *)&pP));
    164:   if (pP) { /* saddle point */
    165:     /* subdomain pressures in global numbering */
    166:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_gP", (PetscObject *)&gP));
    167:     PetscCheck(gP, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gP not present");
    168:     PetscCall(ISGetLocalSize(gP, &nPl));
    169:     PetscCall(VecCreate(PETSC_COMM_SELF, &fetidpmat_ctx->vP));
    170:     PetscCall(VecSetSizes(fetidpmat_ctx->vP, nPl, nPl));
    171:     PetscCall(VecSetType(fetidpmat_ctx->vP, VECSTANDARD));
    172:     PetscCall(VecSetUp(fetidpmat_ctx->vP));

    174:     /* pressure matrix */
    175:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_C", (PetscObject *)&fetidpmat_ctx->C));
    176:     if (!fetidpmat_ctx->C) { /* null pressure block, compute layout and global numbering for pressures */
    177:       IS Pg;

    179:       PetscCall(ISRenumber(gP, NULL, &nPg, &Pg));
    180:       PetscCall(ISLocalToGlobalMappingCreateIS(Pg, &l2gmap_p));
    181:       PetscCall(ISDestroy(&Pg));
    182:       PetscCall(PetscLayoutCreate(comm, &play));
    183:       PetscCall(PetscLayoutSetBlockSize(play, 1));
    184:       PetscCall(PetscLayoutSetSize(play, nPg));
    185:       PetscCall(ISGetLocalSize(pP, &nPgl));
    186:       PetscCall(PetscLayoutSetLocalSize(play, nPgl));
    187:       PetscCall(PetscLayoutSetUp(play));
    188:     } else {
    189:       PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->C));
    190:       PetscCall(MatISGetLocalToGlobalMapping(fetidpmat_ctx->C, &l2gmap_p, NULL));
    191:       PetscCall(PetscObjectReference((PetscObject)l2gmap_p));
    192:       PetscCall(MatGetSize(fetidpmat_ctx->C, &nPg, NULL));
    193:       PetscCall(MatGetLocalSize(fetidpmat_ctx->C, NULL, &nPgl));
    194:       PetscCall(MatGetLayouts(fetidpmat_ctx->C, NULL, &llay));
    195:       PetscCall(PetscLayoutReference(llay, &play));
    196:     }
    197:     PetscCall(VecCreateMPIWithArray(comm, 1, nPgl, nPg, NULL, &fetidpmat_ctx->xPg));
    198:     PetscCall(VecCreateMPIWithArray(comm, 1, nPgl, nPg, NULL, &fetidpmat_ctx->yPg));

    200:     /* import matrices for pressures coupling */
    201:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_B_BI", (PetscObject *)&fetidpmat_ctx->B_BI));
    202:     PetscCheck(fetidpmat_ctx->B_BI, PETSC_COMM_SELF, PETSC_ERR_PLIB, "B_BI not present");
    203:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->B_BI));

    205:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_B_BB", (PetscObject *)&fetidpmat_ctx->B_BB));
    206:     PetscCheck(fetidpmat_ctx->B_BB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "B_BB not present");
    207:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->B_BB));

    209:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_Bt_BI", (PetscObject *)&fetidpmat_ctx->Bt_BI));
    210:     PetscCheck(fetidpmat_ctx->Bt_BI, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bt_BI not present");
    211:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BI));

    213:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_Bt_BB", (PetscObject *)&fetidpmat_ctx->Bt_BB));
    214:     PetscCheck(fetidpmat_ctx->Bt_BB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bt_BB not present");
    215:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BB));

    217:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_flip", (PetscObject *)&fetidpmat_ctx->rhs_flip));
    218:     if (fetidpmat_ctx->rhs_flip) PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->rhs_flip));

    220:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_lP_I", (PetscObject *)&fetidpmat_ctx->lP_I));
    221:     PetscCheck(fetidpmat_ctx->lP_I, PETSC_COMM_SELF, PETSC_ERR_PLIB, "lP_I not present");
    222:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->lP_I));

    224:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_lP_B", (PetscObject *)&fetidpmat_ctx->lP_B));
    225:     PetscCheck(fetidpmat_ctx->lP_B, PETSC_COMM_SELF, PETSC_ERR_PLIB, "lP_B not present");
    226:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->lP_B));
    227:   }

    229:   /* Default type of lagrange multipliers is non-redundant */
    230:   fully_redundant = fetidpmat_ctx->fully_redundant;

    232:   /* Evaluate local and global number of lagrange multipliers */
    233:   PetscCall(VecSet(pcis->vec1_N, 0.0));
    234:   n_local_lambda  = 0;
    235:   partial_sum     = 0;
    236:   n_boundary_dofs = 0;

    238:   /* Get Vertices used to define the BDDC */
    239:   PetscCall(PCBDDCGraphGetCandidatesIS(mat_graph, NULL, NULL, NULL, NULL, &isvert));
    240:   PetscCall(ISGetLocalSize(isvert, &n_vertices));
    241:   PetscCall(ISGetIndices(isvert, &vertex_indices));

    243:   dual_size = pcis->n_B - n_vertices;
    244:   PetscCall(PetscMalloc1(dual_size, &dual_dofs_boundary_indices));
    245:   PetscCall(PetscMalloc1(dual_size, &aux_local_numbering_1));
    246:   PetscCall(PetscMalloc1(dual_size, &aux_local_numbering_2));

    248:   /* the code below does not support multiple subdomains per process
    249:      error out in this case
    250:      TODO: I guess I can use PetscSFGetMultiSF and the code will be easier and more general */
    251:   PetscCall(PetscMalloc2(pcis->n, &count, pcis->n, &neighbours_set));
    252:   for (i = 0, j = 0; i < pcis->n; i++) j += mat_graph->nodes[i].count;
    253:   if (pcis->n) PetscCall(PetscMalloc1(j, &neighbours_set[0]));
    254:   for (i = 0; i < pcis->n; i++) {
    255:     PCBDDCGraphNode *node = &mat_graph->nodes[i];

    257:     count[i] = 0;
    258:     for (j = 0; j < node->count; j++) {
    259:       if (node->neighbours_set[j] == rank) continue;
    260:       neighbours_set[i][count[i]++] = node->neighbours_set[j];
    261:     }
    262:     PetscCheck(count[i] == node->count - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Multiple subdomains per process not supported");
    263:     s = count[i];
    264:     PetscCall(PetscSortRemoveDupsInt(count + i, neighbours_set[i]));
    265:     PetscCheck(s == count[i], PETSC_COMM_SELF, PETSC_ERR_SUP, "Multiple subdomains per process not supported");
    266:     if (i != pcis->n - 1) neighbours_set[i + 1] = neighbours_set[i] + count[i];
    267:   }

    269:   PetscCall(VecGetArray(pcis->vec1_N, &array));
    270:   for (i = 0, s = 0; i < pcis->n; i++) {
    271:     j = count[i]; /* RECALL: count[i] does not count myself */
    272:     if (j > 0) n_boundary_dofs++;
    273:     skip_node = PETSC_FALSE;
    274:     if (s < n_vertices && vertex_indices[s] == i) { /* it works for a sorted set of vertices */
    275:       skip_node = PETSC_TRUE;
    276:       s++;
    277:     }
    278:     if (j < 1) skip_node = PETSC_TRUE;
    279:     if (mat_graph->nodes[i].special_dof == PCBDDCGRAPH_DIRICHLET_MARK) skip_node = PETSC_TRUE;
    280:     if (!skip_node) {
    281:       if (fully_redundant) {
    282:         /* fully redundant set of lagrange multipliers */
    283:         n_lambda_for_dof = (j * (j + 1)) / 2;
    284:       } else {
    285:         n_lambda_for_dof = j;
    286:       }
    287:       n_local_lambda += j;
    288:       /* needed to evaluate global number of lagrange multipliers */
    289:       array[i] = (1.0 * n_lambda_for_dof) / (j + 1.0); /* already scaled for the next global sum */
    290:       /* store some data needed */
    291:       dual_dofs_boundary_indices[partial_sum] = n_boundary_dofs - 1;
    292:       aux_local_numbering_1[partial_sum]      = i;
    293:       aux_local_numbering_2[partial_sum]      = n_lambda_for_dof;
    294:       partial_sum++;
    295:     }
    296:   }
    297:   PetscCall(VecRestoreArray(pcis->vec1_N, &array));
    298:   PetscCall(ISRestoreIndices(isvert, &vertex_indices));
    299:   PetscCall(PCBDDCGraphRestoreCandidatesIS(mat_graph, NULL, NULL, NULL, NULL, &isvert));
    300:   dual_size = partial_sum;

    302:   /* compute global ordering of lagrange multipliers and associate l2g map */
    303:   PetscCall(ISCreateGeneral(comm, partial_sum, aux_local_numbering_1, PETSC_COPY_VALUES, &subset_n));
    304:   PetscCall(ISLocalToGlobalMappingApplyIS(pcis->mapping, subset_n, &subset));
    305:   PetscCall(ISDestroy(&subset_n));
    306:   PetscCall(ISCreateGeneral(comm, partial_sum, aux_local_numbering_2, PETSC_OWN_POINTER, &subset_mult));
    307:   PetscCall(ISRenumber(subset, subset_mult, &fetidpmat_ctx->n_lambda, &subset_n));
    308:   PetscCall(ISDestroy(&subset));

    310:   if (PetscDefined(USE_DEBUG)) {
    311:     PetscCall(VecSet(pcis->vec1_global, 0.0));
    312:     PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    313:     PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    314:     PetscCall(VecSum(pcis->vec1_global, &scalar_value));
    315:     i = (PetscInt)PetscRealPart(scalar_value);
    316:     PetscCheck(i == fetidpmat_ctx->n_lambda, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Global number of multipliers mismatch! (%" PetscInt_FMT " != %" PetscInt_FMT ")", fetidpmat_ctx->n_lambda, i);
    317:   }

    319:   /* init data for scaling factors exchange */
    320:   if (!pcbddc->use_deluxe_scaling) {
    321:     PetscInt    *ptrs_buffer, neigh_position;
    322:     PetscScalar *send_buffer, *recv_buffer;
    323:     MPI_Request *send_reqs, *recv_reqs;
    324:     PetscMPIInt  nreqs;

    326:     partial_sum = 0;
    327:     PetscCall(PetscMalloc1(pcis->n_neigh, &ptrs_buffer));
    328:     PetscCall(PetscMalloc1(PetscMax(pcis->n_neigh - 1, 0), &send_reqs));
    329:     PetscCall(PetscMalloc1(PetscMax(pcis->n_neigh - 1, 0), &recv_reqs));
    330:     PetscCall(PetscMalloc1(pcis->n + 1, &all_factors));
    331:     if (pcis->n_neigh > 0) ptrs_buffer[0] = 0;
    332:     for (i = 1; i < pcis->n_neigh; i++) {
    333:       partial_sum += pcis->n_shared[i];
    334:       ptrs_buffer[i] = ptrs_buffer[i - 1] + pcis->n_shared[i];
    335:     }
    336:     PetscCall(PetscMalloc1(partial_sum, &send_buffer));
    337:     PetscCall(PetscMalloc1(partial_sum, &recv_buffer));
    338:     PetscCall(PetscMalloc1(partial_sum, &all_factors[0]));
    339:     for (i = 0; i < pcis->n - 1; i++) {
    340:       j                  = count[i];
    341:       all_factors[i + 1] = all_factors[i] + j;
    342:     }

    344:     /* scatter B scaling to N vec */
    345:     PetscCall(VecScatterBegin(pcis->N_to_B, pcis->D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
    346:     PetscCall(VecScatterEnd(pcis->N_to_B, pcis->D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
    347:     /* communications */
    348:     PetscCall(VecGetArrayRead(pcis->vec1_N, (const PetscScalar **)&array));
    349:     for (i = 1; i < pcis->n_neigh; i++) {
    350:       for (j = 0; j < pcis->n_shared[i]; j++) send_buffer[ptrs_buffer[i - 1] + j] = array[pcis->shared[i][j]];
    351:       buf_size = ptrs_buffer[i] - ptrs_buffer[i - 1];
    352:       PetscCall(PetscMPIIntCast(pcis->neigh[i], &neigh));
    353:       PetscCallMPI(MPIU_Isend(&send_buffer[ptrs_buffer[i - 1]], buf_size, MPIU_SCALAR, neigh, 0, comm, &send_reqs[i - 1]));
    354:       PetscCallMPI(MPIU_Irecv(&recv_buffer[ptrs_buffer[i - 1]], buf_size, MPIU_SCALAR, neigh, 0, comm, &recv_reqs[i - 1]));
    355:     }
    356:     PetscCall(VecRestoreArrayRead(pcis->vec1_N, (const PetscScalar **)&array));
    357:     PetscCall(PetscMPIIntCast(pcis->n_neigh - 1, &nreqs));
    358:     if (pcis->n_neigh > 0) PetscCallMPI(MPI_Waitall(nreqs, recv_reqs, MPI_STATUSES_IGNORE));
    359:     /* put values in correct places */
    360:     for (i = 1; i < pcis->n_neigh; i++) {
    361:       for (j = 0; j < pcis->n_shared[i]; j++) {
    362:         k              = pcis->shared[i][j];
    363:         neigh_position = 0;
    364:         while (neighbours_set[k][neigh_position] != pcis->neigh[i]) neigh_position++;
    365:         all_factors[k][neigh_position] = recv_buffer[ptrs_buffer[i - 1] + j];
    366:       }
    367:     }
    368:     if (pcis->n_neigh > 0) PetscCallMPI(MPI_Waitall(nreqs, send_reqs, MPI_STATUSES_IGNORE));
    369:     PetscCall(PetscFree(send_reqs));
    370:     PetscCall(PetscFree(recv_reqs));
    371:     PetscCall(PetscFree(send_buffer));
    372:     PetscCall(PetscFree(recv_buffer));
    373:     PetscCall(PetscFree(ptrs_buffer));
    374:   }

    376:   /* Compute B and B_delta (local actions) */
    377:   PetscCall(PetscMalloc1(pcis->n_neigh, &aux_sums));
    378:   PetscCall(PetscMalloc1(n_local_lambda, &l2g_indices));
    379:   PetscCall(PetscMalloc1(n_local_lambda, &vals_B_delta));
    380:   PetscCall(PetscMalloc1(n_local_lambda, &cols_B_delta));
    381:   if (!pcbddc->use_deluxe_scaling) {
    382:     PetscCall(PetscMalloc1(n_local_lambda, &scaling_factors));
    383:   } else {
    384:     scaling_factors = NULL;
    385:     all_factors     = NULL;
    386:   }
    387:   PetscCall(ISGetIndices(subset_n, &aux_global_numbering));
    388:   partial_sum = 0;
    389:   cum         = 0;
    390:   for (i = 0; i < dual_size; i++) {
    391:     n_global_lambda = aux_global_numbering[cum];
    392:     j               = count[aux_local_numbering_1[i]];
    393:     aux_sums[0]     = 0;
    394:     for (s = 1; s < j; s++) aux_sums[s] = aux_sums[s - 1] + j - s + 1;
    395:     if (all_factors) array = all_factors[aux_local_numbering_1[i]];
    396:     n_neg_values = 0;
    397:     while (n_neg_values < j && neighbours_set[aux_local_numbering_1[i]][n_neg_values] < rank) n_neg_values++;
    398:     n_pos_values = j - n_neg_values;
    399:     if (fully_redundant) {
    400:       for (s = 0; s < n_neg_values; s++) {
    401:         l2g_indices[partial_sum + s]  = aux_sums[s] + n_neg_values - s - 1 + n_global_lambda;
    402:         cols_B_delta[partial_sum + s] = dual_dofs_boundary_indices[i];
    403:         vals_B_delta[partial_sum + s] = -1.0;
    404:         if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum + s] = array[s];
    405:       }
    406:       for (s = 0; s < n_pos_values; s++) {
    407:         l2g_indices[partial_sum + s + n_neg_values]  = aux_sums[n_neg_values] + s + n_global_lambda;
    408:         cols_B_delta[partial_sum + s + n_neg_values] = dual_dofs_boundary_indices[i];
    409:         vals_B_delta[partial_sum + s + n_neg_values] = 1.0;
    410:         if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum + s + n_neg_values] = array[s + n_neg_values];
    411:       }
    412:       partial_sum += j;
    413:     } else {
    414:       /* l2g_indices and default cols and vals of B_delta */
    415:       for (s = 0; s < j; s++) {
    416:         l2g_indices[partial_sum + s]  = n_global_lambda + s;
    417:         cols_B_delta[partial_sum + s] = dual_dofs_boundary_indices[i];
    418:         vals_B_delta[partial_sum + s] = 0.0;
    419:       }
    420:       /* B_delta */
    421:       if (n_neg_values > 0) { /* there's a rank next to me to the left */
    422:         vals_B_delta[partial_sum + n_neg_values - 1] = -1.0;
    423:       }
    424:       if (n_neg_values < j) { /* there's a rank next to me to the right */
    425:         vals_B_delta[partial_sum + n_neg_values] = 1.0;
    426:       }
    427:       /* scaling as in Klawonn-Widlund 1999 */
    428:       if (!pcbddc->use_deluxe_scaling) {
    429:         for (s = 0; s < n_neg_values; s++) {
    430:           scalar_value = 0.0;
    431:           for (k = 0; k < s + 1; k++) scalar_value += array[k];
    432:           scaling_factors[partial_sum + s] = -scalar_value;
    433:         }
    434:         for (s = 0; s < n_pos_values; s++) {
    435:           scalar_value = 0.0;
    436:           for (k = s + n_neg_values; k < j; k++) scalar_value += array[k];
    437:           scaling_factors[partial_sum + s + n_neg_values] = scalar_value;
    438:         }
    439:       }
    440:       partial_sum += j;
    441:     }
    442:     cum += aux_local_numbering_2[i];
    443:   }
    444:   PetscCall(ISRestoreIndices(subset_n, &aux_global_numbering));
    445:   PetscCall(ISDestroy(&subset_mult));
    446:   PetscCall(ISDestroy(&subset_n));
    447:   PetscCall(PetscFree(aux_sums));
    448:   PetscCall(PetscFree(aux_local_numbering_1));
    449:   PetscCall(PetscFree(dual_dofs_boundary_indices));
    450:   if (all_factors) {
    451:     PetscCall(PetscFree(all_factors[0]));
    452:     PetscCall(PetscFree(all_factors));
    453:   }
    454:   if (pcis->n) PetscCall(PetscFree(neighbours_set[0]));
    455:   PetscCall(PetscFree2(count, neighbours_set));

    457:   /* Create local part of B_delta */
    458:   PetscCall(MatCreate(PETSC_COMM_SELF, &fetidpmat_ctx->B_delta));
    459:   PetscCall(MatSetSizes(fetidpmat_ctx->B_delta, n_local_lambda, pcis->n_B, n_local_lambda, pcis->n_B));
    460:   PetscCall(MatSetType(fetidpmat_ctx->B_delta, MATSEQAIJ));
    461:   PetscCall(MatSeqAIJSetPreallocation(fetidpmat_ctx->B_delta, 1, NULL));
    462:   PetscCall(MatSetOption(fetidpmat_ctx->B_delta, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
    463:   for (i = 0; i < n_local_lambda; i++) PetscCall(MatSetValue(fetidpmat_ctx->B_delta, i, cols_B_delta[i], vals_B_delta[i], INSERT_VALUES));
    464:   PetscCall(PetscFree(vals_B_delta));
    465:   PetscCall(MatAssemblyBegin(fetidpmat_ctx->B_delta, MAT_FINAL_ASSEMBLY));
    466:   PetscCall(MatAssemblyEnd(fetidpmat_ctx->B_delta, MAT_FINAL_ASSEMBLY));

    468:   BD1 = NULL;
    469:   BD2 = NULL;
    470:   if (fully_redundant) {
    471:     PetscCheck(!pcbddc->use_deluxe_scaling, comm, PETSC_ERR_SUP, "Deluxe FETIDP with fully-redundant multipliers to be implemented");
    472:     PetscCall(MatCreate(PETSC_COMM_SELF, &ScalingMat));
    473:     PetscCall(MatSetSizes(ScalingMat, n_local_lambda, n_local_lambda, n_local_lambda, n_local_lambda));
    474:     PetscCall(MatSetType(ScalingMat, MATSEQAIJ));
    475:     PetscCall(MatSeqAIJSetPreallocation(ScalingMat, 1, NULL));
    476:     for (i = 0; i < n_local_lambda; i++) PetscCall(MatSetValue(ScalingMat, i, i, scaling_factors[i], INSERT_VALUES));
    477:     PetscCall(MatAssemblyBegin(ScalingMat, MAT_FINAL_ASSEMBLY));
    478:     PetscCall(MatAssemblyEnd(ScalingMat, MAT_FINAL_ASSEMBLY));
    479:     PetscCall(MatMatMult(ScalingMat, fetidpmat_ctx->B_delta, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &fetidpmat_ctx->B_Ddelta));
    480:     PetscCall(MatDestroy(&ScalingMat));
    481:   } else {
    482:     PetscCall(MatCreate(PETSC_COMM_SELF, &fetidpmat_ctx->B_Ddelta));
    483:     PetscCall(MatSetSizes(fetidpmat_ctx->B_Ddelta, n_local_lambda, pcis->n_B, n_local_lambda, pcis->n_B));
    484:     if (!pcbddc->use_deluxe_scaling || !pcbddc->sub_schurs) {
    485:       PetscCall(MatSetType(fetidpmat_ctx->B_Ddelta, MATSEQAIJ));
    486:       PetscCall(MatSeqAIJSetPreallocation(fetidpmat_ctx->B_Ddelta, 1, NULL));
    487:       for (i = 0; i < n_local_lambda; i++) PetscCall(MatSetValue(fetidpmat_ctx->B_Ddelta, i, cols_B_delta[i], scaling_factors[i], INSERT_VALUES));
    488:       PetscCall(MatAssemblyBegin(fetidpmat_ctx->B_Ddelta, MAT_FINAL_ASSEMBLY));
    489:       PetscCall(MatAssemblyEnd(fetidpmat_ctx->B_Ddelta, MAT_FINAL_ASSEMBLY));
    490:     } else {
    491:       /* scaling as in Klawonn-Widlund 1999 */
    492:       PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
    493:       PCBDDCSubSchurs     sub_schurs = pcbddc->sub_schurs;
    494:       Mat                 T;
    495:       PetscScalar        *W, lwork, *Bwork;
    496:       const PetscInt     *idxs = NULL;
    497:       PetscInt            cum, mss, *nnz;
    498:       PetscBLASInt       *pivots, B_lwork, B_N, B_ierr;

    500:       PetscCheck(pcbddc->deluxe_singlemat, comm, PETSC_ERR_USER, "Cannot compute B_Ddelta! rerun with -pc_bddc_deluxe_singlemat");
    501:       mss = 0;
    502:       PetscCall(PetscCalloc1(pcis->n_B, &nnz));
    503:       if (sub_schurs->is_Ej_all) {
    504:         PetscCall(ISGetIndices(sub_schurs->is_Ej_all, &idxs));
    505:         for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
    506:           PetscInt subset_size;

    508:           PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
    509:           for (j = 0; j < subset_size; j++) nnz[idxs[j + cum]] = subset_size;
    510:           mss = PetscMax(mss, subset_size);
    511:           cum += subset_size;
    512:         }
    513:       }
    514:       PetscCall(MatCreate(PETSC_COMM_SELF, &T));
    515:       PetscCall(MatSetSizes(T, pcis->n_B, pcis->n_B, pcis->n_B, pcis->n_B));
    516:       PetscCall(MatSetType(T, MATSEQAIJ));
    517:       PetscCall(MatSeqAIJSetPreallocation(T, 0, nnz));
    518:       PetscCall(PetscFree(nnz));

    520:       /* workspace allocation */
    521:       B_lwork = 0;
    522:       if (mss) {
    523:         PetscScalar dummy = 1;

    525:         B_lwork = -1;
    526:         PetscCall(PetscBLASIntCast(mss, &B_N));
    527:         PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    528:         PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&B_N, &dummy, &B_N, &B_N, &lwork, &B_lwork, &B_ierr));
    529:         PetscCall(PetscFPTrapPop());
    530:         PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in query to GETRI Lapack routine %" PetscBLASInt_FMT, B_ierr);
    531:         PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(lwork), &B_lwork));
    532:       }
    533:       PetscCall(PetscMalloc3(mss * mss, &W, mss, &pivots, B_lwork, &Bwork));

    535:       for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
    536:         const PetscScalar *M;
    537:         PetscInt           subset_size;

    539:         PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
    540:         PetscCall(PetscBLASIntCast(subset_size, &B_N));
    541:         PetscCall(MatDenseGetArrayRead(deluxe_ctx->seq_mat[i], &M));
    542:         PetscCall(PetscArraycpy(W, M, subset_size * subset_size));
    543:         PetscCall(MatDenseRestoreArrayRead(deluxe_ctx->seq_mat[i], &M));
    544:         PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    545:         PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&B_N, &B_N, W, &B_N, pivots, &B_ierr));
    546:         PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRF Lapack routine %" PetscBLASInt_FMT, B_ierr);
    547:         PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&B_N, W, &B_N, pivots, Bwork, &B_lwork, &B_ierr));
    548:         PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRI Lapack routine %" PetscBLASInt_FMT, B_ierr);
    549:         PetscCall(PetscFPTrapPop());
    550:         /* silent static analyzer */
    551:         PetscCheck(idxs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "IDXS not present");
    552:         PetscCall(MatSetValues(T, subset_size, idxs + cum, subset_size, idxs + cum, W, INSERT_VALUES));
    553:         cum += subset_size;
    554:       }
    555:       PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
    556:       PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));
    557:       PetscCall(MatMatTransposeMult(T, fetidpmat_ctx->B_delta, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &BD1));
    558:       PetscCall(MatMatMult(fetidpmat_ctx->B_delta, BD1, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &BD2));
    559:       PetscCall(MatDestroy(&T));
    560:       PetscCall(PetscFree3(W, pivots, Bwork));
    561:       if (sub_schurs->is_Ej_all) PetscCall(ISRestoreIndices(sub_schurs->is_Ej_all, &idxs));
    562:     }
    563:   }
    564:   PetscCall(PetscFree(scaling_factors));
    565:   PetscCall(PetscFree(cols_B_delta));

    567:   /* Layout of multipliers */
    568:   PetscCall(PetscLayoutCreate(comm, &llay));
    569:   PetscCall(PetscLayoutSetBlockSize(llay, 1));
    570:   PetscCall(PetscLayoutSetSize(llay, fetidpmat_ctx->n_lambda));
    571:   PetscCall(PetscLayoutSetUp(llay));
    572:   PetscCall(PetscLayoutGetLocalSize(llay, &fetidpmat_ctx->n));

    574:   /* Local work vector of multipliers */
    575:   PetscCall(VecCreate(PETSC_COMM_SELF, &fetidpmat_ctx->lambda_local));
    576:   PetscCall(VecSetSizes(fetidpmat_ctx->lambda_local, n_local_lambda, n_local_lambda));
    577:   PetscCall(VecSetType(fetidpmat_ctx->lambda_local, VECSEQ));

    579:   if (BD2) {
    580:     ISLocalToGlobalMapping l2g;
    581:     Mat                    T, TA, *pT;
    582:     IS                     is;
    583:     PetscInt               nl, N;
    584:     BDdelta_DN             ctx;

    586:     PetscCall(PetscLayoutGetLocalSize(llay, &nl));
    587:     PetscCall(PetscLayoutGetSize(llay, &N));
    588:     PetscCall(MatCreate(comm, &T));
    589:     PetscCall(MatSetSizes(T, nl, nl, N, N));
    590:     PetscCall(MatSetType(T, MATIS));
    591:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &l2g));
    592:     PetscCall(MatSetLocalToGlobalMapping(T, l2g, l2g));
    593:     PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
    594:     PetscCall(MatISSetLocalMat(T, BD2));
    595:     PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
    596:     PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));
    597:     PetscCall(MatDestroy(&BD2));
    598:     PetscCall(MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &TA));
    599:     PetscCall(MatDestroy(&T));
    600:     PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_USE_POINTER, &is));
    601:     PetscCall(MatCreateSubMatrices(TA, 1, &is, &is, MAT_INITIAL_MATRIX, &pT));
    602:     PetscCall(MatDestroy(&TA));
    603:     PetscCall(ISDestroy(&is));
    604:     BD2 = pT[0];
    605:     PetscCall(PetscFree(pT));

    607:     /* B_Ddelta for non-redundant multipliers with deluxe scaling */
    608:     PetscCall(PetscNew(&ctx));
    609:     PetscCall(MatSetType(fetidpmat_ctx->B_Ddelta, MATSHELL));
    610:     PetscCall(MatShellSetContext(fetidpmat_ctx->B_Ddelta, ctx));
    611:     PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT, (PetscErrorCodeFn *)MatMult_BDdelta_deluxe_nonred));
    612:     PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_BDdelta_deluxe_nonred));
    613:     PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_BDdelta_deluxe_nonred));
    614:     PetscCall(MatSetUp(fetidpmat_ctx->B_Ddelta));

    616:     PetscCall(PetscObjectReference((PetscObject)BD1));
    617:     ctx->BD = BD1;
    618:     PetscCall(KSPCreate(PETSC_COMM_SELF, &ctx->kBD));
    619:     PetscCall(KSPSetNestLevel(ctx->kBD, fetidpmat_ctx->pc->kspnestlevel));
    620:     PetscCall(KSPSetOperators(ctx->kBD, BD2, BD2));
    621:     PetscCall(VecDuplicate(fetidpmat_ctx->lambda_local, &ctx->work));
    622:     fetidpmat_ctx->deluxe_nonred = PETSC_TRUE;
    623:   }
    624:   PetscCall(MatDestroy(&BD1));
    625:   PetscCall(MatDestroy(&BD2));

    627:   /* fetidpmat sizes */
    628:   fetidpmat_ctx->n += nPgl;
    629:   fetidpmat_ctx->N = fetidpmat_ctx->n_lambda + nPg;

    631:   /* Global vector for FETI-DP linear system */
    632:   PetscCall(VecCreate(comm, &fetidp_global));
    633:   PetscCall(VecSetSizes(fetidp_global, fetidpmat_ctx->n, fetidpmat_ctx->N));
    634:   PetscCall(VecSetType(fetidp_global, VECMPI));
    635:   PetscCall(VecSetUp(fetidp_global));

    637:   /* Decide layout for fetidp dofs: if it is a saddle point problem
    638:      pressure is ordered first in the local part of the global vector
    639:      of the FETI-DP linear system */
    640:   if (nPg) {
    641:     Vec             v;
    642:     IS              IS_l2g_p, ais;
    643:     PetscLayout     alay;
    644:     const PetscInt *idxs, *pranges, *aranges, *lranges;
    645:     PetscInt       *l2g_indices_p, rst;
    646:     PetscMPIInt     rank;

    648:     PetscCall(PetscMalloc1(nPl, &l2g_indices_p));
    649:     PetscCall(VecGetLayout(fetidp_global, &alay));
    650:     PetscCall(PetscLayoutGetRanges(alay, &aranges));
    651:     PetscCall(PetscLayoutGetRanges(play, &pranges));
    652:     PetscCall(PetscLayoutGetRanges(llay, &lranges));

    654:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fetidp_global), &rank));
    655:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), pranges[rank + 1] - pranges[rank], aranges[rank], 1, &fetidpmat_ctx->pressure));
    656:     PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->pressure, "F_P"));
    657:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), lranges[rank + 1] - lranges[rank], aranges[rank] + pranges[rank + 1] - pranges[rank], 1, &fetidpmat_ctx->lagrange));
    658:     PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->lagrange, "F_L"));
    659:     PetscCall(ISLocalToGlobalMappingGetIndices(l2gmap_p, &idxs));
    660:     /* shift local to global indices for pressure */
    661:     for (i = 0; i < nPl; i++) {
    662:       PetscMPIInt owner;

    664:       PetscCall(PetscLayoutFindOwner(play, idxs[i], &owner));
    665:       l2g_indices_p[i] = idxs[i] - pranges[owner] + aranges[owner];
    666:     }
    667:     PetscCall(ISLocalToGlobalMappingRestoreIndices(l2gmap_p, &idxs));
    668:     PetscCall(ISCreateGeneral(comm, nPl, l2g_indices_p, PETSC_OWN_POINTER, &IS_l2g_p));

    670:     /* local to global scatter for pressure */
    671:     PetscCall(VecScatterCreate(fetidpmat_ctx->vP, NULL, fetidp_global, IS_l2g_p, &fetidpmat_ctx->l2g_p));
    672:     PetscCall(ISDestroy(&IS_l2g_p));

    674:     /* scatter for lagrange multipliers only */
    675:     PetscCall(VecCreate(comm, &v));
    676:     PetscCall(VecSetType(v, VECSTANDARD));
    677:     PetscCall(VecSetLayout(v, llay));
    678:     PetscCall(VecSetUp(v));
    679:     PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &ais));
    680:     PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, v, ais, &fetidpmat_ctx->l2g_lambda_only));
    681:     PetscCall(ISDestroy(&ais));
    682:     PetscCall(VecDestroy(&v));

    684:     /* shift local to global indices for multipliers */
    685:     for (i = 0; i < n_local_lambda; i++) {
    686:       PetscInt    ps;
    687:       PetscMPIInt owner;

    689:       PetscCall(PetscLayoutFindOwner(llay, l2g_indices[i], &owner));
    690:       ps             = pranges[owner + 1] - pranges[owner];
    691:       l2g_indices[i] = l2g_indices[i] - lranges[owner] + aranges[owner] + ps;
    692:     }

    694:     /* scatter from alldofs to pressures global fetidp vector */
    695:     PetscCall(PetscLayoutGetRange(alay, &rst, NULL));
    696:     PetscCall(ISCreateStride(comm, nPgl, rst, 1, &ais));
    697:     PetscCall(VecScatterCreate(pcis->vec1_global, pP, fetidp_global, ais, &fetidpmat_ctx->g2g_p));
    698:     PetscCall(ISDestroy(&ais));
    699:   }
    700:   PetscCall(PetscLayoutDestroy(&llay));
    701:   PetscCall(PetscLayoutDestroy(&play));
    702:   PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_OWN_POINTER, &IS_l2g_lambda));

    704:   /* scatter from local to global multipliers */
    705:   PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, fetidp_global, IS_l2g_lambda, &fetidpmat_ctx->l2g_lambda));
    706:   PetscCall(ISDestroy(&IS_l2g_lambda));
    707:   PetscCall(ISLocalToGlobalMappingDestroy(&l2gmap_p));
    708:   PetscCall(VecDestroy(&fetidp_global));

    710:   /* Create some work vectors needed by fetidp */
    711:   PetscCall(VecDuplicate(pcis->vec1_B, &fetidpmat_ctx->temp_solution_B));
    712:   PetscCall(VecDuplicate(pcis->vec1_D, &fetidpmat_ctx->temp_solution_D));
    713:   PetscFunctionReturn(PETSC_SUCCESS);
    714: }

    716: PetscErrorCode PCBDDCSetupFETIDPPCContext(Mat fetimat, FETIDPPC_ctx fetidppc_ctx)
    717: {
    718:   FETIDPMat_ctx mat_ctx;
    719:   PC_BDDC      *pcbddc = (PC_BDDC *)fetidppc_ctx->pc->data;
    720:   PC_IS        *pcis   = (PC_IS *)fetidppc_ctx->pc->data;
    721:   PetscBool     lumped = PETSC_FALSE;

    723:   PetscFunctionBegin;
    724:   PetscCall(MatShellGetContext(fetimat, &mat_ctx));
    725:   /* get references from objects created when setting up feti mat context */
    726:   PetscCall(PetscObjectReference((PetscObject)mat_ctx->lambda_local));
    727:   fetidppc_ctx->lambda_local = mat_ctx->lambda_local;
    728:   PetscCall(PetscObjectReference((PetscObject)mat_ctx->B_Ddelta));
    729:   fetidppc_ctx->B_Ddelta = mat_ctx->B_Ddelta;
    730:   if (mat_ctx->deluxe_nonred) {
    731:     PC            pc, mpc;
    732:     BDdelta_DN    ctx;
    733:     MatSolverType solver;
    734:     const char   *prefix;

    736:     PetscCall(MatShellGetContext(mat_ctx->B_Ddelta, &ctx));
    737:     PetscCall(KSPSetType(ctx->kBD, KSPPREONLY));
    738:     PetscCall(KSPGetPC(ctx->kBD, &mpc));
    739:     PetscCall(KSPGetPC(pcbddc->ksp_D, &pc));
    740:     PetscCall(PCSetType(mpc, PCLU));
    741:     PetscCall(PCFactorGetMatSolverType(pc, &solver));
    742:     if (solver) PetscCall(PCFactorSetMatSolverType(mpc, solver));
    743:     PetscCall(MatGetOptionsPrefix(fetimat, &prefix));
    744:     PetscCall(KSPSetOptionsPrefix(ctx->kBD, prefix));
    745:     PetscCall(KSPAppendOptionsPrefix(ctx->kBD, "bddelta_"));
    746:     PetscCall(KSPSetFromOptions(ctx->kBD));
    747:   }

    749:   if (mat_ctx->l2g_lambda_only) {
    750:     PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda_only));
    751:     fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda_only;
    752:   } else {
    753:     PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda));
    754:     fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda;
    755:   }
    756:   /* Dirichlet preconditioner */
    757:   PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_lumped", &lumped, NULL));
    758:   if (!lumped) {
    759:     IS        iV;
    760:     PetscBool discrete_harmonic = PETSC_FALSE;

    762:     PetscCall(PetscObjectQuery((PetscObject)fetidppc_ctx->pc, "__KSPFETIDP_iV", (PetscObject *)&iV));
    763:     if (iV) PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_discrete_harmonic", &discrete_harmonic, NULL));
    764:     if (discrete_harmonic) {
    765:       KSP             sksp;
    766:       PC              pc;
    767:       PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs;
    768:       Mat             A_II, A_IB, A_BI;
    769:       IS              iP = NULL;
    770:       PetscBool       isshell, reuse = PETSC_FALSE;
    771:       KSPType         ksptype;
    772:       const char     *prefix;

    774:       /*
    775:         We constructs a Schur complement for

    777:         | A_II A_ID |
    778:         | A_DI A_DD |

    780:         instead of

    782:         | A_II  B^t_II A_ID |
    783:         | B_II -C_II   B_ID |
    784:         | A_DI  B^t_ID A_DD |

    786:       */
    787:       if (sub_schurs && sub_schurs->reuse_solver) {
    788:         PetscCall(PetscObjectQuery((PetscObject)sub_schurs->A, "__KSPFETIDP_iP", (PetscObject *)&iP));
    789:         if (iP) reuse = PETSC_TRUE;
    790:       }
    791:       if (!reuse) {
    792:         IS       aB;
    793:         PetscInt nb;
    794:         PetscCall(ISGetLocalSize(pcis->is_B_local, &nb));
    795:         PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pcis->A_II), nb, 0, 1, &aB));
    796:         PetscCall(MatCreateSubMatrix(pcis->A_II, iV, iV, MAT_INITIAL_MATRIX, &A_II));
    797:         PetscCall(MatCreateSubMatrix(pcis->A_IB, iV, aB, MAT_INITIAL_MATRIX, &A_IB));
    798:         PetscCall(MatCreateSubMatrix(pcis->A_BI, aB, iV, MAT_INITIAL_MATRIX, &A_BI));
    799:         PetscCall(ISDestroy(&aB));
    800:       } else {
    801:         PetscCall(MatCreateSubMatrix(sub_schurs->A, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &A_IB));
    802:         PetscCall(MatCreateSubMatrix(sub_schurs->A, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &A_BI));
    803:         PetscCall(PetscObjectReference((PetscObject)pcis->A_II));
    804:         A_II = pcis->A_II;
    805:       }
    806:       PetscCall(MatCreateSchurComplement(A_II, A_II, A_IB, A_BI, pcis->A_BB, &fetidppc_ctx->S_j));

    808:       /* propagate settings of solver */
    809:       PetscCall(MatSchurComplementGetKSP(fetidppc_ctx->S_j, &sksp));
    810:       PetscCall(KSPGetType(pcis->ksp_D, &ksptype));
    811:       PetscCall(KSPSetType(sksp, ksptype));
    812:       PetscCall(KSPGetPC(pcis->ksp_D, &pc));
    813:       PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCSHELL, &isshell));
    814:       if (!isshell) {
    815:         MatSolverType solver;
    816:         PCType        pctype;

    818:         PetscCall(PCGetType(pc, &pctype));
    819:         PetscCall(PCFactorGetMatSolverType(pc, &solver));
    820:         PetscCall(KSPGetPC(sksp, &pc));
    821:         PetscCall(PCSetType(pc, pctype));
    822:         if (solver) PetscCall(PCFactorSetMatSolverType(pc, solver));
    823:       } else {
    824:         PetscCall(KSPGetPC(sksp, &pc));
    825:         PetscCall(PCSetType(pc, PCLU));
    826:       }
    827:       PetscCall(MatDestroy(&A_II));
    828:       PetscCall(MatDestroy(&A_IB));
    829:       PetscCall(MatDestroy(&A_BI));
    830:       PetscCall(MatGetOptionsPrefix(fetimat, &prefix));
    831:       PetscCall(KSPSetOptionsPrefix(sksp, prefix));
    832:       PetscCall(KSPAppendOptionsPrefix(sksp, "harmonic_"));
    833:       PetscCall(KSPSetFromOptions(sksp));
    834:       if (reuse) {
    835:         PetscCall(KSPSetPC(sksp, sub_schurs->reuse_solver->interior_solver));
    836:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)sub_schurs->reuse_solver->interior_solver, (PetscObject)sksp, 0));
    837:       }
    838:     } else { /* default Dirichlet preconditioner is pde-harmonic */
    839:       PetscCall(MatCreateSchurComplement(pcis->A_II, pcis->A_II, pcis->A_IB, pcis->A_BI, pcis->A_BB, &fetidppc_ctx->S_j));
    840:       PetscCall(MatSchurComplementSetKSP(fetidppc_ctx->S_j, pcis->ksp_D));
    841:     }
    842:   } else {
    843:     PetscCall(PetscObjectReference((PetscObject)pcis->A_BB));
    844:     fetidppc_ctx->S_j = pcis->A_BB;
    845:   }
    846:   /* saddle-point */
    847:   if (mat_ctx->xPg) {
    848:     PetscCall(PetscObjectReference((PetscObject)mat_ctx->xPg));
    849:     fetidppc_ctx->xPg = mat_ctx->xPg;
    850:     PetscCall(PetscObjectReference((PetscObject)mat_ctx->yPg));
    851:     fetidppc_ctx->yPg = mat_ctx->yPg;
    852:   }
    853:   PetscFunctionReturn(PETSC_SUCCESS);
    854: }

    856: static PetscErrorCode FETIDPMatMult_Kernel(Mat fetimat, Vec x, Vec y, PetscBool trans)
    857: {
    858:   FETIDPMat_ctx mat_ctx;
    859:   PC_BDDC      *pcbddc;
    860:   PC_IS        *pcis;

    862:   PetscFunctionBegin;
    863:   PetscCall(MatShellGetContext(fetimat, &mat_ctx));
    864:   pcis   = (PC_IS *)mat_ctx->pc->data;
    865:   pcbddc = (PC_BDDC *)mat_ctx->pc->data;
    866:   /* Application of B_delta^T */
    867:   PetscCall(VecSet(pcis->vec1_B, 0.));
    868:   PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    869:   PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    870:   PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B));

    872:   /* Add contribution from saddle point */
    873:   if (mat_ctx->l2g_p) {
    874:     PetscCall(VecScatterBegin(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
    875:     PetscCall(VecScatterEnd(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
    876:     if (pcbddc->switch_static) {
    877:       if (trans) {
    878:         PetscCall(MatMultTranspose(mat_ctx->B_BI, mat_ctx->vP, pcis->vec1_D));
    879:       } else {
    880:         PetscCall(MatMult(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D));
    881:       }
    882:     }
    883:     if (trans) {
    884:       PetscCall(MatMultTransposeAdd(mat_ctx->B_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
    885:     } else {
    886:       PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
    887:     }
    888:   } else {
    889:     if (pcbddc->switch_static) PetscCall(VecSet(pcis->vec1_D, 0.0));
    890:   }
    891:   /* Application of \widetilde{S}^-1 */
    892:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
    893:   PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, trans));
    894:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
    895:   PetscCall(VecSet(y, 0.0));
    896:   /* Application of B_delta */
    897:   PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local));
    898:   /* Contribution from boundary pressures */
    899:   if (mat_ctx->C) {
    900:     const PetscScalar *lx;
    901:     PetscScalar       *ly;

    903:     /* pressure ordered first in the local part of x and y */
    904:     PetscCall(VecGetArrayRead(x, &lx));
    905:     PetscCall(VecGetArray(y, &ly));
    906:     PetscCall(VecPlaceArray(mat_ctx->xPg, lx));
    907:     PetscCall(VecPlaceArray(mat_ctx->yPg, ly));
    908:     if (trans) {
    909:       PetscCall(MatMultTranspose(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg));
    910:     } else {
    911:       PetscCall(MatMult(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg));
    912:     }
    913:     PetscCall(VecResetArray(mat_ctx->xPg));
    914:     PetscCall(VecResetArray(mat_ctx->yPg));
    915:     PetscCall(VecRestoreArrayRead(x, &lx));
    916:     PetscCall(VecRestoreArray(y, &ly));
    917:   }
    918:   /* Add contribution from saddle point */
    919:   if (mat_ctx->l2g_p) {
    920:     PetscCall(VecISSet(pcis->vec1_B, mat_ctx->lP_B, 0));
    921:     if (trans) {
    922:       PetscCall(MatMultTranspose(mat_ctx->Bt_BB, pcis->vec1_B, mat_ctx->vP));
    923:     } else {
    924:       PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP));
    925:     }
    926:     if (pcbddc->switch_static) {
    927:       PetscCall(VecISSet(pcis->vec1_D, mat_ctx->lP_I, 0));
    928:       if (trans) {
    929:         PetscCall(MatMultTransposeAdd(mat_ctx->Bt_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
    930:       } else {
    931:         PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
    932:       }
    933:     }
    934:     PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD));
    935:     PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD));
    936:   }
    937:   PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD));
    938:   PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD));
    939:   PetscFunctionReturn(PETSC_SUCCESS);
    940: }

    942: PetscErrorCode FETIDPMatMult(Mat fetimat, Vec x, Vec y)
    943: {
    944:   PetscFunctionBegin;
    945:   PetscCall(FETIDPMatMult_Kernel(fetimat, x, y, PETSC_FALSE));
    946:   PetscFunctionReturn(PETSC_SUCCESS);
    947: }

    949: PetscErrorCode FETIDPMatMultTranspose(Mat fetimat, Vec x, Vec y)
    950: {
    951:   PetscFunctionBegin;
    952:   PetscCall(FETIDPMatMult_Kernel(fetimat, x, y, PETSC_TRUE));
    953:   PetscFunctionReturn(PETSC_SUCCESS);
    954: }

    956: static PetscErrorCode FETIDPPCApply_Kernel(PC fetipc, Vec x, Vec y, PetscBool trans)
    957: {
    958:   FETIDPPC_ctx pc_ctx;
    959:   PC_IS       *pcis;

    961:   PetscFunctionBegin;
    962:   PetscCall(PCShellGetContext(fetipc, &pc_ctx));
    963:   pcis = (PC_IS *)pc_ctx->pc->data;
    964:   /* Application of B_Ddelta^T */
    965:   PetscCall(VecScatterBegin(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    966:   PetscCall(VecScatterEnd(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    967:   PetscCall(VecSet(pcis->vec2_B, 0.0));
    968:   PetscCall(MatMultTranspose(pc_ctx->B_Ddelta, pc_ctx->lambda_local, pcis->vec2_B));
    969:   /* Application of local Schur complement */
    970:   if (trans) {
    971:     PetscCall(MatMultTranspose(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B));
    972:   } else {
    973:     PetscCall(MatMult(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B));
    974:   }
    975:   /* Application of B_Ddelta */
    976:   PetscCall(MatMult(pc_ctx->B_Ddelta, pcis->vec1_B, pc_ctx->lambda_local));
    977:   PetscCall(VecSet(y, 0.0));
    978:   PetscCall(VecScatterBegin(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD));
    979:   PetscCall(VecScatterEnd(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD));
    980:   PetscFunctionReturn(PETSC_SUCCESS);
    981: }

    983: PetscErrorCode FETIDPPCApply(PC pc, Vec x, Vec y)
    984: {
    985:   PetscFunctionBegin;
    986:   PetscCall(FETIDPPCApply_Kernel(pc, x, y, PETSC_FALSE));
    987:   PetscFunctionReturn(PETSC_SUCCESS);
    988: }

    990: PetscErrorCode FETIDPPCApplyTranspose(PC pc, Vec x, Vec y)
    991: {
    992:   PetscFunctionBegin;
    993:   PetscCall(FETIDPPCApply_Kernel(pc, x, y, PETSC_TRUE));
    994:   PetscFunctionReturn(PETSC_SUCCESS);
    995: }

    997: PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer)
    998: {
    999:   FETIDPPC_ctx pc_ctx;
    1000:   PetscBool    isascii;
    1001:   PetscViewer  sviewer;

    1003:   PetscFunctionBegin;
    1004:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    1005:   if (isascii) {
    1006:     PetscMPIInt rank;
    1007:     PetscBool   isschur, isshell;

    1009:     PetscCall(PCShellGetContext(pc, &pc_ctx));
    1010:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
    1011:     PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->S_j, MATSCHURCOMPLEMENT, &isschur));
    1012:     if (isschur) {
    1013:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Dirichlet preconditioner (just from rank 0)\n"));
    1014:     } else {
    1015:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Lumped preconditioner (just from rank 0)\n"));
    1016:     }
    1017:     PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer));
    1018:     if (rank == 0) {
    1019:       PetscCall(PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO));
    1020:       PetscCall(PetscViewerASCIIPushTab(sviewer));
    1021:       PetscCall(MatView(pc_ctx->S_j, sviewer));
    1022:       PetscCall(PetscViewerASCIIPopTab(sviewer));
    1023:       PetscCall(PetscViewerPopFormat(sviewer));
    1024:     }
    1025:     PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer));
    1026:     PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->B_Ddelta, MATSHELL, &isshell));
    1027:     if (isshell) {
    1028:       BDdelta_DN ctx;
    1029:       PetscCall(PetscViewerASCIIPrintf(viewer, "  FETI-DP BDdelta: DB^t * (B D^-1 B^t)^-1 for deluxe scaling (just from rank 0)\n"));
    1030:       PetscCall(MatShellGetContext(pc_ctx->B_Ddelta, &ctx));
    1031:       PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer));
    1032:       if (rank == 0) {
    1033:         PetscInt tl;

    1035:         PetscCall(PetscViewerASCIIGetTab(sviewer, &tl));
    1036:         PetscCall(PetscObjectSetTabLevel((PetscObject)ctx->kBD, tl));
    1037:         PetscCall(KSPView(ctx->kBD, sviewer));
    1038:         PetscCall(PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO));
    1039:         PetscCall(MatView(ctx->BD, sviewer));
    1040:         PetscCall(PetscViewerPopFormat(sviewer));
    1041:       }
    1042:       PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer));
    1043:     }
    1044:     PetscCall(PetscViewerFlush(viewer));
    1045:   }
    1046:   PetscFunctionReturn(PETSC_SUCCESS);
    1047: }