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:   MatShellGetContext(A, &ctx);
 10:   MatMultTranspose(ctx->BD, x, ctx->work);
 11:   KSPSolveTranspose(ctx->kBD, ctx->work, y);
 12:   /* No PC so cannot propagate up failure in KSPSolveTranspose() */
 13:   return 0;
 14: }

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

 20:   MatShellGetContext(A, &ctx);
 21:   KSPSolve(ctx->kBD, x, ctx->work);
 22:   /* No PC so cannot propagate up failure in KSPSolve() */
 23:   MatMult(ctx->BD, ctx->work, y);
 24:   return 0;
 25: }

 27: static PetscErrorCode MatDestroy_BDdelta_deluxe_nonred(Mat A)
 28: {
 29:   BDdelta_DN ctx;

 31:   MatShellGetContext(A, &ctx);
 32:   MatDestroy(&ctx->BD);
 33:   KSPDestroy(&ctx->kBD);
 34:   VecDestroy(&ctx->work);
 35:   PetscFree(ctx);
 36:   return 0;
 37: }

 39: PetscErrorCode PCBDDCCreateFETIDPMatContext(PC pc, FETIDPMat_ctx *fetidpmat_ctx)
 40: {
 41:   FETIDPMat_ctx newctx;

 43:   PetscNew(&newctx);
 44:   /* increase the reference count for BDDC preconditioner */
 45:   PetscObjectReference((PetscObject)pc);
 46:   newctx->pc     = pc;
 47:   *fetidpmat_ctx = newctx;
 48:   return 0;
 49: }

 51: PetscErrorCode PCBDDCCreateFETIDPPCContext(PC pc, FETIDPPC_ctx *fetidppc_ctx)
 52: {
 53:   FETIDPPC_ctx newctx;

 55:   PetscNew(&newctx);
 56:   /* increase the reference count for BDDC preconditioner */
 57:   PetscObjectReference((PetscObject)pc);
 58:   newctx->pc    = pc;
 59:   *fetidppc_ctx = newctx;
 60:   return 0;
 61: }

 63: PetscErrorCode PCBDDCDestroyFETIDPMat(Mat A)
 64: {
 65:   FETIDPMat_ctx mat_ctx;

 67:   MatShellGetContext(A, &mat_ctx);
 68:   VecDestroy(&mat_ctx->lambda_local);
 69:   VecDestroy(&mat_ctx->temp_solution_D);
 70:   VecDestroy(&mat_ctx->temp_solution_B);
 71:   MatDestroy(&mat_ctx->B_delta);
 72:   MatDestroy(&mat_ctx->B_Ddelta);
 73:   MatDestroy(&mat_ctx->B_BB);
 74:   MatDestroy(&mat_ctx->B_BI);
 75:   MatDestroy(&mat_ctx->Bt_BB);
 76:   MatDestroy(&mat_ctx->Bt_BI);
 77:   MatDestroy(&mat_ctx->C);
 78:   VecDestroy(&mat_ctx->rhs_flip);
 79:   VecDestroy(&mat_ctx->vP);
 80:   VecDestroy(&mat_ctx->xPg);
 81:   VecDestroy(&mat_ctx->yPg);
 82:   VecScatterDestroy(&mat_ctx->l2g_lambda);
 83:   VecScatterDestroy(&mat_ctx->l2g_lambda_only);
 84:   VecScatterDestroy(&mat_ctx->l2g_p);
 85:   VecScatterDestroy(&mat_ctx->g2g_p);
 86:   PCDestroy(&mat_ctx->pc); /* decrease PCBDDC reference count */
 87:   ISDestroy(&mat_ctx->pressure);
 88:   ISDestroy(&mat_ctx->lagrange);
 89:   PetscFree(mat_ctx);
 90:   return 0;
 91: }

 93: PetscErrorCode PCBDDCDestroyFETIDPPC(PC pc)
 94: {
 95:   FETIDPPC_ctx pc_ctx;

 97:   PCShellGetContext(pc, &pc_ctx);
 98:   VecDestroy(&pc_ctx->lambda_local);
 99:   MatDestroy(&pc_ctx->B_Ddelta);
100:   VecScatterDestroy(&pc_ctx->l2g_lambda);
101:   MatDestroy(&pc_ctx->S_j);
102:   PCDestroy(&pc_ctx->pc); /* decrease PCBDDC reference count */
103:   VecDestroy(&pc_ctx->xPg);
104:   VecDestroy(&pc_ctx->yPg);
105:   PetscFree(pc_ctx);
106:   return 0;
107: }

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

134:   /* saddlepoint */
135:   ISLocalToGlobalMapping l2gmap_p;
136:   PetscLayout            play;
137:   IS                     gP, pP;
138:   PetscInt               nPl, nPg, nPgl;

140:   PetscObjectGetComm((PetscObject)(fetidpmat_ctx->pc), &comm);
141:   MPI_Comm_rank(comm, &rank);
142:   MPI_Comm_size(comm, &size);

144:   /* saddlepoint */
145:   nPl      = 0;
146:   nPg      = 0;
147:   nPgl     = 0;
148:   gP       = NULL;
149:   pP       = NULL;
150:   l2gmap_p = NULL;
151:   play     = NULL;
152:   PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_pP", (PetscObject *)&pP);
153:   if (pP) { /* saddle point */
154:     /* subdomain pressures in global numbering */
155:     PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_gP", (PetscObject *)&gP);
157:     ISGetLocalSize(gP, &nPl);
158:     VecCreate(PETSC_COMM_SELF, &fetidpmat_ctx->vP);
159:     VecSetSizes(fetidpmat_ctx->vP, nPl, nPl);
160:     VecSetType(fetidpmat_ctx->vP, VECSTANDARD);
161:     VecSetUp(fetidpmat_ctx->vP);

163:     /* pressure matrix */
164:     PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_C", (PetscObject *)&fetidpmat_ctx->C);
165:     if (!fetidpmat_ctx->C) { /* null pressure block, compute layout and global numbering for pressures */
166:       IS Pg;

168:       ISRenumber(gP, NULL, &nPg, &Pg);
169:       ISLocalToGlobalMappingCreateIS(Pg, &l2gmap_p);
170:       ISDestroy(&Pg);
171:       PetscLayoutCreate(comm, &play);
172:       PetscLayoutSetBlockSize(play, 1);
173:       PetscLayoutSetSize(play, nPg);
174:       ISGetLocalSize(pP, &nPgl);
175:       PetscLayoutSetLocalSize(play, nPgl);
176:       PetscLayoutSetUp(play);
177:     } else {
178:       PetscObjectReference((PetscObject)fetidpmat_ctx->C);
179:       MatISGetLocalToGlobalMapping(fetidpmat_ctx->C, &l2gmap_p, NULL);
180:       PetscObjectReference((PetscObject)l2gmap_p);
181:       MatGetSize(fetidpmat_ctx->C, &nPg, NULL);
182:       MatGetLocalSize(fetidpmat_ctx->C, NULL, &nPgl);
183:       MatGetLayouts(fetidpmat_ctx->C, NULL, &llay);
184:       PetscLayoutReference(llay, &play);
185:     }
186:     VecCreateMPIWithArray(comm, 1, nPgl, nPg, NULL, &fetidpmat_ctx->xPg);
187:     VecCreateMPIWithArray(comm, 1, nPgl, nPg, NULL, &fetidpmat_ctx->yPg);

189:     /* import matrices for pressures coupling */
190:     PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_B_BI", (PetscObject *)&fetidpmat_ctx->B_BI);
192:     PetscObjectReference((PetscObject)fetidpmat_ctx->B_BI);

194:     PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_B_BB", (PetscObject *)&fetidpmat_ctx->B_BB);
196:     PetscObjectReference((PetscObject)fetidpmat_ctx->B_BB);

198:     PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_Bt_BI", (PetscObject *)&fetidpmat_ctx->Bt_BI);
200:     PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BI);

202:     PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_Bt_BB", (PetscObject *)&fetidpmat_ctx->Bt_BB);
204:     PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BB);

206:     PetscObjectQuery((PetscObject)fetidpmat_ctx->pc, "__KSPFETIDP_flip", (PetscObject *)&fetidpmat_ctx->rhs_flip);
207:     if (fetidpmat_ctx->rhs_flip) PetscObjectReference((PetscObject)fetidpmat_ctx->rhs_flip);
208:   }

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

213:   /* Evaluate local and global number of lagrange multipliers */
214:   VecSet(pcis->vec1_N, 0.0);
215:   n_local_lambda  = 0;
216:   partial_sum     = 0;
217:   n_boundary_dofs = 0;
218:   s               = 0;

220:   /* Get Vertices used to define the BDDC */
221:   PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &isvert);
222:   ISGetLocalSize(isvert, &n_vertices);
223:   ISGetIndices(isvert, &vertex_indices);

225:   dual_size = pcis->n_B - n_vertices;
226:   PetscMalloc1(dual_size, &dual_dofs_boundary_indices);
227:   PetscMalloc1(dual_size, &aux_local_numbering_1);
228:   PetscMalloc1(dual_size, &aux_local_numbering_2);

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

263:   /* compute global ordering of lagrange multipliers and associate l2g map */
264:   ISCreateGeneral(comm, partial_sum, aux_local_numbering_1, PETSC_COPY_VALUES, &subset_n);
265:   ISLocalToGlobalMappingApplyIS(pcis->mapping, subset_n, &subset);
266:   ISDestroy(&subset_n);
267:   ISCreateGeneral(comm, partial_sum, aux_local_numbering_2, PETSC_OWN_POINTER, &subset_mult);
268:   ISRenumber(subset, subset_mult, &fetidpmat_ctx->n_lambda, &subset_n);
269:   ISDestroy(&subset);

271:   if (PetscDefined(USE_DEBUG)) {
272:     VecSet(pcis->vec1_global, 0.0);
273:     VecScatterBegin(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE);
274:     VecScatterEnd(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE);
275:     VecSum(pcis->vec1_global, &scalar_value);
276:     i = (PetscInt)PetscRealPart(scalar_value);
278:   }

280:   /* init data for scaling factors exchange */
281:   if (!pcbddc->use_deluxe_scaling) {
282:     PetscInt    *ptrs_buffer, neigh_position;
283:     PetscScalar *send_buffer, *recv_buffer;
284:     MPI_Request *send_reqs, *recv_reqs;

286:     partial_sum = 0;
287:     PetscMalloc1(pcis->n_neigh, &ptrs_buffer);
288:     PetscMalloc1(PetscMax(pcis->n_neigh - 1, 0), &send_reqs);
289:     PetscMalloc1(PetscMax(pcis->n_neigh - 1, 0), &recv_reqs);
290:     PetscMalloc1(pcis->n + 1, &all_factors);
291:     if (pcis->n_neigh > 0) ptrs_buffer[0] = 0;
292:     for (i = 1; i < pcis->n_neigh; i++) {
293:       partial_sum += pcis->n_shared[i];
294:       ptrs_buffer[i] = ptrs_buffer[i - 1] + pcis->n_shared[i];
295:     }
296:     PetscMalloc1(partial_sum, &send_buffer);
297:     PetscMalloc1(partial_sum, &recv_buffer);
298:     PetscMalloc1(partial_sum, &all_factors[0]);
299:     for (i = 0; i < pcis->n - 1; i++) {
300:       j                  = mat_graph->count[i];
301:       all_factors[i + 1] = all_factors[i] + j;
302:     }

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

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

414:   /* Create local part of B_delta */
415:   MatCreate(PETSC_COMM_SELF, &fetidpmat_ctx->B_delta);
416:   MatSetSizes(fetidpmat_ctx->B_delta, n_local_lambda, pcis->n_B, n_local_lambda, pcis->n_B);
417:   MatSetType(fetidpmat_ctx->B_delta, MATSEQAIJ);
418:   MatSeqAIJSetPreallocation(fetidpmat_ctx->B_delta, 1, NULL);
419:   MatSetOption(fetidpmat_ctx->B_delta, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
420:   for (i = 0; i < n_local_lambda; i++) MatSetValue(fetidpmat_ctx->B_delta, i, cols_B_delta[i], vals_B_delta[i], INSERT_VALUES);
421:   PetscFree(vals_B_delta);
422:   MatAssemblyBegin(fetidpmat_ctx->B_delta, MAT_FINAL_ASSEMBLY);
423:   MatAssemblyEnd(fetidpmat_ctx->B_delta, MAT_FINAL_ASSEMBLY);

425:   BD1 = NULL;
426:   BD2 = NULL;
427:   if (fully_redundant) {
429:     MatCreate(PETSC_COMM_SELF, &ScalingMat);
430:     MatSetSizes(ScalingMat, n_local_lambda, n_local_lambda, n_local_lambda, n_local_lambda);
431:     MatSetType(ScalingMat, MATSEQAIJ);
432:     MatSeqAIJSetPreallocation(ScalingMat, 1, NULL);
433:     for (i = 0; i < n_local_lambda; i++) MatSetValue(ScalingMat, i, i, scaling_factors[i], INSERT_VALUES);
434:     MatAssemblyBegin(ScalingMat, MAT_FINAL_ASSEMBLY);
435:     MatAssemblyEnd(ScalingMat, MAT_FINAL_ASSEMBLY);
436:     MatMatMult(ScalingMat, fetidpmat_ctx->B_delta, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &fetidpmat_ctx->B_Ddelta);
437:     MatDestroy(&ScalingMat);
438:   } else {
439:     MatCreate(PETSC_COMM_SELF, &fetidpmat_ctx->B_Ddelta);
440:     MatSetSizes(fetidpmat_ctx->B_Ddelta, n_local_lambda, pcis->n_B, n_local_lambda, pcis->n_B);
441:     if (!pcbddc->use_deluxe_scaling || !pcbddc->sub_schurs) {
442:       MatSetType(fetidpmat_ctx->B_Ddelta, MATSEQAIJ);
443:       MatSeqAIJSetPreallocation(fetidpmat_ctx->B_Ddelta, 1, NULL);
444:       for (i = 0; i < n_local_lambda; i++) MatSetValue(fetidpmat_ctx->B_Ddelta, i, cols_B_delta[i], scaling_factors[i], INSERT_VALUES);
445:       MatAssemblyBegin(fetidpmat_ctx->B_Ddelta, MAT_FINAL_ASSEMBLY);
446:       MatAssemblyEnd(fetidpmat_ctx->B_Ddelta, MAT_FINAL_ASSEMBLY);
447:     } else {
448:       /* scaling as in Klawonn-Widlund 1999 */
449:       PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
450:       PCBDDCSubSchurs     sub_schurs = pcbddc->sub_schurs;
451:       Mat                 T;
452:       PetscScalar        *W, lwork, *Bwork;
453:       const PetscInt     *idxs = NULL;
454:       PetscInt            cum, mss, *nnz;
455:       PetscBLASInt       *pivots, B_lwork, B_N, B_ierr;

458:       mss = 0;
459:       PetscCalloc1(pcis->n_B, &nnz);
460:       if (sub_schurs->is_Ej_all) {
461:         ISGetIndices(sub_schurs->is_Ej_all, &idxs);
462:         for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
463:           PetscInt subset_size;

465:           ISGetLocalSize(sub_schurs->is_subs[i], &subset_size);
466:           for (j = 0; j < subset_size; j++) nnz[idxs[j + cum]] = subset_size;
467:           mss = PetscMax(mss, subset_size);
468:           cum += subset_size;
469:         }
470:       }
471:       MatCreate(PETSC_COMM_SELF, &T);
472:       MatSetSizes(T, pcis->n_B, pcis->n_B, pcis->n_B, pcis->n_B);
473:       MatSetType(T, MATSEQAIJ);
474:       MatSeqAIJSetPreallocation(T, 0, nnz);
475:       PetscFree(nnz);

477:       /* workspace allocation */
478:       B_lwork = 0;
479:       if (mss) {
480:         PetscScalar dummy = 1;

482:         B_lwork = -1;
483:         PetscBLASIntCast(mss, &B_N);
484:         PetscFPTrapPush(PETSC_FP_TRAP_OFF);
485:         PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&B_N, &dummy, &B_N, &B_N, &lwork, &B_lwork, &B_ierr));
486:         PetscFPTrapPop();
488:         PetscBLASIntCast((PetscInt)PetscRealPart(lwork), &B_lwork);
489:       }
490:       PetscMalloc3(mss * mss, &W, mss, &pivots, B_lwork, &Bwork);

492:       for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
493:         const PetscScalar *M;
494:         PetscInt           subset_size;

496:         ISGetLocalSize(sub_schurs->is_subs[i], &subset_size);
497:         PetscBLASIntCast(subset_size, &B_N);
498:         MatDenseGetArrayRead(deluxe_ctx->seq_mat[i], &M);
499:         PetscArraycpy(W, M, subset_size * subset_size);
500:         MatDenseRestoreArrayRead(deluxe_ctx->seq_mat[i], &M);
501:         PetscFPTrapPush(PETSC_FP_TRAP_OFF);
502:         PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&B_N, &B_N, W, &B_N, pivots, &B_ierr));
504:         PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&B_N, W, &B_N, pivots, Bwork, &B_lwork, &B_ierr));
506:         PetscFPTrapPop();
507:         /* silent static analyzer */
509:         MatSetValues(T, subset_size, idxs + cum, subset_size, idxs + cum, W, INSERT_VALUES);
510:         cum += subset_size;
511:       }
512:       MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY);
513:       MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY);
514:       MatMatTransposeMult(T, fetidpmat_ctx->B_delta, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &BD1);
515:       MatMatMult(fetidpmat_ctx->B_delta, BD1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &BD2);
516:       MatDestroy(&T);
517:       PetscFree3(W, pivots, Bwork);
518:       if (sub_schurs->is_Ej_all) ISRestoreIndices(sub_schurs->is_Ej_all, &idxs);
519:     }
520:   }
521:   PetscFree(scaling_factors);
522:   PetscFree(cols_B_delta);

524:   /* Layout of multipliers */
525:   PetscLayoutCreate(comm, &llay);
526:   PetscLayoutSetBlockSize(llay, 1);
527:   PetscLayoutSetSize(llay, fetidpmat_ctx->n_lambda);
528:   PetscLayoutSetUp(llay);
529:   PetscLayoutGetLocalSize(llay, &fetidpmat_ctx->n);

531:   /* Local work vector of multipliers */
532:   VecCreate(PETSC_COMM_SELF, &fetidpmat_ctx->lambda_local);
533:   VecSetSizes(fetidpmat_ctx->lambda_local, n_local_lambda, n_local_lambda);
534:   VecSetType(fetidpmat_ctx->lambda_local, VECSEQ);

536:   if (BD2) {
537:     ISLocalToGlobalMapping l2g;
538:     Mat                    T, TA, *pT;
539:     IS                     is;
540:     PetscInt               nl, N;
541:     BDdelta_DN             ctx;

543:     PetscLayoutGetLocalSize(llay, &nl);
544:     PetscLayoutGetSize(llay, &N);
545:     MatCreate(comm, &T);
546:     MatSetSizes(T, nl, nl, N, N);
547:     MatSetType(T, MATIS);
548:     ISLocalToGlobalMappingCreate(comm, 1, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &l2g);
549:     MatSetLocalToGlobalMapping(T, l2g, l2g);
550:     ISLocalToGlobalMappingDestroy(&l2g);
551:     MatISSetLocalMat(T, BD2);
552:     MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY);
553:     MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY);
554:     MatDestroy(&BD2);
555:     MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &TA);
556:     MatDestroy(&T);
557:     ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_USE_POINTER, &is);
558:     MatCreateSubMatrices(TA, 1, &is, &is, MAT_INITIAL_MATRIX, &pT);
559:     MatDestroy(&TA);
560:     ISDestroy(&is);
561:     BD2 = pT[0];
562:     PetscFree(pT);

564:     /* B_Ddelta for non-redundant multipliers with deluxe scaling */
565:     PetscNew(&ctx);
566:     MatSetType(fetidpmat_ctx->B_Ddelta, MATSHELL);
567:     MatShellSetContext(fetidpmat_ctx->B_Ddelta, ctx);
568:     MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT, (void (*)(void))MatMult_BDdelta_deluxe_nonred);
569:     MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_BDdelta_deluxe_nonred);
570:     MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_DESTROY, (void (*)(void))MatDestroy_BDdelta_deluxe_nonred);
571:     MatSetUp(fetidpmat_ctx->B_Ddelta);

573:     PetscObjectReference((PetscObject)BD1);
574:     ctx->BD = BD1;
575:     KSPCreate(PETSC_COMM_SELF, &ctx->kBD);
576:     KSPSetOperators(ctx->kBD, BD2, BD2);
577:     VecDuplicate(fetidpmat_ctx->lambda_local, &ctx->work);
578:     fetidpmat_ctx->deluxe_nonred = PETSC_TRUE;
579:   }
580:   MatDestroy(&BD1);
581:   MatDestroy(&BD2);

583:   /* fetidpmat sizes */
584:   fetidpmat_ctx->n += nPgl;
585:   fetidpmat_ctx->N = fetidpmat_ctx->n_lambda + nPg;

587:   /* Global vector for FETI-DP linear system */
588:   VecCreate(comm, &fetidp_global);
589:   VecSetSizes(fetidp_global, fetidpmat_ctx->n, fetidpmat_ctx->N);
590:   VecSetType(fetidp_global, VECMPI);
591:   VecSetUp(fetidp_global);

593:   /* Decide layout for fetidp dofs: if it is a saddle point problem
594:      pressure is ordered first in the local part of the global vector
595:      of the FETI-DP linear system */
596:   if (nPg) {
597:     Vec             v;
598:     IS              IS_l2g_p, ais;
599:     PetscLayout     alay;
600:     const PetscInt *idxs, *pranges, *aranges, *lranges;
601:     PetscInt       *l2g_indices_p, rst;
602:     PetscMPIInt     rank;

604:     PetscMalloc1(nPl, &l2g_indices_p);
605:     VecGetLayout(fetidp_global, &alay);
606:     PetscLayoutGetRanges(alay, &aranges);
607:     PetscLayoutGetRanges(play, &pranges);
608:     PetscLayoutGetRanges(llay, &lranges);

610:     MPI_Comm_rank(PetscObjectComm((PetscObject)fetidp_global), &rank);
611:     ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), pranges[rank + 1] - pranges[rank], aranges[rank], 1, &fetidpmat_ctx->pressure);
612:     PetscObjectSetName((PetscObject)fetidpmat_ctx->pressure, "F_P");
613:     ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), lranges[rank + 1] - lranges[rank], aranges[rank] + pranges[rank + 1] - pranges[rank], 1, &fetidpmat_ctx->lagrange);
614:     PetscObjectSetName((PetscObject)fetidpmat_ctx->lagrange, "F_L");
615:     ISLocalToGlobalMappingGetIndices(l2gmap_p, &idxs);
616:     /* shift local to global indices for pressure */
617:     for (i = 0; i < nPl; i++) {
618:       PetscMPIInt owner;

620:       PetscLayoutFindOwner(play, idxs[i], &owner);
621:       l2g_indices_p[i] = idxs[i] - pranges[owner] + aranges[owner];
622:     }
623:     ISLocalToGlobalMappingRestoreIndices(l2gmap_p, &idxs);
624:     ISCreateGeneral(comm, nPl, l2g_indices_p, PETSC_OWN_POINTER, &IS_l2g_p);

626:     /* local to global scatter for pressure */
627:     VecScatterCreate(fetidpmat_ctx->vP, NULL, fetidp_global, IS_l2g_p, &fetidpmat_ctx->l2g_p);
628:     ISDestroy(&IS_l2g_p);

630:     /* scatter for lagrange multipliers only */
631:     VecCreate(comm, &v);
632:     VecSetType(v, VECSTANDARD);
633:     VecSetLayout(v, llay);
634:     VecSetUp(v);
635:     ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &ais);
636:     VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, v, ais, &fetidpmat_ctx->l2g_lambda_only);
637:     ISDestroy(&ais);
638:     VecDestroy(&v);

640:     /* shift local to global indices for multipliers */
641:     for (i = 0; i < n_local_lambda; i++) {
642:       PetscInt    ps;
643:       PetscMPIInt owner;

645:       PetscLayoutFindOwner(llay, l2g_indices[i], &owner);
646:       ps             = pranges[owner + 1] - pranges[owner];
647:       l2g_indices[i] = l2g_indices[i] - lranges[owner] + aranges[owner] + ps;
648:     }

650:     /* scatter from alldofs to pressures global fetidp vector */
651:     PetscLayoutGetRange(alay, &rst, NULL);
652:     ISCreateStride(comm, nPgl, rst, 1, &ais);
653:     VecScatterCreate(pcis->vec1_global, pP, fetidp_global, ais, &fetidpmat_ctx->g2g_p);
654:     ISDestroy(&ais);
655:   }
656:   PetscLayoutDestroy(&llay);
657:   PetscLayoutDestroy(&play);
658:   ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_OWN_POINTER, &IS_l2g_lambda);

660:   /* scatter from local to global multipliers */
661:   VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, fetidp_global, IS_l2g_lambda, &fetidpmat_ctx->l2g_lambda);
662:   ISDestroy(&IS_l2g_lambda);
663:   ISLocalToGlobalMappingDestroy(&l2gmap_p);
664:   VecDestroy(&fetidp_global);

666:   /* Create some work vectors needed by fetidp */
667:   VecDuplicate(pcis->vec1_B, &fetidpmat_ctx->temp_solution_B);
668:   VecDuplicate(pcis->vec1_D, &fetidpmat_ctx->temp_solution_D);
669:   return 0;
670: }

672: PetscErrorCode PCBDDCSetupFETIDPPCContext(Mat fetimat, FETIDPPC_ctx fetidppc_ctx)
673: {
674:   FETIDPMat_ctx mat_ctx;
675:   PC_BDDC      *pcbddc = (PC_BDDC *)fetidppc_ctx->pc->data;
676:   PC_IS        *pcis   = (PC_IS *)fetidppc_ctx->pc->data;
677:   PetscBool     lumped = PETSC_FALSE;

679:   MatShellGetContext(fetimat, &mat_ctx);
680:   /* get references from objects created when setting up feti mat context */
681:   PetscObjectReference((PetscObject)mat_ctx->lambda_local);
682:   fetidppc_ctx->lambda_local = mat_ctx->lambda_local;
683:   PetscObjectReference((PetscObject)mat_ctx->B_Ddelta);
684:   fetidppc_ctx->B_Ddelta = mat_ctx->B_Ddelta;
685:   if (mat_ctx->deluxe_nonred) {
686:     PC            pc, mpc;
687:     BDdelta_DN    ctx;
688:     MatSolverType solver;
689:     const char   *prefix;

691:     MatShellGetContext(mat_ctx->B_Ddelta, &ctx);
692:     KSPSetType(ctx->kBD, KSPPREONLY);
693:     KSPGetPC(ctx->kBD, &mpc);
694:     KSPGetPC(pcbddc->ksp_D, &pc);
695:     PCSetType(mpc, PCLU);
696:     PCFactorGetMatSolverType(pc, (MatSolverType *)&solver);
697:     if (solver) PCFactorSetMatSolverType(mpc, solver);
698:     MatGetOptionsPrefix(fetimat, &prefix);
699:     KSPSetOptionsPrefix(ctx->kBD, prefix);
700:     KSPAppendOptionsPrefix(ctx->kBD, "bddelta_");
701:     KSPSetFromOptions(ctx->kBD);
702:   }

704:   if (mat_ctx->l2g_lambda_only) {
705:     PetscObjectReference((PetscObject)mat_ctx->l2g_lambda_only);
706:     fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda_only;
707:   } else {
708:     PetscObjectReference((PetscObject)mat_ctx->l2g_lambda);
709:     fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda;
710:   }
711:   /* Dirichlet preconditioner */
712:   PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_lumped", &lumped, NULL);
713:   if (!lumped) {
714:     IS        iV;
715:     PetscBool discrete_harmonic = PETSC_FALSE;

717:     PetscObjectQuery((PetscObject)fetidppc_ctx->pc, "__KSPFETIDP_iV", (PetscObject *)&iV);
718:     if (iV) PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_discrete_harmonic", &discrete_harmonic, NULL);
719:     if (discrete_harmonic) {
720:       KSP             sksp;
721:       PC              pc;
722:       PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs;
723:       Mat             A_II, A_IB, A_BI;
724:       IS              iP = NULL;
725:       PetscBool       isshell, reuse = PETSC_FALSE;
726:       KSPType         ksptype;
727:       const char     *prefix;

729:       /*
730:         We constructs a Schur complement for

732:         | A_II A_ID |
733:         | A_DI A_DD |

735:         instead of

737:         | A_II  B^t_II A_ID |
738:         | B_II -C_II   B_ID |
739:         | A_DI  B^t_ID A_DD |

741:       */
742:       if (sub_schurs && sub_schurs->reuse_solver) {
743:         PetscObjectQuery((PetscObject)sub_schurs->A, "__KSPFETIDP_iP", (PetscObject *)&iP);
744:         if (iP) reuse = PETSC_TRUE;
745:       }
746:       if (!reuse) {
747:         IS       aB;
748:         PetscInt nb;
749:         ISGetLocalSize(pcis->is_B_local, &nb);
750:         ISCreateStride(PetscObjectComm((PetscObject)pcis->A_II), nb, 0, 1, &aB);
751:         MatCreateSubMatrix(pcis->A_II, iV, iV, MAT_INITIAL_MATRIX, &A_II);
752:         MatCreateSubMatrix(pcis->A_IB, iV, aB, MAT_INITIAL_MATRIX, &A_IB);
753:         MatCreateSubMatrix(pcis->A_BI, aB, iV, MAT_INITIAL_MATRIX, &A_BI);
754:         ISDestroy(&aB);
755:       } else {
756:         MatCreateSubMatrix(sub_schurs->A, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &A_IB);
757:         MatCreateSubMatrix(sub_schurs->A, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &A_BI);
758:         PetscObjectReference((PetscObject)pcis->A_II);
759:         A_II = pcis->A_II;
760:       }
761:       MatCreateSchurComplement(A_II, A_II, A_IB, A_BI, pcis->A_BB, &fetidppc_ctx->S_j);

763:       /* propagate settings of solver */
764:       MatSchurComplementGetKSP(fetidppc_ctx->S_j, &sksp);
765:       KSPGetType(pcis->ksp_D, &ksptype);
766:       KSPSetType(sksp, ksptype);
767:       KSPGetPC(pcis->ksp_D, &pc);
768:       PetscObjectTypeCompare((PetscObject)pc, PCSHELL, &isshell);
769:       if (!isshell) {
770:         MatSolverType solver;
771:         PCType        pctype;

773:         PCGetType(pc, &pctype);
774:         PCFactorGetMatSolverType(pc, (MatSolverType *)&solver);
775:         KSPGetPC(sksp, &pc);
776:         PCSetType(pc, pctype);
777:         if (solver) PCFactorSetMatSolverType(pc, solver);
778:       } else {
779:         KSPGetPC(sksp, &pc);
780:         PCSetType(pc, PCLU);
781:       }
782:       MatDestroy(&A_II);
783:       MatDestroy(&A_IB);
784:       MatDestroy(&A_BI);
785:       MatGetOptionsPrefix(fetimat, &prefix);
786:       KSPSetOptionsPrefix(sksp, prefix);
787:       KSPAppendOptionsPrefix(sksp, "harmonic_");
788:       KSPSetFromOptions(sksp);
789:       if (reuse) {
790:         KSPSetPC(sksp, sub_schurs->reuse_solver->interior_solver);
791:         PetscObjectIncrementTabLevel((PetscObject)sub_schurs->reuse_solver->interior_solver, (PetscObject)sksp, 0);
792:       }
793:     } else { /* default Dirichlet preconditioner is pde-harmonic */
794:       MatCreateSchurComplement(pcis->A_II, pcis->A_II, pcis->A_IB, pcis->A_BI, pcis->A_BB, &fetidppc_ctx->S_j);
795:       MatSchurComplementSetKSP(fetidppc_ctx->S_j, pcis->ksp_D);
796:     }
797:   } else {
798:     PetscObjectReference((PetscObject)pcis->A_BB);
799:     fetidppc_ctx->S_j = pcis->A_BB;
800:   }
801:   /* saddle-point */
802:   if (mat_ctx->xPg) {
803:     PetscObjectReference((PetscObject)mat_ctx->xPg);
804:     fetidppc_ctx->xPg = mat_ctx->xPg;
805:     PetscObjectReference((PetscObject)mat_ctx->yPg);
806:     fetidppc_ctx->yPg = mat_ctx->yPg;
807:   }
808:   return 0;
809: }

811: PetscErrorCode FETIDPMatMult_Kernel(Mat fetimat, Vec x, Vec y, PetscBool trans)
812: {
813:   FETIDPMat_ctx mat_ctx;
814:   PC_BDDC      *pcbddc;
815:   PC_IS        *pcis;

817:   MatShellGetContext(fetimat, &mat_ctx);
818:   pcis   = (PC_IS *)mat_ctx->pc->data;
819:   pcbddc = (PC_BDDC *)mat_ctx->pc->data;
820:   /* Application of B_delta^T */
821:   VecSet(pcis->vec1_B, 0.);
822:   VecScatterBegin(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE);
823:   VecScatterEnd(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE);
824:   MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B);

826:   /* Add contribution from saddle point */
827:   if (mat_ctx->l2g_p) {
828:     VecScatterBegin(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE);
829:     VecScatterEnd(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE);
830:     if (pcbddc->switch_static) {
831:       if (trans) {
832:         MatMultTranspose(mat_ctx->B_BI, mat_ctx->vP, pcis->vec1_D);
833:       } else {
834:         MatMult(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D);
835:       }
836:     }
837:     if (trans) {
838:       MatMultTransposeAdd(mat_ctx->B_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B);
839:     } else {
840:       MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B);
841:     }
842:   } else {
843:     if (pcbddc->switch_static) VecSet(pcis->vec1_D, 0.0);
844:   }
845:   /* Application of \widetilde{S}^-1 */
846:   PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n);
847:   PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, trans);
848:   PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n);
849:   VecSet(y, 0.0);
850:   /* Application of B_delta */
851:   MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local);
852:   /* Contribution from boundary pressures */
853:   if (mat_ctx->C) {
854:     const PetscScalar *lx;
855:     PetscScalar       *ly;

857:     /* pressure ordered first in the local part of x and y */
858:     VecGetArrayRead(x, &lx);
859:     VecGetArray(y, &ly);
860:     VecPlaceArray(mat_ctx->xPg, lx);
861:     VecPlaceArray(mat_ctx->yPg, ly);
862:     if (trans) {
863:       MatMultTranspose(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg);
864:     } else {
865:       MatMult(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg);
866:     }
867:     VecResetArray(mat_ctx->xPg);
868:     VecResetArray(mat_ctx->yPg);
869:     VecRestoreArrayRead(x, &lx);
870:     VecRestoreArray(y, &ly);
871:   }
872:   /* Add contribution from saddle point */
873:   if (mat_ctx->l2g_p) {
874:     if (trans) {
875:       MatMultTranspose(mat_ctx->Bt_BB, pcis->vec1_B, mat_ctx->vP);
876:     } else {
877:       MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP);
878:     }
879:     if (pcbddc->switch_static) {
880:       if (trans) {
881:         MatMultTransposeAdd(mat_ctx->Bt_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP);
882:       } else {
883:         MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP);
884:       }
885:     }
886:     VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD);
887:     VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD);
888:   }
889:   VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD);
890:   VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD);
891:   return 0;
892: }

894: PetscErrorCode FETIDPMatMult(Mat fetimat, Vec x, Vec y)
895: {
896:   FETIDPMatMult_Kernel(fetimat, x, y, PETSC_FALSE);
897:   return 0;
898: }

900: PetscErrorCode FETIDPMatMultTranspose(Mat fetimat, Vec x, Vec y)
901: {
902:   FETIDPMatMult_Kernel(fetimat, x, y, PETSC_TRUE);
903:   return 0;
904: }

906: PetscErrorCode FETIDPPCApply_Kernel(PC fetipc, Vec x, Vec y, PetscBool trans)
907: {
908:   FETIDPPC_ctx pc_ctx;
909:   PC_IS       *pcis;

911:   PCShellGetContext(fetipc, &pc_ctx);
912:   pcis = (PC_IS *)pc_ctx->pc->data;
913:   /* Application of B_Ddelta^T */
914:   VecScatterBegin(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE);
915:   VecScatterEnd(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE);
916:   VecSet(pcis->vec2_B, 0.0);
917:   MatMultTranspose(pc_ctx->B_Ddelta, pc_ctx->lambda_local, pcis->vec2_B);
918:   /* Application of local Schur complement */
919:   if (trans) {
920:     MatMultTranspose(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B);
921:   } else {
922:     MatMult(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B);
923:   }
924:   /* Application of B_Ddelta */
925:   MatMult(pc_ctx->B_Ddelta, pcis->vec1_B, pc_ctx->lambda_local);
926:   VecSet(y, 0.0);
927:   VecScatterBegin(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD);
928:   VecScatterEnd(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD);
929:   return 0;
930: }

932: PetscErrorCode FETIDPPCApply(PC pc, Vec x, Vec y)
933: {
934:   FETIDPPCApply_Kernel(pc, x, y, PETSC_FALSE);
935:   return 0;
936: }

938: PetscErrorCode FETIDPPCApplyTranspose(PC pc, Vec x, Vec y)
939: {
940:   FETIDPPCApply_Kernel(pc, x, y, PETSC_TRUE);
941:   return 0;
942: }

944: PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer)
945: {
946:   FETIDPPC_ctx pc_ctx;
947:   PetscBool    iascii;
948:   PetscViewer  sviewer;

950:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
951:   if (iascii) {
952:     PetscMPIInt rank;
953:     PetscBool   isschur, isshell;

955:     PCShellGetContext(pc, &pc_ctx);
956:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
957:     PetscObjectTypeCompare((PetscObject)pc_ctx->S_j, MATSCHURCOMPLEMENT, &isschur);
958:     if (isschur) {
959:       PetscViewerASCIIPrintf(viewer, "  Dirichlet preconditioner (just from rank 0)\n");
960:     } else {
961:       PetscViewerASCIIPrintf(viewer, "  Lumped preconditioner (just from rank 0)\n");
962:     }
963:     PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer);
964:     if (rank == 0) {
965:       PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO);
966:       PetscViewerASCIIPushTab(sviewer);
967:       MatView(pc_ctx->S_j, sviewer);
968:       PetscViewerASCIIPopTab(sviewer);
969:       PetscViewerPopFormat(sviewer);
970:     }
971:     PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer);
972:     PetscObjectTypeCompare((PetscObject)pc_ctx->B_Ddelta, MATSHELL, &isshell);
973:     if (isshell) {
974:       BDdelta_DN ctx;
975:       PetscViewerASCIIPrintf(viewer, "  FETI-DP BDdelta: DB^t * (B D^-1 B^t)^-1 for deluxe scaling (just from rank 0)\n");
976:       MatShellGetContext(pc_ctx->B_Ddelta, &ctx);
977:       PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer);
978:       if (rank == 0) {
979:         PetscInt tl;

981:         PetscViewerASCIIGetTab(sviewer, &tl);
982:         PetscObjectSetTabLevel((PetscObject)ctx->kBD, tl);
983:         KSPView(ctx->kBD, sviewer);
984:         PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO);
985:         MatView(ctx->BD, sviewer);
986:         PetscViewerPopFormat(sviewer);
987:       }
988:       PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer);
989:     }
990:     PetscViewerFlush(viewer);
991:   }
992:   return 0;
993: }