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(MatDestroy(&mat_ctx->B_BB));
 80:   PetscCall(MatDestroy(&mat_ctx->B_BI));
 81:   PetscCall(MatDestroy(&mat_ctx->Bt_BB));
 82:   PetscCall(MatDestroy(&mat_ctx->Bt_BI));
 83:   PetscCall(MatDestroy(&mat_ctx->C));
 84:   PetscCall(VecDestroy(&mat_ctx->rhs_flip));
 85:   PetscCall(VecDestroy(&mat_ctx->vP));
 86:   PetscCall(VecDestroy(&mat_ctx->xPg));
 87:   PetscCall(VecDestroy(&mat_ctx->yPg));
 88:   PetscCall(VecScatterDestroy(&mat_ctx->l2g_lambda));
 89:   PetscCall(VecScatterDestroy(&mat_ctx->l2g_lambda_only));
 90:   PetscCall(VecScatterDestroy(&mat_ctx->l2g_p));
 91:   PetscCall(VecScatterDestroy(&mat_ctx->g2g_p));
 92:   PetscCall(PCDestroy(&mat_ctx->pc)); /* decrease PCBDDC reference count */
 93:   PetscCall(ISDestroy(&mat_ctx->pressure));
 94:   PetscCall(ISDestroy(&mat_ctx->lagrange));
 95:   PetscCall(PetscFree(mat_ctx));
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

 99: PetscErrorCode PCBDDCDestroyFETIDPPC(PC pc)
100: {
101:   FETIDPPC_ctx pc_ctx;

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

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

141:   /* saddlepoint */
142:   ISLocalToGlobalMapping l2gmap_p;
143:   PetscLayout            play;
144:   IS                     gP, pP;
145:   PetscInt               nPl, nPg, nPgl;

147:   PetscFunctionBegin;
148:   PetscCall(PetscObjectGetComm((PetscObject)(fetidpmat_ctx->pc), &comm));
149:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
150:   PetscCallMPI(MPI_Comm_size(comm, &size));

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

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

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

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

202:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_B_BB", (PetscObject *)&fetidpmat_ctx->B_BB));
203:     PetscCheck(fetidpmat_ctx->B_BB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "B_BB not present");
204:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->B_BB));

206:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_Bt_BI", (PetscObject *)&fetidpmat_ctx->Bt_BI));
207:     PetscCheck(fetidpmat_ctx->Bt_BI, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bt_BI not present");
208:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BI));

210:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_Bt_BB", (PetscObject *)&fetidpmat_ctx->Bt_BB));
211:     PetscCheck(fetidpmat_ctx->Bt_BB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bt_BB not present");
212:     PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BB));

214:     PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_flip", (PetscObject *)&fetidpmat_ctx->rhs_flip));
215:     if (fetidpmat_ctx->rhs_flip) PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->rhs_flip));
216:   }

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

221:   /* Evaluate local and global number of lagrange multipliers */
222:   PetscCall(VecSet(pcis->vec1_N, 0.0));
223:   n_local_lambda  = 0;
224:   partial_sum     = 0;
225:   n_boundary_dofs = 0;
226:   s               = 0;

228:   /* Get Vertices used to define the BDDC */
229:   PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &isvert));
230:   PetscCall(ISGetLocalSize(isvert, &n_vertices));
231:   PetscCall(ISGetIndices(isvert, &vertex_indices));

233:   dual_size = pcis->n_B - n_vertices;
234:   PetscCall(PetscMalloc1(dual_size, &dual_dofs_boundary_indices));
235:   PetscCall(PetscMalloc1(dual_size, &aux_local_numbering_1));
236:   PetscCall(PetscMalloc1(dual_size, &aux_local_numbering_2));

238:   PetscCall(VecGetArray(pcis->vec1_N, &array));
239:   for (i = 0; i < pcis->n; i++) {
240:     j = mat_graph->count[i]; /* RECALL: mat_graph->count[i] does not count myself */
241:     if (j > 0) n_boundary_dofs++;
242:     skip_node = PETSC_FALSE;
243:     if (s < n_vertices && vertex_indices[s] == i) { /* it works for a sorted set of vertices */
244:       skip_node = PETSC_TRUE;
245:       s++;
246:     }
247:     if (j < 1) skip_node = PETSC_TRUE;
248:     if (mat_graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) skip_node = PETSC_TRUE;
249:     if (!skip_node) {
250:       if (fully_redundant) {
251:         /* fully redundant set of lagrange multipliers */
252:         n_lambda_for_dof = (j * (j + 1)) / 2;
253:       } else {
254:         n_lambda_for_dof = j;
255:       }
256:       n_local_lambda += j;
257:       /* needed to evaluate global number of lagrange multipliers */
258:       array[i] = (1.0 * n_lambda_for_dof) / (j + 1.0); /* already scaled for the next global sum */
259:       /* store some data needed */
260:       dual_dofs_boundary_indices[partial_sum] = n_boundary_dofs - 1;
261:       aux_local_numbering_1[partial_sum]      = i;
262:       aux_local_numbering_2[partial_sum]      = n_lambda_for_dof;
263:       partial_sum++;
264:     }
265:   }
266:   PetscCall(VecRestoreArray(pcis->vec1_N, &array));
267:   PetscCall(ISRestoreIndices(isvert, &vertex_indices));
268:   PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &isvert));
269:   dual_size = partial_sum;

271:   /* compute global ordering of lagrange multipliers and associate l2g map */
272:   PetscCall(ISCreateGeneral(comm, partial_sum, aux_local_numbering_1, PETSC_COPY_VALUES, &subset_n));
273:   PetscCall(ISLocalToGlobalMappingApplyIS(pcis->mapping, subset_n, &subset));
274:   PetscCall(ISDestroy(&subset_n));
275:   PetscCall(ISCreateGeneral(comm, partial_sum, aux_local_numbering_2, PETSC_OWN_POINTER, &subset_mult));
276:   PetscCall(ISRenumber(subset, subset_mult, &fetidpmat_ctx->n_lambda, &subset_n));
277:   PetscCall(ISDestroy(&subset));

279:   if (PetscDefined(USE_DEBUG)) {
280:     PetscCall(VecSet(pcis->vec1_global, 0.0));
281:     PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
282:     PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
283:     PetscCall(VecSum(pcis->vec1_global, &scalar_value));
284:     i = (PetscInt)PetscRealPart(scalar_value);
285:     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);
286:   }

288:   /* init data for scaling factors exchange */
289:   if (!pcbddc->use_deluxe_scaling) {
290:     PetscInt    *ptrs_buffer, neigh_position;
291:     PetscScalar *send_buffer, *recv_buffer;
292:     MPI_Request *send_reqs, *recv_reqs;

294:     partial_sum = 0;
295:     PetscCall(PetscMalloc1(pcis->n_neigh, &ptrs_buffer));
296:     PetscCall(PetscMalloc1(PetscMax(pcis->n_neigh - 1, 0), &send_reqs));
297:     PetscCall(PetscMalloc1(PetscMax(pcis->n_neigh - 1, 0), &recv_reqs));
298:     PetscCall(PetscMalloc1(pcis->n + 1, &all_factors));
299:     if (pcis->n_neigh > 0) ptrs_buffer[0] = 0;
300:     for (i = 1; i < pcis->n_neigh; i++) {
301:       partial_sum += pcis->n_shared[i];
302:       ptrs_buffer[i] = ptrs_buffer[i - 1] + pcis->n_shared[i];
303:     }
304:     PetscCall(PetscMalloc1(partial_sum, &send_buffer));
305:     PetscCall(PetscMalloc1(partial_sum, &recv_buffer));
306:     PetscCall(PetscMalloc1(partial_sum, &all_factors[0]));
307:     for (i = 0; i < pcis->n - 1; i++) {
308:       j                  = mat_graph->count[i];
309:       all_factors[i + 1] = all_factors[i] + j;
310:     }

312:     /* scatter B scaling to N vec */
313:     PetscCall(VecScatterBegin(pcis->N_to_B, pcis->D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
314:     PetscCall(VecScatterEnd(pcis->N_to_B, pcis->D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
315:     /* communications */
316:     PetscCall(VecGetArrayRead(pcis->vec1_N, (const PetscScalar **)&array));
317:     for (i = 1; i < pcis->n_neigh; i++) {
318:       for (j = 0; j < pcis->n_shared[i]; j++) send_buffer[ptrs_buffer[i - 1] + j] = array[pcis->shared[i][j]];
319:       PetscCall(PetscMPIIntCast(ptrs_buffer[i] - ptrs_buffer[i - 1], &buf_size));
320:       PetscCall(PetscMPIIntCast(pcis->neigh[i], &neigh));
321:       PetscCallMPI(MPI_Isend(&send_buffer[ptrs_buffer[i - 1]], buf_size, MPIU_SCALAR, neigh, 0, comm, &send_reqs[i - 1]));
322:       PetscCallMPI(MPI_Irecv(&recv_buffer[ptrs_buffer[i - 1]], buf_size, MPIU_SCALAR, neigh, 0, comm, &recv_reqs[i - 1]));
323:     }
324:     PetscCall(VecRestoreArrayRead(pcis->vec1_N, (const PetscScalar **)&array));
325:     if (pcis->n_neigh > 0) PetscCallMPI(MPI_Waitall(pcis->n_neigh - 1, recv_reqs, MPI_STATUSES_IGNORE));
326:     /* put values in correct places */
327:     for (i = 1; i < pcis->n_neigh; i++) {
328:       for (j = 0; j < pcis->n_shared[i]; j++) {
329:         k              = pcis->shared[i][j];
330:         neigh_position = 0;
331:         while (mat_graph->neighbours_set[k][neigh_position] != pcis->neigh[i]) neigh_position++;
332:         all_factors[k][neigh_position] = recv_buffer[ptrs_buffer[i - 1] + j];
333:       }
334:     }
335:     if (pcis->n_neigh > 0) PetscCallMPI(MPI_Waitall(pcis->n_neigh - 1, send_reqs, MPI_STATUSES_IGNORE));
336:     PetscCall(PetscFree(send_reqs));
337:     PetscCall(PetscFree(recv_reqs));
338:     PetscCall(PetscFree(send_buffer));
339:     PetscCall(PetscFree(recv_buffer));
340:     PetscCall(PetscFree(ptrs_buffer));
341:   }

343:   /* Compute B and B_delta (local actions) */
344:   PetscCall(PetscMalloc1(pcis->n_neigh, &aux_sums));
345:   PetscCall(PetscMalloc1(n_local_lambda, &l2g_indices));
346:   PetscCall(PetscMalloc1(n_local_lambda, &vals_B_delta));
347:   PetscCall(PetscMalloc1(n_local_lambda, &cols_B_delta));
348:   if (!pcbddc->use_deluxe_scaling) {
349:     PetscCall(PetscMalloc1(n_local_lambda, &scaling_factors));
350:   } else {
351:     scaling_factors = NULL;
352:     all_factors     = NULL;
353:   }
354:   PetscCall(ISGetIndices(subset_n, &aux_global_numbering));
355:   partial_sum = 0;
356:   cum         = 0;
357:   for (i = 0; i < dual_size; i++) {
358:     n_global_lambda = aux_global_numbering[cum];
359:     j               = mat_graph->count[aux_local_numbering_1[i]];
360:     aux_sums[0]     = 0;
361:     for (s = 1; s < j; s++) aux_sums[s] = aux_sums[s - 1] + j - s + 1;
362:     if (all_factors) array = all_factors[aux_local_numbering_1[i]];
363:     n_neg_values = 0;
364:     while (n_neg_values < j && mat_graph->neighbours_set[aux_local_numbering_1[i]][n_neg_values] < rank) n_neg_values++;
365:     n_pos_values = j - n_neg_values;
366:     if (fully_redundant) {
367:       for (s = 0; s < n_neg_values; s++) {
368:         l2g_indices[partial_sum + s]  = aux_sums[s] + n_neg_values - s - 1 + n_global_lambda;
369:         cols_B_delta[partial_sum + s] = dual_dofs_boundary_indices[i];
370:         vals_B_delta[partial_sum + s] = -1.0;
371:         if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum + s] = array[s];
372:       }
373:       for (s = 0; s < n_pos_values; s++) {
374:         l2g_indices[partial_sum + s + n_neg_values]  = aux_sums[n_neg_values] + s + n_global_lambda;
375:         cols_B_delta[partial_sum + s + n_neg_values] = dual_dofs_boundary_indices[i];
376:         vals_B_delta[partial_sum + s + n_neg_values] = 1.0;
377:         if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum + s + n_neg_values] = array[s + n_neg_values];
378:       }
379:       partial_sum += j;
380:     } else {
381:       /* l2g_indices and default cols and vals of B_delta */
382:       for (s = 0; s < j; s++) {
383:         l2g_indices[partial_sum + s]  = n_global_lambda + s;
384:         cols_B_delta[partial_sum + s] = dual_dofs_boundary_indices[i];
385:         vals_B_delta[partial_sum + s] = 0.0;
386:       }
387:       /* B_delta */
388:       if (n_neg_values > 0) { /* there's a rank next to me to the left */
389:         vals_B_delta[partial_sum + n_neg_values - 1] = -1.0;
390:       }
391:       if (n_neg_values < j) { /* there's a rank next to me to the right */
392:         vals_B_delta[partial_sum + n_neg_values] = 1.0;
393:       }
394:       /* scaling as in Klawonn-Widlund 1999 */
395:       if (!pcbddc->use_deluxe_scaling) {
396:         for (s = 0; s < n_neg_values; s++) {
397:           scalar_value = 0.0;
398:           for (k = 0; k < s + 1; k++) scalar_value += array[k];
399:           scaling_factors[partial_sum + s] = -scalar_value;
400:         }
401:         for (s = 0; s < n_pos_values; s++) {
402:           scalar_value = 0.0;
403:           for (k = s + n_neg_values; k < j; k++) scalar_value += array[k];
404:           scaling_factors[partial_sum + s + n_neg_values] = scalar_value;
405:         }
406:       }
407:       partial_sum += j;
408:     }
409:     cum += aux_local_numbering_2[i];
410:   }
411:   PetscCall(ISRestoreIndices(subset_n, &aux_global_numbering));
412:   PetscCall(ISDestroy(&subset_mult));
413:   PetscCall(ISDestroy(&subset_n));
414:   PetscCall(PetscFree(aux_sums));
415:   PetscCall(PetscFree(aux_local_numbering_1));
416:   PetscCall(PetscFree(dual_dofs_boundary_indices));
417:   if (all_factors) {
418:     PetscCall(PetscFree(all_factors[0]));
419:     PetscCall(PetscFree(all_factors));
420:   }

422:   /* Create local part of B_delta */
423:   PetscCall(MatCreate(PETSC_COMM_SELF, &fetidpmat_ctx->B_delta));
424:   PetscCall(MatSetSizes(fetidpmat_ctx->B_delta, n_local_lambda, pcis->n_B, n_local_lambda, pcis->n_B));
425:   PetscCall(MatSetType(fetidpmat_ctx->B_delta, MATSEQAIJ));
426:   PetscCall(MatSeqAIJSetPreallocation(fetidpmat_ctx->B_delta, 1, NULL));
427:   PetscCall(MatSetOption(fetidpmat_ctx->B_delta, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
428:   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));
429:   PetscCall(PetscFree(vals_B_delta));
430:   PetscCall(MatAssemblyBegin(fetidpmat_ctx->B_delta, MAT_FINAL_ASSEMBLY));
431:   PetscCall(MatAssemblyEnd(fetidpmat_ctx->B_delta, MAT_FINAL_ASSEMBLY));

433:   BD1 = NULL;
434:   BD2 = NULL;
435:   if (fully_redundant) {
436:     PetscCheck(!pcbddc->use_deluxe_scaling, comm, PETSC_ERR_SUP, "Deluxe FETIDP with fully-redundant multipliers to be implemented");
437:     PetscCall(MatCreate(PETSC_COMM_SELF, &ScalingMat));
438:     PetscCall(MatSetSizes(ScalingMat, n_local_lambda, n_local_lambda, n_local_lambda, n_local_lambda));
439:     PetscCall(MatSetType(ScalingMat, MATSEQAIJ));
440:     PetscCall(MatSeqAIJSetPreallocation(ScalingMat, 1, NULL));
441:     for (i = 0; i < n_local_lambda; i++) PetscCall(MatSetValue(ScalingMat, i, i, scaling_factors[i], INSERT_VALUES));
442:     PetscCall(MatAssemblyBegin(ScalingMat, MAT_FINAL_ASSEMBLY));
443:     PetscCall(MatAssemblyEnd(ScalingMat, MAT_FINAL_ASSEMBLY));
444:     PetscCall(MatMatMult(ScalingMat, fetidpmat_ctx->B_delta, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &fetidpmat_ctx->B_Ddelta));
445:     PetscCall(MatDestroy(&ScalingMat));
446:   } else {
447:     PetscCall(MatCreate(PETSC_COMM_SELF, &fetidpmat_ctx->B_Ddelta));
448:     PetscCall(MatSetSizes(fetidpmat_ctx->B_Ddelta, n_local_lambda, pcis->n_B, n_local_lambda, pcis->n_B));
449:     if (!pcbddc->use_deluxe_scaling || !pcbddc->sub_schurs) {
450:       PetscCall(MatSetType(fetidpmat_ctx->B_Ddelta, MATSEQAIJ));
451:       PetscCall(MatSeqAIJSetPreallocation(fetidpmat_ctx->B_Ddelta, 1, NULL));
452:       for (i = 0; i < n_local_lambda; i++) PetscCall(MatSetValue(fetidpmat_ctx->B_Ddelta, i, cols_B_delta[i], scaling_factors[i], INSERT_VALUES));
453:       PetscCall(MatAssemblyBegin(fetidpmat_ctx->B_Ddelta, MAT_FINAL_ASSEMBLY));
454:       PetscCall(MatAssemblyEnd(fetidpmat_ctx->B_Ddelta, MAT_FINAL_ASSEMBLY));
455:     } else {
456:       /* scaling as in Klawonn-Widlund 1999 */
457:       PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
458:       PCBDDCSubSchurs     sub_schurs = pcbddc->sub_schurs;
459:       Mat                 T;
460:       PetscScalar        *W, lwork, *Bwork;
461:       const PetscInt     *idxs = NULL;
462:       PetscInt            cum, mss, *nnz;
463:       PetscBLASInt       *pivots, B_lwork, B_N, B_ierr;

465:       PetscCheck(pcbddc->deluxe_singlemat, comm, PETSC_ERR_USER, "Cannot compute B_Ddelta! rerun with -pc_bddc_deluxe_singlemat");
466:       mss = 0;
467:       PetscCall(PetscCalloc1(pcis->n_B, &nnz));
468:       if (sub_schurs->is_Ej_all) {
469:         PetscCall(ISGetIndices(sub_schurs->is_Ej_all, &idxs));
470:         for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
471:           PetscInt subset_size;

473:           PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
474:           for (j = 0; j < subset_size; j++) nnz[idxs[j + cum]] = subset_size;
475:           mss = PetscMax(mss, subset_size);
476:           cum += subset_size;
477:         }
478:       }
479:       PetscCall(MatCreate(PETSC_COMM_SELF, &T));
480:       PetscCall(MatSetSizes(T, pcis->n_B, pcis->n_B, pcis->n_B, pcis->n_B));
481:       PetscCall(MatSetType(T, MATSEQAIJ));
482:       PetscCall(MatSeqAIJSetPreallocation(T, 0, nnz));
483:       PetscCall(PetscFree(nnz));

485:       /* workspace allocation */
486:       B_lwork = 0;
487:       if (mss) {
488:         PetscScalar dummy = 1;

490:         B_lwork = -1;
491:         PetscCall(PetscBLASIntCast(mss, &B_N));
492:         PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
493:         PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&B_N, &dummy, &B_N, &B_N, &lwork, &B_lwork, &B_ierr));
494:         PetscCall(PetscFPTrapPop());
495:         PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in query to GETRI Lapack routine %d", (int)B_ierr);
496:         PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(lwork), &B_lwork));
497:       }
498:       PetscCall(PetscMalloc3(mss * mss, &W, mss, &pivots, B_lwork, &Bwork));

500:       for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
501:         const PetscScalar *M;
502:         PetscInt           subset_size;

504:         PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
505:         PetscCall(PetscBLASIntCast(subset_size, &B_N));
506:         PetscCall(MatDenseGetArrayRead(deluxe_ctx->seq_mat[i], &M));
507:         PetscCall(PetscArraycpy(W, M, subset_size * subset_size));
508:         PetscCall(MatDenseRestoreArrayRead(deluxe_ctx->seq_mat[i], &M));
509:         PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
510:         PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&B_N, &B_N, W, &B_N, pivots, &B_ierr));
511:         PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRF Lapack routine %d", (int)B_ierr);
512:         PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&B_N, W, &B_N, pivots, Bwork, &B_lwork, &B_ierr));
513:         PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRI Lapack routine %d", (int)B_ierr);
514:         PetscCall(PetscFPTrapPop());
515:         /* silent static analyzer */
516:         PetscCheck(idxs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "IDXS not present");
517:         PetscCall(MatSetValues(T, subset_size, idxs + cum, subset_size, idxs + cum, W, INSERT_VALUES));
518:         cum += subset_size;
519:       }
520:       PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
521:       PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));
522:       PetscCall(MatMatTransposeMult(T, fetidpmat_ctx->B_delta, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &BD1));
523:       PetscCall(MatMatMult(fetidpmat_ctx->B_delta, BD1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &BD2));
524:       PetscCall(MatDestroy(&T));
525:       PetscCall(PetscFree3(W, pivots, Bwork));
526:       if (sub_schurs->is_Ej_all) PetscCall(ISRestoreIndices(sub_schurs->is_Ej_all, &idxs));
527:     }
528:   }
529:   PetscCall(PetscFree(scaling_factors));
530:   PetscCall(PetscFree(cols_B_delta));

532:   /* Layout of multipliers */
533:   PetscCall(PetscLayoutCreate(comm, &llay));
534:   PetscCall(PetscLayoutSetBlockSize(llay, 1));
535:   PetscCall(PetscLayoutSetSize(llay, fetidpmat_ctx->n_lambda));
536:   PetscCall(PetscLayoutSetUp(llay));
537:   PetscCall(PetscLayoutGetLocalSize(llay, &fetidpmat_ctx->n));

539:   /* Local work vector of multipliers */
540:   PetscCall(VecCreate(PETSC_COMM_SELF, &fetidpmat_ctx->lambda_local));
541:   PetscCall(VecSetSizes(fetidpmat_ctx->lambda_local, n_local_lambda, n_local_lambda));
542:   PetscCall(VecSetType(fetidpmat_ctx->lambda_local, VECSEQ));

544:   if (BD2) {
545:     ISLocalToGlobalMapping l2g;
546:     Mat                    T, TA, *pT;
547:     IS                     is;
548:     PetscInt               nl, N;
549:     BDdelta_DN             ctx;

551:     PetscCall(PetscLayoutGetLocalSize(llay, &nl));
552:     PetscCall(PetscLayoutGetSize(llay, &N));
553:     PetscCall(MatCreate(comm, &T));
554:     PetscCall(MatSetSizes(T, nl, nl, N, N));
555:     PetscCall(MatSetType(T, MATIS));
556:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &l2g));
557:     PetscCall(MatSetLocalToGlobalMapping(T, l2g, l2g));
558:     PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
559:     PetscCall(MatISSetLocalMat(T, BD2));
560:     PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
561:     PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));
562:     PetscCall(MatDestroy(&BD2));
563:     PetscCall(MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &TA));
564:     PetscCall(MatDestroy(&T));
565:     PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_USE_POINTER, &is));
566:     PetscCall(MatCreateSubMatrices(TA, 1, &is, &is, MAT_INITIAL_MATRIX, &pT));
567:     PetscCall(MatDestroy(&TA));
568:     PetscCall(ISDestroy(&is));
569:     BD2 = pT[0];
570:     PetscCall(PetscFree(pT));

572:     /* B_Ddelta for non-redundant multipliers with deluxe scaling */
573:     PetscCall(PetscNew(&ctx));
574:     PetscCall(MatSetType(fetidpmat_ctx->B_Ddelta, MATSHELL));
575:     PetscCall(MatShellSetContext(fetidpmat_ctx->B_Ddelta, ctx));
576:     PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT, (void (*)(void))MatMult_BDdelta_deluxe_nonred));
577:     PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_BDdelta_deluxe_nonred));
578:     PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_DESTROY, (void (*)(void))MatDestroy_BDdelta_deluxe_nonred));
579:     PetscCall(MatSetUp(fetidpmat_ctx->B_Ddelta));

581:     PetscCall(PetscObjectReference((PetscObject)BD1));
582:     ctx->BD = BD1;
583:     PetscCall(KSPCreate(PETSC_COMM_SELF, &ctx->kBD));
584:     PetscCall(KSPSetNestLevel(ctx->kBD, fetidpmat_ctx->pc->kspnestlevel));
585:     PetscCall(KSPSetOperators(ctx->kBD, BD2, BD2));
586:     PetscCall(VecDuplicate(fetidpmat_ctx->lambda_local, &ctx->work));
587:     fetidpmat_ctx->deluxe_nonred = PETSC_TRUE;
588:   }
589:   PetscCall(MatDestroy(&BD1));
590:   PetscCall(MatDestroy(&BD2));

592:   /* fetidpmat sizes */
593:   fetidpmat_ctx->n += nPgl;
594:   fetidpmat_ctx->N = fetidpmat_ctx->n_lambda + nPg;

596:   /* Global vector for FETI-DP linear system */
597:   PetscCall(VecCreate(comm, &fetidp_global));
598:   PetscCall(VecSetSizes(fetidp_global, fetidpmat_ctx->n, fetidpmat_ctx->N));
599:   PetscCall(VecSetType(fetidp_global, VECMPI));
600:   PetscCall(VecSetUp(fetidp_global));

602:   /* Decide layout for fetidp dofs: if it is a saddle point problem
603:      pressure is ordered first in the local part of the global vector
604:      of the FETI-DP linear system */
605:   if (nPg) {
606:     Vec             v;
607:     IS              IS_l2g_p, ais;
608:     PetscLayout     alay;
609:     const PetscInt *idxs, *pranges, *aranges, *lranges;
610:     PetscInt       *l2g_indices_p, rst;
611:     PetscMPIInt     rank;

613:     PetscCall(PetscMalloc1(nPl, &l2g_indices_p));
614:     PetscCall(VecGetLayout(fetidp_global, &alay));
615:     PetscCall(PetscLayoutGetRanges(alay, &aranges));
616:     PetscCall(PetscLayoutGetRanges(play, &pranges));
617:     PetscCall(PetscLayoutGetRanges(llay, &lranges));

619:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fetidp_global), &rank));
620:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), pranges[rank + 1] - pranges[rank], aranges[rank], 1, &fetidpmat_ctx->pressure));
621:     PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->pressure, "F_P"));
622:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), lranges[rank + 1] - lranges[rank], aranges[rank] + pranges[rank + 1] - pranges[rank], 1, &fetidpmat_ctx->lagrange));
623:     PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->lagrange, "F_L"));
624:     PetscCall(ISLocalToGlobalMappingGetIndices(l2gmap_p, &idxs));
625:     /* shift local to global indices for pressure */
626:     for (i = 0; i < nPl; i++) {
627:       PetscMPIInt owner;

629:       PetscCall(PetscLayoutFindOwner(play, idxs[i], &owner));
630:       l2g_indices_p[i] = idxs[i] - pranges[owner] + aranges[owner];
631:     }
632:     PetscCall(ISLocalToGlobalMappingRestoreIndices(l2gmap_p, &idxs));
633:     PetscCall(ISCreateGeneral(comm, nPl, l2g_indices_p, PETSC_OWN_POINTER, &IS_l2g_p));

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

639:     /* scatter for lagrange multipliers only */
640:     PetscCall(VecCreate(comm, &v));
641:     PetscCall(VecSetType(v, VECSTANDARD));
642:     PetscCall(VecSetLayout(v, llay));
643:     PetscCall(VecSetUp(v));
644:     PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &ais));
645:     PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, v, ais, &fetidpmat_ctx->l2g_lambda_only));
646:     PetscCall(ISDestroy(&ais));
647:     PetscCall(VecDestroy(&v));

649:     /* shift local to global indices for multipliers */
650:     for (i = 0; i < n_local_lambda; i++) {
651:       PetscInt    ps;
652:       PetscMPIInt owner;

654:       PetscCall(PetscLayoutFindOwner(llay, l2g_indices[i], &owner));
655:       ps             = pranges[owner + 1] - pranges[owner];
656:       l2g_indices[i] = l2g_indices[i] - lranges[owner] + aranges[owner] + ps;
657:     }

659:     /* scatter from alldofs to pressures global fetidp vector */
660:     PetscCall(PetscLayoutGetRange(alay, &rst, NULL));
661:     PetscCall(ISCreateStride(comm, nPgl, rst, 1, &ais));
662:     PetscCall(VecScatterCreate(pcis->vec1_global, pP, fetidp_global, ais, &fetidpmat_ctx->g2g_p));
663:     PetscCall(ISDestroy(&ais));
664:   }
665:   PetscCall(PetscLayoutDestroy(&llay));
666:   PetscCall(PetscLayoutDestroy(&play));
667:   PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_OWN_POINTER, &IS_l2g_lambda));

669:   /* scatter from local to global multipliers */
670:   PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, fetidp_global, IS_l2g_lambda, &fetidpmat_ctx->l2g_lambda));
671:   PetscCall(ISDestroy(&IS_l2g_lambda));
672:   PetscCall(ISLocalToGlobalMappingDestroy(&l2gmap_p));
673:   PetscCall(VecDestroy(&fetidp_global));

675:   /* Create some work vectors needed by fetidp */
676:   PetscCall(VecDuplicate(pcis->vec1_B, &fetidpmat_ctx->temp_solution_B));
677:   PetscCall(VecDuplicate(pcis->vec1_D, &fetidpmat_ctx->temp_solution_D));
678:   PetscFunctionReturn(PETSC_SUCCESS);
679: }

681: PetscErrorCode PCBDDCSetupFETIDPPCContext(Mat fetimat, FETIDPPC_ctx fetidppc_ctx)
682: {
683:   FETIDPMat_ctx mat_ctx;
684:   PC_BDDC      *pcbddc = (PC_BDDC *)fetidppc_ctx->pc->data;
685:   PC_IS        *pcis   = (PC_IS *)fetidppc_ctx->pc->data;
686:   PetscBool     lumped = PETSC_FALSE;

688:   PetscFunctionBegin;
689:   PetscCall(MatShellGetContext(fetimat, &mat_ctx));
690:   /* get references from objects created when setting up feti mat context */
691:   PetscCall(PetscObjectReference((PetscObject)mat_ctx->lambda_local));
692:   fetidppc_ctx->lambda_local = mat_ctx->lambda_local;
693:   PetscCall(PetscObjectReference((PetscObject)mat_ctx->B_Ddelta));
694:   fetidppc_ctx->B_Ddelta = mat_ctx->B_Ddelta;
695:   if (mat_ctx->deluxe_nonred) {
696:     PC            pc, mpc;
697:     BDdelta_DN    ctx;
698:     MatSolverType solver;
699:     const char   *prefix;

701:     PetscCall(MatShellGetContext(mat_ctx->B_Ddelta, &ctx));
702:     PetscCall(KSPSetType(ctx->kBD, KSPPREONLY));
703:     PetscCall(KSPGetPC(ctx->kBD, &mpc));
704:     PetscCall(KSPGetPC(pcbddc->ksp_D, &pc));
705:     PetscCall(PCSetType(mpc, PCLU));
706:     PetscCall(PCFactorGetMatSolverType(pc, (MatSolverType *)&solver));
707:     if (solver) PetscCall(PCFactorSetMatSolverType(mpc, solver));
708:     PetscCall(MatGetOptionsPrefix(fetimat, &prefix));
709:     PetscCall(KSPSetOptionsPrefix(ctx->kBD, prefix));
710:     PetscCall(KSPAppendOptionsPrefix(ctx->kBD, "bddelta_"));
711:     PetscCall(KSPSetFromOptions(ctx->kBD));
712:   }

714:   if (mat_ctx->l2g_lambda_only) {
715:     PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda_only));
716:     fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda_only;
717:   } else {
718:     PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda));
719:     fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda;
720:   }
721:   /* Dirichlet preconditioner */
722:   PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_lumped", &lumped, NULL));
723:   if (!lumped) {
724:     IS        iV;
725:     PetscBool discrete_harmonic = PETSC_FALSE;

727:     PetscCall(PetscObjectQuery((PetscObject)fetidppc_ctx->pc, "__KSPFETIDP_iV", (PetscObject *)&iV));
728:     if (iV) PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_discrete_harmonic", &discrete_harmonic, NULL));
729:     if (discrete_harmonic) {
730:       KSP             sksp;
731:       PC              pc;
732:       PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs;
733:       Mat             A_II, A_IB, A_BI;
734:       IS              iP = NULL;
735:       PetscBool       isshell, reuse = PETSC_FALSE;
736:       KSPType         ksptype;
737:       const char     *prefix;

739:       /*
740:         We constructs a Schur complement for

742:         | A_II A_ID |
743:         | A_DI A_DD |

745:         instead of

747:         | A_II  B^t_II A_ID |
748:         | B_II -C_II   B_ID |
749:         | A_DI  B^t_ID A_DD |

751:       */
752:       if (sub_schurs && sub_schurs->reuse_solver) {
753:         PetscCall(PetscObjectQuery((PetscObject)sub_schurs->A, "__KSPFETIDP_iP", (PetscObject *)&iP));
754:         if (iP) reuse = PETSC_TRUE;
755:       }
756:       if (!reuse) {
757:         IS       aB;
758:         PetscInt nb;
759:         PetscCall(ISGetLocalSize(pcis->is_B_local, &nb));
760:         PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pcis->A_II), nb, 0, 1, &aB));
761:         PetscCall(MatCreateSubMatrix(pcis->A_II, iV, iV, MAT_INITIAL_MATRIX, &A_II));
762:         PetscCall(MatCreateSubMatrix(pcis->A_IB, iV, aB, MAT_INITIAL_MATRIX, &A_IB));
763:         PetscCall(MatCreateSubMatrix(pcis->A_BI, aB, iV, MAT_INITIAL_MATRIX, &A_BI));
764:         PetscCall(ISDestroy(&aB));
765:       } else {
766:         PetscCall(MatCreateSubMatrix(sub_schurs->A, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &A_IB));
767:         PetscCall(MatCreateSubMatrix(sub_schurs->A, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &A_BI));
768:         PetscCall(PetscObjectReference((PetscObject)pcis->A_II));
769:         A_II = pcis->A_II;
770:       }
771:       PetscCall(MatCreateSchurComplement(A_II, A_II, A_IB, A_BI, pcis->A_BB, &fetidppc_ctx->S_j));

773:       /* propagate settings of solver */
774:       PetscCall(MatSchurComplementGetKSP(fetidppc_ctx->S_j, &sksp));
775:       PetscCall(KSPGetType(pcis->ksp_D, &ksptype));
776:       PetscCall(KSPSetType(sksp, ksptype));
777:       PetscCall(KSPGetPC(pcis->ksp_D, &pc));
778:       PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCSHELL, &isshell));
779:       if (!isshell) {
780:         MatSolverType solver;
781:         PCType        pctype;

783:         PetscCall(PCGetType(pc, &pctype));
784:         PetscCall(PCFactorGetMatSolverType(pc, (MatSolverType *)&solver));
785:         PetscCall(KSPGetPC(sksp, &pc));
786:         PetscCall(PCSetType(pc, pctype));
787:         if (solver) PetscCall(PCFactorSetMatSolverType(pc, solver));
788:       } else {
789:         PetscCall(KSPGetPC(sksp, &pc));
790:         PetscCall(PCSetType(pc, PCLU));
791:       }
792:       PetscCall(MatDestroy(&A_II));
793:       PetscCall(MatDestroy(&A_IB));
794:       PetscCall(MatDestroy(&A_BI));
795:       PetscCall(MatGetOptionsPrefix(fetimat, &prefix));
796:       PetscCall(KSPSetOptionsPrefix(sksp, prefix));
797:       PetscCall(KSPAppendOptionsPrefix(sksp, "harmonic_"));
798:       PetscCall(KSPSetFromOptions(sksp));
799:       if (reuse) {
800:         PetscCall(KSPSetPC(sksp, sub_schurs->reuse_solver->interior_solver));
801:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)sub_schurs->reuse_solver->interior_solver, (PetscObject)sksp, 0));
802:       }
803:     } else { /* default Dirichlet preconditioner is pde-harmonic */
804:       PetscCall(MatCreateSchurComplement(pcis->A_II, pcis->A_II, pcis->A_IB, pcis->A_BI, pcis->A_BB, &fetidppc_ctx->S_j));
805:       PetscCall(MatSchurComplementSetKSP(fetidppc_ctx->S_j, pcis->ksp_D));
806:     }
807:   } else {
808:     PetscCall(PetscObjectReference((PetscObject)pcis->A_BB));
809:     fetidppc_ctx->S_j = pcis->A_BB;
810:   }
811:   /* saddle-point */
812:   if (mat_ctx->xPg) {
813:     PetscCall(PetscObjectReference((PetscObject)mat_ctx->xPg));
814:     fetidppc_ctx->xPg = mat_ctx->xPg;
815:     PetscCall(PetscObjectReference((PetscObject)mat_ctx->yPg));
816:     fetidppc_ctx->yPg = mat_ctx->yPg;
817:   }
818:   PetscFunctionReturn(PETSC_SUCCESS);
819: }

821: static PetscErrorCode FETIDPMatMult_Kernel(Mat fetimat, Vec x, Vec y, PetscBool trans)
822: {
823:   FETIDPMat_ctx mat_ctx;
824:   PC_BDDC      *pcbddc;
825:   PC_IS        *pcis;

827:   PetscFunctionBegin;
828:   PetscCall(MatShellGetContext(fetimat, &mat_ctx));
829:   pcis   = (PC_IS *)mat_ctx->pc->data;
830:   pcbddc = (PC_BDDC *)mat_ctx->pc->data;
831:   /* Application of B_delta^T */
832:   PetscCall(VecSet(pcis->vec1_B, 0.));
833:   PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
834:   PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
835:   PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B));

837:   /* Add contribution from saddle point */
838:   if (mat_ctx->l2g_p) {
839:     PetscCall(VecScatterBegin(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
840:     PetscCall(VecScatterEnd(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
841:     if (pcbddc->switch_static) {
842:       if (trans) {
843:         PetscCall(MatMultTranspose(mat_ctx->B_BI, mat_ctx->vP, pcis->vec1_D));
844:       } else {
845:         PetscCall(MatMult(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D));
846:       }
847:     }
848:     if (trans) {
849:       PetscCall(MatMultTransposeAdd(mat_ctx->B_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
850:     } else {
851:       PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
852:     }
853:   } else {
854:     if (pcbddc->switch_static) PetscCall(VecSet(pcis->vec1_D, 0.0));
855:   }
856:   /* Application of \widetilde{S}^-1 */
857:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
858:   PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, trans));
859:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
860:   PetscCall(VecSet(y, 0.0));
861:   /* Application of B_delta */
862:   PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local));
863:   /* Contribution from boundary pressures */
864:   if (mat_ctx->C) {
865:     const PetscScalar *lx;
866:     PetscScalar       *ly;

868:     /* pressure ordered first in the local part of x and y */
869:     PetscCall(VecGetArrayRead(x, &lx));
870:     PetscCall(VecGetArray(y, &ly));
871:     PetscCall(VecPlaceArray(mat_ctx->xPg, lx));
872:     PetscCall(VecPlaceArray(mat_ctx->yPg, ly));
873:     if (trans) {
874:       PetscCall(MatMultTranspose(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg));
875:     } else {
876:       PetscCall(MatMult(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg));
877:     }
878:     PetscCall(VecResetArray(mat_ctx->xPg));
879:     PetscCall(VecResetArray(mat_ctx->yPg));
880:     PetscCall(VecRestoreArrayRead(x, &lx));
881:     PetscCall(VecRestoreArray(y, &ly));
882:   }
883:   /* Add contribution from saddle point */
884:   if (mat_ctx->l2g_p) {
885:     if (trans) {
886:       PetscCall(MatMultTranspose(mat_ctx->Bt_BB, pcis->vec1_B, mat_ctx->vP));
887:     } else {
888:       PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP));
889:     }
890:     if (pcbddc->switch_static) {
891:       if (trans) {
892:         PetscCall(MatMultTransposeAdd(mat_ctx->Bt_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
893:       } else {
894:         PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
895:       }
896:     }
897:     PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD));
898:     PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD));
899:   }
900:   PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD));
901:   PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD));
902:   PetscFunctionReturn(PETSC_SUCCESS);
903: }

905: PetscErrorCode FETIDPMatMult(Mat fetimat, Vec x, Vec y)
906: {
907:   PetscFunctionBegin;
908:   PetscCall(FETIDPMatMult_Kernel(fetimat, x, y, PETSC_FALSE));
909:   PetscFunctionReturn(PETSC_SUCCESS);
910: }

912: PetscErrorCode FETIDPMatMultTranspose(Mat fetimat, Vec x, Vec y)
913: {
914:   PetscFunctionBegin;
915:   PetscCall(FETIDPMatMult_Kernel(fetimat, x, y, PETSC_TRUE));
916:   PetscFunctionReturn(PETSC_SUCCESS);
917: }

919: static PetscErrorCode FETIDPPCApply_Kernel(PC fetipc, Vec x, Vec y, PetscBool trans)
920: {
921:   FETIDPPC_ctx pc_ctx;
922:   PC_IS       *pcis;

924:   PetscFunctionBegin;
925:   PetscCall(PCShellGetContext(fetipc, &pc_ctx));
926:   pcis = (PC_IS *)pc_ctx->pc->data;
927:   /* Application of B_Ddelta^T */
928:   PetscCall(VecScatterBegin(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
929:   PetscCall(VecScatterEnd(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
930:   PetscCall(VecSet(pcis->vec2_B, 0.0));
931:   PetscCall(MatMultTranspose(pc_ctx->B_Ddelta, pc_ctx->lambda_local, pcis->vec2_B));
932:   /* Application of local Schur complement */
933:   if (trans) {
934:     PetscCall(MatMultTranspose(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B));
935:   } else {
936:     PetscCall(MatMult(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B));
937:   }
938:   /* Application of B_Ddelta */
939:   PetscCall(MatMult(pc_ctx->B_Ddelta, pcis->vec1_B, pc_ctx->lambda_local));
940:   PetscCall(VecSet(y, 0.0));
941:   PetscCall(VecScatterBegin(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD));
942:   PetscCall(VecScatterEnd(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD));
943:   PetscFunctionReturn(PETSC_SUCCESS);
944: }

946: PetscErrorCode FETIDPPCApply(PC pc, Vec x, Vec y)
947: {
948:   PetscFunctionBegin;
949:   PetscCall(FETIDPPCApply_Kernel(pc, x, y, PETSC_FALSE));
950:   PetscFunctionReturn(PETSC_SUCCESS);
951: }

953: PetscErrorCode FETIDPPCApplyTranspose(PC pc, Vec x, Vec y)
954: {
955:   PetscFunctionBegin;
956:   PetscCall(FETIDPPCApply_Kernel(pc, x, y, PETSC_TRUE));
957:   PetscFunctionReturn(PETSC_SUCCESS);
958: }

960: PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer)
961: {
962:   FETIDPPC_ctx pc_ctx;
963:   PetscBool    iascii;
964:   PetscViewer  sviewer;

966:   PetscFunctionBegin;
967:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
968:   if (iascii) {
969:     PetscMPIInt rank;
970:     PetscBool   isschur, isshell;

972:     PetscCall(PCShellGetContext(pc, &pc_ctx));
973:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
974:     PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->S_j, MATSCHURCOMPLEMENT, &isschur));
975:     if (isschur) {
976:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Dirichlet preconditioner (just from rank 0)\n"));
977:     } else {
978:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Lumped preconditioner (just from rank 0)\n"));
979:     }
980:     PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer));
981:     if (rank == 0) {
982:       PetscCall(PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO));
983:       PetscCall(PetscViewerASCIIPushTab(sviewer));
984:       PetscCall(MatView(pc_ctx->S_j, sviewer));
985:       PetscCall(PetscViewerASCIIPopTab(sviewer));
986:       PetscCall(PetscViewerPopFormat(sviewer));
987:     }
988:     PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer));
989:     PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->B_Ddelta, MATSHELL, &isshell));
990:     if (isshell) {
991:       BDdelta_DN ctx;
992:       PetscCall(PetscViewerASCIIPrintf(viewer, "  FETI-DP BDdelta: DB^t * (B D^-1 B^t)^-1 for deluxe scaling (just from rank 0)\n"));
993:       PetscCall(MatShellGetContext(pc_ctx->B_Ddelta, &ctx));
994:       PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer));
995:       if (rank == 0) {
996:         PetscInt tl;

998:         PetscCall(PetscViewerASCIIGetTab(sviewer, &tl));
999:         PetscCall(PetscObjectSetTabLevel((PetscObject)ctx->kBD, tl));
1000:         PetscCall(KSPView(ctx->kBD, sviewer));
1001:         PetscCall(PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO));
1002:         PetscCall(MatView(ctx->BD, sviewer));
1003:         PetscCall(PetscViewerPopFormat(sviewer));
1004:       }
1005:       PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer));
1006:     }
1007:     PetscCall(PetscViewerFlush(viewer));
1008:   }
1009:   PetscFunctionReturn(PETSC_SUCCESS);
1010: }