Actual source code: bddc.c

petsc-dev 2014-08-28
Report Typos and Errors
  1: /* TODOLIST

  3:    ConstraintsSetup
  4:    - tolerances for constraints as an option (take care of single precision!)

  6:    Solvers
  7:    - Add support for cholesky for coarse solver (similar to local solvers)
  8:    - Propagate ksp prefixes for solvers to mat objects?
  9:    - Propagate nearnullspace info among levels

 11:    User interface
 12:    - ** DofSplitting and DM attached to pc?

 14:    Debugging output
 15:    - * Better management of verbosity levels of debugging output

 17:    Build
 18:    - make runexe59

 20:    Extra
 21:    - ** GetRid of PCBDDCApplySchur, use MatSchur instead
 22:    - *** Is it possible to work with PCBDDCGraph on boundary indices only (less memory consumed)?
 23:    - add support for computing h,H and related using coordinates?
 24:    - Change of basis approach does not work with my nonlinear mechanics example. why? (seems not an issue with l2gmap)
 25:    - Better management in PCIS code
 26:    - BDDC with MG framework?

 28:    FETIDP
 29:    - Move FETIDP code to its own classes

 31:    MATIS related operations contained in BDDC code
 32:    - Provide general case for subassembling
 33:    - *** Preallocation routines in MatISGetMPIAXAIJ

 35: */

 37: /* ----------------------------------------------------------------------------------------------------------------------------------------------
 38:    Implementation of BDDC preconditioner based on:
 39:    C. Dohrmann "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007
 40:    ---------------------------------------------------------------------------------------------------------------------------------------------- */

 42: #include <../src/ksp/pc/impls/bddc/bddc.h> /*I "petscpc.h" I*/  /* includes for fortran wrappers */
 43: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
 44: #include <petscblaslapack.h>

 46: /* -------------------------------------------------------------------------- */
 49: PetscErrorCode PCSetFromOptions_BDDC(PC pc)
 50: {
 51:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;

 55:   PetscOptionsHead("BDDC options");
 56:   /* Verbose debugging */
 57:   PetscOptionsInt("-pc_bddc_check_level","Verbose output for PCBDDC (intended for debug)","none",pcbddc->dbg_flag,&pcbddc->dbg_flag,NULL);
 58:   /* Primal space cumstomization */
 59:   PetscOptionsBool("-pc_bddc_use_local_mat_graph","Use or not adjacency graph of local mat for interface analysis","none",pcbddc->use_local_adj,&pcbddc->use_local_adj,NULL);
 60:   PetscOptionsBool("-pc_bddc_use_vertices","Use or not corner dofs in coarse space","none",pcbddc->use_vertices,&pcbddc->use_vertices,NULL);
 61:   PetscOptionsBool("-pc_bddc_use_edges","Use or not edge constraints in coarse space","none",pcbddc->use_edges,&pcbddc->use_edges,NULL);
 62:   PetscOptionsBool("-pc_bddc_use_faces","Use or not face constraints in coarse space","none",pcbddc->use_faces,&pcbddc->use_faces,NULL);
 63:   /* Change of basis */
 64:   PetscOptionsBool("-pc_bddc_use_change_of_basis","Use or not internal change of basis on local edge nodes","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);
 65:   PetscOptionsBool("-pc_bddc_use_change_on_faces","Use or not internal change of basis on local face nodes","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);
 66:   if (!pcbddc->use_change_of_basis) {
 67:     pcbddc->use_change_on_faces = PETSC_FALSE;
 68:   }
 69:   /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */
 70:   PetscOptionsBool("-pc_bddc_switch_static","Switch on static condensation ops around the interface preconditioner","none",pcbddc->switch_static,&pcbddc->switch_static,NULL);
 71:   PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);
 72:   PetscOptionsInt("-pc_bddc_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);
 73:   PetscOptionsBool("-pc_bddc_use_coarse_estimates","Use estimated eigenvalues for coarse problem","none",pcbddc->use_coarse_estimates,&pcbddc->use_coarse_estimates,NULL);
 74:   PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);
 75:   PetscOptionsInt("-pc_bddc_deluxe_threshold","Deluxe subproblems (Schur principal minors) smaller than this value are explicilty computed (-1 computes all)","none",pcbddc->deluxe_threshold,&pcbddc->deluxe_threshold,NULL);
 76:   PetscOptionsBool("-pc_bddc_deluxe_rebuild","Whether or not the interface graph for deluxe has to be rebuilt (i.e. use a standard definition of the interface)","none",pcbddc->deluxe_rebuild,&pcbddc->deluxe_rebuild,NULL);
 77:   PetscOptionsInt("-pc_bddc_deluxe_layers","Number of dofs' layers for the application of deluxe cheap version (i.e. -1 uses all dofs)","none",pcbddc->deluxe_layers,&pcbddc->deluxe_layers,NULL);
 78:   PetscOptionsBool("-pc_bddc_deluxe_use_useradj","Whether or not the CSR graph specified by the user should be used for computing layers (default is to use adj of local mat)","none",pcbddc->deluxe_use_useradj,&pcbddc->deluxe_use_useradj,NULL);
 79:   PetscOptionsTail();
 80:   return(0);
 81: }
 82: /* -------------------------------------------------------------------------- */
 85: static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change)
 86: {
 87:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;

 91:   MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);
 92:   PetscObjectReference((PetscObject)change);
 93:   pcbddc->user_ChangeOfBasisMatrix = change;
 94:   return(0);
 95: }
 98: /*@
 99:  PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs

101:    Collective on PC

103:    Input Parameters:
104: +  pc - the preconditioning context
105: -  change - the change of basis matrix

107:    Level: intermediate

109:    Notes:

111: .seealso: PCBDDC
112: @*/
113: PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change)
114: {

121:   if (pc->mat) {
122:     PetscInt rows_c,cols_c,rows,cols;
123:     MatGetSize(pc->mat,&rows,&cols);
124:     MatGetSize(change,&rows_c,&cols_c);
125:     if (rows_c != rows) {
126:       SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of rows for change of basis matrix! %d != %d",rows_c,rows);
127:     }
128:     if (cols_c != cols) {
129:       SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of columns for change of basis matrix! %d != %d",cols_c,cols);
130:     }
131:     MatGetLocalSize(pc->mat,&rows,&cols);
132:     MatGetLocalSize(change,&rows_c,&cols_c);
133:     if (rows_c != rows) {
134:       SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local rows for change of basis matrix! %d != %d",rows_c,rows);
135:     }
136:     if (cols_c != cols) {
137:       SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local columns for change of basis matrix! %d != %d",cols_c,cols);
138:     }
139:   }
140:   PetscTryMethod(pc,"PCBDDCSetChangeOfBasisMat_C",(PC,Mat),(pc,change));
141:   return(0);
142: }
143: /* -------------------------------------------------------------------------- */
146: static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
147: {
148:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;

152:   ISDestroy(&pcbddc->user_primal_vertices);
153:   PetscObjectReference((PetscObject)PrimalVertices);
154:   pcbddc->user_primal_vertices = PrimalVertices;
155:   return(0);
156: }
159: /*@
160:  PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC

162:    Not collective

164:    Input Parameters:
165: +  pc - the preconditioning context
166: -  PrimalVertices - index set of primal vertices in local numbering

168:    Level: intermediate

170:    Notes:

172: .seealso: PCBDDC
173: @*/
174: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
175: {

181:   PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));
182:   return(0);
183: }
184: /* -------------------------------------------------------------------------- */
187: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k)
188: {
189:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

192:   pcbddc->coarsening_ratio = k;
193:   return(0);
194: }

198: /*@
199:  PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel

201:    Logically collective on PC

203:    Input Parameters:
204: +  pc - the preconditioning context
205: -  k - coarsening ratio (H/h at the coarser level)

207:    Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level

209:    Level: intermediate

211:    Notes:

213: .seealso: PCBDDC
214: @*/
215: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k)
216: {

222:   PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));
223:   return(0);
224: }

226: /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
229: static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc,PetscBool flg)
230: {
231:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

234:   pcbddc->use_exact_dirichlet_trick = flg;
235:   return(0);
236: }

240: PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg)
241: {

247:   PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));
248:   return(0);
249: }

253: static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level)
254: {
255:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

258:   pcbddc->current_level = level;
259:   return(0);
260: }

264: PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level)
265: {

271:   PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));
272:   return(0);
273: }

277: static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels)
278: {
279:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

282:   pcbddc->max_levels = levels;
283:   return(0);
284: }

288: /*@
289:  PCBDDCSetLevels - Sets the maximum number of levels for multilevel

291:    Logically collective on PC

293:    Input Parameters:
294: +  pc - the preconditioning context
295: -  levels - the maximum number of levels (max 9)

297:    Default value is 0, i.e. traditional one-level BDDC

299:    Level: intermediate

301:    Notes:

303: .seealso: PCBDDC
304: @*/
305: PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels)
306: {

312:   if (levels > 99) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Maximum number of levels for bddc is 99\n");
313:   PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));
314:   return(0);
315: }
316: /* -------------------------------------------------------------------------- */

320: static PetscErrorCode PCBDDCSetNullSpace_BDDC(PC pc,MatNullSpace NullSpace)
321: {
322:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

326:   PetscObjectReference((PetscObject)NullSpace);
327:   MatNullSpaceDestroy(&pcbddc->NullSpace);
328:   pcbddc->NullSpace = NullSpace;
329:   return(0);
330: }

334: /*@
335:  PCBDDCSetNullSpace - Set nullspace for BDDC operator

337:    Logically collective on PC and MatNullSpace

339:    Input Parameters:
340: +  pc - the preconditioning context
341: -  NullSpace - Null space of the linear operator to be preconditioned (Pmat)

343:    Level: intermediate

345:    Notes:

347: .seealso: PCBDDC
348: @*/
349: PetscErrorCode PCBDDCSetNullSpace(PC pc,MatNullSpace NullSpace)
350: {

357:   PetscTryMethod(pc,"PCBDDCSetNullSpace_C",(PC,MatNullSpace),(pc,NullSpace));
358:   return(0);
359: }
360: /* -------------------------------------------------------------------------- */

364: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries)
365: {
366:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

370:   /* last user setting takes precendence -> destroy any other customization */
371:   ISDestroy(&pcbddc->DirichletBoundariesLocal);
372:   ISDestroy(&pcbddc->DirichletBoundaries);
373:   PetscObjectReference((PetscObject)DirichletBoundaries);
374:   pcbddc->DirichletBoundaries = DirichletBoundaries;
375:   pcbddc->recompute_topography = PETSC_TRUE;
376:   return(0);
377: }

381: /*@
382:  PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem.

384:    Collective

386:    Input Parameters:
387: +  pc - the preconditioning context
388: -  DirichletBoundaries - parallel IS defining the Dirichlet boundaries

390:    Level: intermediate

392:    Notes: Any process can list any global node

394: .seealso: PCBDDC
395: @*/
396: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries)
397: {

404:   PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));
405:   return(0);
406: }
407: /* -------------------------------------------------------------------------- */

411: static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc,IS DirichletBoundaries)
412: {
413:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

417:   /* last user setting takes precendence -> destroy any other customization */
418:   ISDestroy(&pcbddc->DirichletBoundariesLocal);
419:   ISDestroy(&pcbddc->DirichletBoundaries);
420:   PetscObjectReference((PetscObject)DirichletBoundaries);
421:   pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
422:   pcbddc->recompute_topography = PETSC_TRUE;
423:   return(0);
424: }

428: /*@
429:  PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering.

431:    Collective

433:    Input Parameters:
434: +  pc - the preconditioning context
435: -  DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering)

437:    Level: intermediate

439:    Notes:

441: .seealso: PCBDDC
442: @*/
443: PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries)
444: {

451:   PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));
452:   return(0);
453: }
454: /* -------------------------------------------------------------------------- */

458: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries)
459: {
460:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

464:   /* last user setting takes precendence -> destroy any other customization */
465:   ISDestroy(&pcbddc->NeumannBoundariesLocal);
466:   ISDestroy(&pcbddc->NeumannBoundaries);
467:   PetscObjectReference((PetscObject)NeumannBoundaries);
468:   pcbddc->NeumannBoundaries = NeumannBoundaries;
469:   pcbddc->recompute_topography = PETSC_TRUE;
470:   return(0);
471: }

475: /*@
476:  PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem.

478:    Collective

480:    Input Parameters:
481: +  pc - the preconditioning context
482: -  NeumannBoundaries - parallel IS defining the Neumann boundaries

484:    Level: intermediate

486:    Notes: Any process can list any global node

488: .seealso: PCBDDC
489: @*/
490: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries)
491: {

498:   PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));
499:   return(0);
500: }
501: /* -------------------------------------------------------------------------- */

505: static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries)
506: {
507:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

511:   /* last user setting takes precendence -> destroy any other customization */
512:   ISDestroy(&pcbddc->NeumannBoundariesLocal);
513:   ISDestroy(&pcbddc->NeumannBoundaries);
514:   PetscObjectReference((PetscObject)NeumannBoundaries);
515:   pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
516:   pcbddc->recompute_topography = PETSC_TRUE;
517:   return(0);
518: }

522: /*@
523:  PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering.

525:    Collective

527:    Input Parameters:
528: +  pc - the preconditioning context
529: -  NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering)

531:    Level: intermediate

533:    Notes:

535: .seealso: PCBDDC
536: @*/
537: PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries)
538: {

545:   PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));
546:   return(0);
547: }
548: /* -------------------------------------------------------------------------- */

552: static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries)
553: {
554:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

557:   *DirichletBoundaries = pcbddc->DirichletBoundaries;
558:   return(0);
559: }

563: /*@
564:  PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries

566:    Collective

568:    Input Parameters:
569: .  pc - the preconditioning context

571:    Output Parameters:
572: .  DirichletBoundaries - index set defining the Dirichlet boundaries

574:    Level: intermediate

576:    Notes: The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries

578: .seealso: PCBDDC
579: @*/
580: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries)
581: {

586:   PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));
587:   return(0);
588: }
589: /* -------------------------------------------------------------------------- */

593: static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries)
594: {
595:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

598:   *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
599:   return(0);
600: }

604: /*@
605:  PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering)

607:    Collective

609:    Input Parameters:
610: .  pc - the preconditioning context

612:    Output Parameters:
613: .  DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries

615:    Level: intermediate

617:    Notes:

619: .seealso: PCBDDC
620: @*/
621: PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries)
622: {

627:   PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));
628:   return(0);
629: }
630: /* -------------------------------------------------------------------------- */

634: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries)
635: {
636:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

639:   *NeumannBoundaries = pcbddc->NeumannBoundaries;
640:   return(0);
641: }

645: /*@
646:  PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries

648:    Collective

650:    Input Parameters:
651: .  pc - the preconditioning context

653:    Output Parameters:
654: .  NeumannBoundaries - index set defining the Neumann boundaries

656:    Level: intermediate

658:    Notes: The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries

660: .seealso: PCBDDC
661: @*/
662: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries)
663: {

668:   PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));
669:   return(0);
670: }
671: /* -------------------------------------------------------------------------- */

675: static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries)
676: {
677:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

680:   *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
681:   return(0);
682: }

686: /*@
687:  PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering)

689:    Collective

691:    Input Parameters:
692: .  pc - the preconditioning context

694:    Output Parameters:
695: .  NeumannBoundaries - index set defining the subdomain part of Neumann boundaries

697:    Level: intermediate

699:    Notes:

701: .seealso: PCBDDC
702: @*/
703: PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries)
704: {

709:   PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));
710:   return(0);
711: }
712: /* -------------------------------------------------------------------------- */

716: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
717: {
718:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
719:   PCBDDCGraph    mat_graph = pcbddc->mat_graph;

723:   /* free old CSR */
724:   PCBDDCGraphResetCSR(mat_graph);
725:   /* TODO: PCBDDCGraphSetAdjacency */
726:   /* get CSR into graph structure */
727:   if (copymode == PETSC_COPY_VALUES) {
728:     PetscMalloc1((nvtxs+1),&mat_graph->xadj);
729:     PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);
730:     PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));
731:     PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));
732:   } else if (copymode == PETSC_OWN_POINTER) {
733:     mat_graph->xadj = (PetscInt*)xadj;
734:     mat_graph->adjncy = (PetscInt*)adjncy;
735:   }
736:   mat_graph->nvtxs_csr = nvtxs;
737:   return(0);
738: }

742: /*@
743:  PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local Neumann matrix

745:    Not collective

747:    Input Parameters:
748: +  pc - the preconditioning context
749: .  nvtxs - number of local vertices of the graph (i.e., the local size of your problem)
750: .  xadj, adjncy - the CSR graph
751: -  copymode - either PETSC_COPY_VALUES or PETSC_OWN_POINTER.

753:    Level: intermediate

755:    Notes:

757: .seealso: PCBDDC,PetscCopyMode
758: @*/
759: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
760: {
761:   void (*f)(void) = 0;

768:   if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) {
769:     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d in %s\n",copymode,__FUNCT__);
770:   }
771:   PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));
772:   /* free arrays if PCBDDC is not the PC type */
773:   PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);
774:   if (!f && copymode == PETSC_OWN_POINTER) {
775:     PetscFree(xadj);
776:     PetscFree(adjncy);
777:   }
778:   return(0);
779: }
780: /* -------------------------------------------------------------------------- */

784: static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[])
785: {
786:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
787:   PetscInt i;

791:   /* Destroy ISes if they were already set */
792:   for (i=0;i<pcbddc->n_ISForDofsLocal;i++) {
793:     ISDestroy(&pcbddc->ISForDofsLocal[i]);
794:   }
795:   PetscFree(pcbddc->ISForDofsLocal);
796:   /* last user setting takes precendence -> destroy any other customization */
797:   for (i=0;i<pcbddc->n_ISForDofs;i++) {
798:     ISDestroy(&pcbddc->ISForDofs[i]);
799:   }
800:   PetscFree(pcbddc->ISForDofs);
801:   pcbddc->n_ISForDofs = 0;
802:   /* allocate space then set */
803:   PetscMalloc(n_is*sizeof(IS),&pcbddc->ISForDofsLocal);
804:   for (i=0;i<n_is;i++) {
805:     PetscObjectReference((PetscObject)ISForDofs[i]);
806:     pcbddc->ISForDofsLocal[i]=ISForDofs[i];
807:   }
808:   pcbddc->n_ISForDofsLocal=n_is;
809:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
810:   pcbddc->recompute_topography = PETSC_TRUE;
811:   return(0);
812: }

816: /*@
817:  PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix

819:    Collective

821:    Input Parameters:
822: +  pc - the preconditioning context
823: -  n_is - number of index sets defining the fields
824: .  ISForDofs - array of IS describing the fields in local ordering

826:    Level: intermediate

828:    Notes: n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to a different field.

830: .seealso: PCBDDC
831: @*/
832: PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[])
833: {
834:   PetscInt       i;

840:   for (i=0;i<n_is;i++) {
843:   }
844:   PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));
845:   return(0);
846: }
847: /* -------------------------------------------------------------------------- */

851: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[])
852: {
853:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
854:   PetscInt i;

858:   /* Destroy ISes if they were already set */
859:   for (i=0;i<pcbddc->n_ISForDofs;i++) {
860:     ISDestroy(&pcbddc->ISForDofs[i]);
861:   }
862:   PetscFree(pcbddc->ISForDofs);
863:   /* last user setting takes precendence -> destroy any other customization */
864:   for (i=0;i<pcbddc->n_ISForDofsLocal;i++) {
865:     ISDestroy(&pcbddc->ISForDofsLocal[i]);
866:   }
867:   PetscFree(pcbddc->ISForDofsLocal);
868:   pcbddc->n_ISForDofsLocal = 0;
869:   /* allocate space then set */
870:   PetscMalloc1(n_is,&pcbddc->ISForDofs);
871:   for (i=0;i<n_is;i++) {
872:     PetscObjectReference((PetscObject)ISForDofs[i]);
873:     pcbddc->ISForDofs[i]=ISForDofs[i];
874:   }
875:   pcbddc->n_ISForDofs=n_is;
876:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
877:   pcbddc->recompute_topography = PETSC_TRUE;
878:   return(0);
879: }

883: /*@
884:  PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix

886:    Collective

888:    Input Parameters:
889: +  pc - the preconditioning context
890: -  n_is - number of index sets defining the fields
891: .  ISForDofs - array of IS describing the fields in global ordering

893:    Level: intermediate

895:    Notes: Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to a different field.

897: .seealso: PCBDDC
898: @*/
899: PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[])
900: {
901:   PetscInt       i;

907:   for (i=0;i<n_is;i++) {
910:   }
911:   PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));
912:   return(0);
913: }

915: /* -------------------------------------------------------------------------- */
918: /* -------------------------------------------------------------------------- */
919: /*
920:    PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial
921:                      guess if a transformation of basis approach has been selected.

923:    Input Parameter:
924: +  pc - the preconditioner contex

926:    Application Interface Routine: PCPreSolve()

928:    Notes:
929:    The interface routine PCPreSolve() is not usually called directly by
930:    the user, but instead is called by KSPSolve().
931: */
932: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
933: {
935:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
936:   PC_IS          *pcis = (PC_IS*)(pc->data);
937:   IS             dirIS;
938:   Vec            used_vec;

941:   /* if we are working with cg, one dirichlet solve can be avoided during Krylov iterations */
942:   if (ksp) {
943:     PetscBool iscg;
944:     PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);
945:     if (!iscg) {
946:       PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);
947:     }
948:   }
949:   /* Creates parallel work vectors used in presolve */
950:   if (!pcbddc->original_rhs) {
951:     VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);
952:   }
953:   if (!pcbddc->temp_solution) {
954:     VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);
955:   }
956:   if (x) {
957:     PetscObjectReference((PetscObject)x);
958:     used_vec = x;
959:   } else {
960:     PetscObjectReference((PetscObject)pcbddc->temp_solution);
961:     used_vec = pcbddc->temp_solution;
962:     VecSet(used_vec,0.0);
963:   }
964:   /* hack into ksp data structure PCPreSolve comes earlier in src/ksp/ksp/interface/itfunc.c */
965:   if (ksp) {
966:     PetscBool guess_nonzero;
967:     KSPGetInitialGuessNonzero(ksp,&guess_nonzero);
968:     if (!guess_nonzero) {
969:       VecSet(used_vec,0.0);
970:     }
971:   }

973:   /* store the original rhs */
974:   VecCopy(rhs,pcbddc->original_rhs);

976:   /* Take into account zeroed rows -> change rhs and store solution removed */
977:   /* note that Dirichlet boundaries in global ordering (if any) has already been translated into local ordering in PCBDDCAnalyzeInterface */
978:   PCBDDCGetDirichletBoundariesLocal(pc,&dirIS);
979:   if (rhs && dirIS) {
980:     Mat_IS      *matis = (Mat_IS*)pc->pmat->data;
981:     PetscInt    dirsize,i,*is_indices;
982:     PetscScalar *array_x,*array_diagonal;

984:     MatGetDiagonal(pc->pmat,pcis->vec1_global);
985:     VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);
986:     VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);
987:     VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);
988:     VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
989:     VecScatterEnd(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
990:     ISGetLocalSize(dirIS,&dirsize);
991:     VecGetArray(pcis->vec1_N,&array_x);
992:     VecGetArray(pcis->vec2_N,&array_diagonal);
993:     ISGetIndices(dirIS,(const PetscInt**)&is_indices);
994:     for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
995:     ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);
996:     VecRestoreArray(pcis->vec2_N,&array_diagonal);
997:     VecRestoreArray(pcis->vec1_N,&array_x);
998:     VecScatterBegin(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);
999:     VecScatterEnd(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);

1001:     /* remove the computed solution from the rhs */
1002:     VecScale(used_vec,-1.0);
1003:     MatMultAdd(pc->pmat,used_vec,rhs,rhs);
1004:     VecScale(used_vec,-1.0);
1005:   }

1007:   /* store partially computed solution and set initial guess */
1008:   if (x) {
1009:     VecCopy(used_vec,pcbddc->temp_solution);
1010:     VecSet(used_vec,0.0);
1011:     if (pcbddc->use_exact_dirichlet_trick) {
1012:       VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1013:       VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1014:       KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);
1015:       VecScatterBegin(pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);
1016:       VecScatterEnd(pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);
1017:       if (ksp) {
1018:         KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
1019:       }
1020:     }
1021:   }

1023:   if (pcbddc->ChangeOfBasisMatrix) {
1024:     PCBDDCChange_ctx change_ctx;

1026:     /* get change ctx */
1027:     MatShellGetContext(pcbddc->new_global_mat,&change_ctx);

1029:     /* set current iteration matrix inside change context (change of basis has been already set into the ctx during PCSetUp) */
1030:     MatDestroy(&change_ctx->original_mat);
1031:     PetscObjectReference((PetscObject)pc->mat);
1032:     change_ctx->original_mat = pc->mat;

1034:     /* change iteration matrix */
1035:     MatDestroy(&pc->mat);
1036:     PetscObjectReference((PetscObject)pcbddc->new_global_mat);
1037:     pc->mat = pcbddc->new_global_mat;

1039:     /* change rhs */
1040:     MatMultTranspose(change_ctx->global_change,rhs,pcis->vec1_global);
1041:     VecCopy(pcis->vec1_global,rhs);
1042:   }

1044:   /* remove nullspace if present */
1045:   if (ksp && pcbddc->NullSpace) {
1046:     MatNullSpaceRemove(pcbddc->NullSpace,used_vec);
1047:     MatNullSpaceRemove(pcbddc->NullSpace,rhs);
1048:   }
1049:   VecDestroy(&used_vec);
1050:   return(0);
1051: }

1053: /* -------------------------------------------------------------------------- */
1056: /* -------------------------------------------------------------------------- */
1057: /*
1058:    PCPostSolve_BDDC - Changes the computed solution if a transformation of basis
1059:                      approach has been selected. Also, restores rhs to its original state.

1061:    Input Parameter:
1062: +  pc - the preconditioner contex

1064:    Application Interface Routine: PCPostSolve()

1066:    Notes:
1067:    The interface routine PCPostSolve() is not usually called directly by
1068:    the user, but instead is called by KSPSolve().
1069: */
1070: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1071: {
1073:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;

1076:   if (pcbddc->ChangeOfBasisMatrix) {
1077:     PCBDDCChange_ctx change_ctx;

1079:     /* get change ctx */
1080:     MatShellGetContext(pcbddc->new_global_mat,&change_ctx);

1082:     /* restore iteration matrix */
1083:     MatDestroy(&pc->mat);
1084:     PetscObjectReference((PetscObject)change_ctx->original_mat);
1085:     pc->mat = change_ctx->original_mat;

1087:     /* get solution in original basis */
1088:     if (x) {
1089:       PC_IS *pcis = (PC_IS*)(pc->data);
1090:       MatMult(change_ctx->global_change,x,pcis->vec1_global);
1091:       VecCopy(pcis->vec1_global,x);
1092:     }
1093:   }

1095:   /* add solution removed in presolve */
1096:   if (x) {
1097:     VecAXPY(x,1.0,pcbddc->temp_solution);
1098:   }

1100:   /* restore rhs to its original state */
1101:   if (rhs) {
1102:     VecCopy(pcbddc->original_rhs,rhs);
1103:   }
1104:   return(0);
1105: }
1106: /* -------------------------------------------------------------------------- */
1109: /* -------------------------------------------------------------------------- */
1110: /*
1111:    PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner
1112:                   by setting data structures and options.

1114:    Input Parameter:
1115: +  pc - the preconditioner context

1117:    Application Interface Routine: PCSetUp()

1119:    Notes:
1120:    The interface routine PCSetUp() is not usually called directly by
1121:    the user, but instead is called by PCApply() if necessary.
1122: */
1123: PetscErrorCode PCSetUp_BDDC(PC pc)
1124: {
1125:   PetscErrorCode   ierr;
1126:   PC_BDDC*         pcbddc = (PC_BDDC*)pc->data;
1127:   MatNullSpace     nearnullspace;
1128:   PetscBool        computeis,computetopography,computesolvers;
1129:   PetscBool        new_nearnullspace_provided;

1132:   /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
1133:   /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
1134:      Also, BDDC directly build the Dirichlet problem */
1135:   /* split work */
1136:   if (pc->setupcalled) {
1137:     computeis = PETSC_FALSE;
1138:     if (pc->flag == SAME_NONZERO_PATTERN) {
1139:       computetopography = PETSC_FALSE;
1140:       computesolvers = PETSC_TRUE;
1141:     } else { /* DIFFERENT_NONZERO_PATTERN */
1142:       computetopography = PETSC_TRUE;
1143:       computesolvers = PETSC_TRUE;
1144:     }
1145:   } else {
1146:     computeis = PETSC_TRUE;
1147:     computetopography = PETSC_TRUE;
1148:     computesolvers = PETSC_TRUE;
1149:   }
1150:   if (pcbddc->recompute_topography) {
1151:     computetopography = PETSC_TRUE;
1152:   }

1154:   /* Get stdout for dbg */
1155:   if (pcbddc->dbg_flag) {
1156:     if (!pcbddc->dbg_viewer) {
1157:       pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
1158:       PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);
1159:     }
1160:     PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);
1161:   }

1163:   /* Set up all the "iterative substructuring" common block without computing solvers */
1164:   if (computeis) {
1165:     /* HACK INTO PCIS */
1166:     PC_IS* pcis = (PC_IS*)pc->data;
1167:     pcis->computesolvers = PETSC_FALSE;
1168:     PCISSetUp(pc);
1169:     ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcbddc->BtoNmap);
1170:   }

1172:   /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
1173:   if (pcbddc->user_ChangeOfBasisMatrix) {
1174:     pcbddc->use_change_of_basis = PETSC_FALSE;
1175:   }

1177:   /* Analyze interface */
1178:   if (computetopography) {
1179:     PCBDDCAnalyzeInterface(pc);
1180:     /* Schurs on subsets should be reset */
1181:     if (pcbddc->deluxe_ctx) {
1182:       PCBDDCSubSchursReset(pcbddc->deluxe_ctx->sub_schurs);
1183:     }
1184:   }

1186:   /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
1187:   new_nearnullspace_provided = PETSC_FALSE;
1188:   MatGetNearNullSpace(pc->pmat,&nearnullspace);
1189:   if (pcbddc->onearnullspace) { /* already used nearnullspace */
1190:     if (!nearnullspace) { /* near null space attached to mat has been destroyed */
1191:       new_nearnullspace_provided = PETSC_TRUE;
1192:     } else {
1193:       /* determine if the two nullspaces are different (should be lightweight) */
1194:       if (nearnullspace != pcbddc->onearnullspace) {
1195:         new_nearnullspace_provided = PETSC_TRUE;
1196:       } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
1197:         PetscInt         i;
1198:         const Vec        *nearnullvecs;
1199:         PetscObjectState state;
1200:         PetscInt         nnsp_size;
1201:         MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);
1202:         for (i=0;i<nnsp_size;i++) {
1203:           PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);
1204:           if (pcbddc->onearnullvecs_state[i] != state) {
1205:             new_nearnullspace_provided = PETSC_TRUE;
1206:             break;
1207:           }
1208:         }
1209:       }
1210:     }
1211:   } else {
1212:     if (!nearnullspace) { /* both nearnullspaces are null */
1213:       new_nearnullspace_provided = PETSC_FALSE;
1214:     } else { /* nearnullspace attached later */
1215:       new_nearnullspace_provided = PETSC_TRUE;
1216:     }
1217:   }

1219:   /* Setup constraints and related work vectors */
1220:   /* reset primal space flags */
1221:   pcbddc->new_primal_space = PETSC_FALSE;
1222:   pcbddc->new_primal_space_local = PETSC_FALSE;
1223:   if (computetopography || new_nearnullspace_provided) {
1224:     /* It also sets the primal space flags */
1225:     PCBDDCConstraintsSetUp(pc);
1226:     /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
1227:     PCBDDCSetUpLocalWorkVectors(pc);
1228:     /* Schurs on subsets should be reset */
1229:     if (pcbddc->deluxe_ctx) {
1230:       PCBDDCSubSchursReset(pcbddc->deluxe_ctx->sub_schurs);
1231:     }
1232:   }

1234:   if (computesolvers || pcbddc->new_primal_space) {
1235:     /* Create coarse and local stuffs */
1236:     PCBDDCSetUpSolvers(pc);
1237:     /* Create scaling operators */
1238:     PCBDDCScalingSetUp(pc);
1239:   }

1241:   if (pcbddc->dbg_flag) {
1242:     PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);
1243:   }
1244:   return(0);
1245: }

1247: /* -------------------------------------------------------------------------- */
1248: /*
1249:    PCApply_BDDC - Applies the BDDC operator to a vector.

1251:    Input Parameters:
1252: .  pc - the preconditioner context
1253: .  r - input vector (global)

1255:    Output Parameter:
1256: .  z - output vector (global)

1258:    Application Interface Routine: PCApply()
1259:  */
1262: PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z)
1263: {
1264:   PC_IS             *pcis = (PC_IS*)(pc->data);
1265:   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
1266:   PetscErrorCode    ierr;
1267:   const PetscScalar one = 1.0;
1268:   const PetscScalar m_one = -1.0;
1269:   const PetscScalar zero = 0.0;

1271: /* This code is similar to that provided in nn.c for PCNN
1272:    NN interface preconditioner changed to BDDC
1273:    Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static = PETSC_TRUE) */

1276:   if (!pcbddc->use_exact_dirichlet_trick) {
1277:     /* First Dirichlet solve */
1278:     VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1279:     VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1280:     KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);
1281:     /*
1282:       Assembling right hand side for BDDC operator
1283:       - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE)
1284:       - pcis->vec1_B the interface part of the global vector z
1285:     */
1286:     VecScale(pcis->vec2_D,m_one);
1287:     MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);
1288:     if (pcbddc->switch_static) { MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D); }
1289:     VecScale(pcis->vec2_D,m_one);
1290:     VecCopy(r,z);
1291:     VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);
1292:     VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);
1293:     PCBDDCScalingRestriction(pc,z,pcis->vec1_B);
1294:   } else {
1295:     VecSet(pcis->vec1_D,zero);
1296:     VecSet(pcis->vec2_D,zero);
1297:     PCBDDCScalingRestriction(pc,r,pcis->vec1_B);
1298:   }

1300:   /* Apply interface preconditioner
1301:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1302:   PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);

1304:   /* Apply transpose of partition of unity operator */
1305:   PCBDDCScalingExtension(pc,pcis->vec1_B,z);

1307:   /* Second Dirichlet solve and assembling of output */
1308:   VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1309:   VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1310:   MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);
1311:   if (pcbddc->switch_static) { MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D); }
1312:   KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);
1313:   VecScale(pcis->vec4_D,m_one);
1314:   if (pcbddc->switch_static) { VecAXPY(pcis->vec4_D,one,pcis->vec1_D); }
1315:   VecAXPY(pcis->vec2_D,one,pcis->vec4_D);
1316:   VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
1317:   VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
1318:   return(0);
1319: }

1321: /* -------------------------------------------------------------------------- */
1322: /*
1323:    PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector.

1325:    Input Parameters:
1326: .  pc - the preconditioner context
1327: .  r - input vector (global)

1329:    Output Parameter:
1330: .  z - output vector (global)

1332:    Application Interface Routine: PCApplyTranspose()
1333:  */
1336: PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z)
1337: {
1338:   PC_IS             *pcis = (PC_IS*)(pc->data);
1339:   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
1340:   PetscErrorCode    ierr;
1341:   const PetscScalar one = 1.0;
1342:   const PetscScalar m_one = -1.0;
1343:   const PetscScalar zero = 0.0;

1346:   if (!pcbddc->use_exact_dirichlet_trick) {
1347:     /* First Dirichlet solve */
1348:     VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1349:     VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1350:     KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);
1351:     /*
1352:       Assembling right hand side for BDDC operator
1353:       - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE)
1354:       - pcis->vec1_B the interface part of the global vector z
1355:     */
1356:     VecScale(pcis->vec2_D,m_one);
1357:     MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);
1358:     if (pcbddc->switch_static) { MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D); }
1359:     VecScale(pcis->vec2_D,m_one);
1360:     VecCopy(r,z);
1361:     VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);
1362:     VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);
1363:     PCBDDCScalingRestriction(pc,z,pcis->vec1_B);
1364:   } else {
1365:     VecSet(pcis->vec1_D,zero);
1366:     VecSet(pcis->vec2_D,zero);
1367:     PCBDDCScalingRestriction(pc,r,pcis->vec1_B);
1368:   }

1370:   /* Apply interface preconditioner
1371:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1372:   PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);

1374:   /* Apply transpose of partition of unity operator */
1375:   PCBDDCScalingExtension(pc,pcis->vec1_B,z);

1377:   /* Second Dirichlet solve and assembling of output */
1378:   VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1379:   VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1380:   MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);
1381:   if (pcbddc->switch_static) { MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D); }
1382:   KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);
1383:   VecScale(pcis->vec4_D,m_one);
1384:   if (pcbddc->switch_static) { VecAXPY(pcis->vec4_D,one,pcis->vec1_D); }
1385:   VecAXPY(pcis->vec2_D,one,pcis->vec4_D);
1386:   VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
1387:   VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
1388:   return(0);
1389: }
1390: /* -------------------------------------------------------------------------- */

1394: PetscErrorCode PCDestroy_BDDC(PC pc)
1395: {
1396:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;

1400:   /* free data created by PCIS */
1401:   PCISDestroy(pc);
1402:   /* free BDDC custom data  */
1403:   PCBDDCResetCustomization(pc);
1404:   /* destroy objects related to topography */
1405:   PCBDDCResetTopography(pc);
1406:   /* free allocated graph structure */
1407:   PetscFree(pcbddc->mat_graph);
1408:   /* destroy objects for scaling operator */
1409:   PCBDDCScalingDestroy(pc);
1410:   PetscFree(pcbddc->deluxe_ctx);
1411:   /* free solvers stuff */
1412:   PCBDDCResetSolvers(pc);
1413:   /* free global vectors needed in presolve */
1414:   VecDestroy(&pcbddc->temp_solution);
1415:   VecDestroy(&pcbddc->original_rhs);
1416:   /* free stuff for change of basis hooks */
1417:   if (pcbddc->new_global_mat) {
1418:     PCBDDCChange_ctx change_ctx;
1419:     MatShellGetContext(pcbddc->new_global_mat,&change_ctx);
1420:     MatDestroy(&change_ctx->original_mat);
1421:     MatDestroy(&change_ctx->global_change);
1422:     VecDestroyVecs(2,&change_ctx->work);
1423:     PetscFree(change_ctx);
1424:   }
1425:   MatDestroy(&pcbddc->new_global_mat);
1426:   /* remove map from local boundary to local numbering */
1427:   ISLocalToGlobalMappingDestroy(&pcbddc->BtoNmap);
1428:   /* remove functions */
1429:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);
1430:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);
1431:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);
1432:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);
1433:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);
1434:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);
1435:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);
1436:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);
1437:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);
1438:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);
1439:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);
1440:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);
1441:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);
1442:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);
1443:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);
1444:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);
1445:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);
1446:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);
1447:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);
1448:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);
1449:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);
1450:   /* Free the private data structure */
1451:   PetscFree(pc->data);
1452:   return(0);
1453: }
1454: /* -------------------------------------------------------------------------- */

1458: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
1459: {
1460:   FETIDPMat_ctx  mat_ctx;
1461:   PC_IS*         pcis;
1462:   PC_BDDC*       pcbddc;

1466:   MatShellGetContext(fetidp_mat,&mat_ctx);
1467:   pcis = (PC_IS*)mat_ctx->pc->data;
1468:   pcbddc = (PC_BDDC*)mat_ctx->pc->data;

1470:   /* change of basis for physical rhs if needed
1471:      It also changes the rhs in case of dirichlet boundaries */
1472:   PCPreSolve_BDDC(mat_ctx->pc,NULL,standard_rhs,NULL);
1473:   /* store vectors for computation of fetidp final solution */
1474:   VecScatterBegin(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);
1475:   VecScatterEnd(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);
1476:   /* scale rhs since it should be unassembled */
1477:   /* TODO use counter scaling? (also below) */
1478:   VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
1479:   VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
1480:   /* Apply partition of unity */
1481:   VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);
1482:   /* PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B); */
1483:   if (!pcbddc->switch_static) {
1484:     /* compute partially subassembled Schur complement right-hand side */
1485:     KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);
1486:     MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);
1487:     VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);
1488:     VecSet(standard_rhs,0.0);
1489:     VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);
1490:     VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);
1491:     /* PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B); */
1492:     VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
1493:     VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
1494:     VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);
1495:   }
1496:   /* BDDC rhs */
1497:   VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);
1498:   if (pcbddc->switch_static) {
1499:     VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);
1500:   }
1501:   /* apply BDDC */
1502:   PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);
1503:   /* Application of B_delta and assembling of rhs for fetidp fluxes */
1504:   VecSet(fetidp_flux_rhs,0.0);
1505:   MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);
1506:   VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);
1507:   VecScatterEnd  (mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);
1508:   /* restore original rhs */
1509:   VecCopy(pcbddc->original_rhs,standard_rhs);
1510:   return(0);
1511: }

1515: /*@
1516:  PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETIDP linear system

1518:    Collective

1520:    Input Parameters:
1521: +  fetidp_mat   - the FETIDP matrix object obtained by calling PCBDDCCreateFETIDPOperators
1522: .  standard_rhs - the right-hand side for your linear system

1524:    Output Parameters:
1525: -  fetidp_flux_rhs   - the right-hand side for the FETIDP linear system

1527:    Level: developer

1529:    Notes:

1531: .seealso: PCBDDC,PCBDDCCreateFETIDPOperators
1532: @*/
1533: PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
1534: {
1535:   FETIDPMat_ctx  mat_ctx;

1539:   MatShellGetContext(fetidp_mat,&mat_ctx);
1540:   PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));
1541:   return(0);
1542: }
1543: /* -------------------------------------------------------------------------- */

1547: static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
1548: {
1549:   FETIDPMat_ctx  mat_ctx;
1550:   PC_IS*         pcis;
1551:   PC_BDDC*       pcbddc;

1555:   MatShellGetContext(fetidp_mat,&mat_ctx);
1556:   pcis = (PC_IS*)mat_ctx->pc->data;
1557:   pcbddc = (PC_BDDC*)mat_ctx->pc->data;

1559:   /* apply B_delta^T */
1560:   VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
1561:   VecScatterEnd  (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
1562:   MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);
1563:   /* compute rhs for BDDC application */
1564:   VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);
1565:   if (pcbddc->switch_static) {
1566:     VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);
1567:   }
1568:   /* apply BDDC */
1569:   PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);
1570:   /* put values into standard global vector */
1571:   VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
1572:   VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
1573:   if (!pcbddc->switch_static) {
1574:     /* compute values into the interior if solved for the partially subassembled Schur complement */
1575:     MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);
1576:     VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);
1577:     KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);
1578:   }
1579:   VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
1580:   VecScatterEnd  (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
1581:   /* final change of basis if needed
1582:      Is also sums the dirichlet part removed during RHS assembling */
1583:   PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);
1584:   return(0);
1585: }

1589: /*@
1590:  PCBDDCMatFETIDPGetSolution - Compute the physical solution from the solution of the FETIDP linear system

1592:    Collective

1594:    Input Parameters:
1595: +  fetidp_mat        - the FETIDP matrix obtained by calling PCBDDCCreateFETIDPOperators
1596: .  fetidp_flux_sol - the solution of the FETIDP linear system

1598:    Output Parameters:
1599: -  standard_sol      - the solution defined on the physical domain

1601:    Level: developer

1603:    Notes:

1605: .seealso: PCBDDC,PCBDDCCreateFETIDPOperators
1606: @*/
1607: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
1608: {
1609:   FETIDPMat_ctx  mat_ctx;

1613:   MatShellGetContext(fetidp_mat,&mat_ctx);
1614:   PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));
1615:   return(0);
1616: }
1617: /* -------------------------------------------------------------------------- */

1619: extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec);
1620: extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec);
1621: extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat);
1622: extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec);
1623: extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec);
1624: extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC);

1628: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
1629: {

1631:   FETIDPMat_ctx  fetidpmat_ctx;
1632:   Mat            newmat;
1633:   FETIDPPC_ctx   fetidppc_ctx;
1634:   PC             newpc;
1635:   MPI_Comm       comm;

1639:   PetscObjectGetComm((PetscObject)pc,&comm);
1640:   /* FETIDP linear matrix */
1641:   PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);
1642:   PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);
1643:   MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);
1644:   MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);
1645:   MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);
1646:   MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);
1647:   MatSetUp(newmat);
1648:   /* FETIDP preconditioner */
1649:   PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);
1650:   PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);
1651:   PCCreate(comm,&newpc);
1652:   PCSetType(newpc,PCSHELL);
1653:   PCShellSetContext(newpc,fetidppc_ctx);
1654:   PCShellSetApply(newpc,FETIDPPCApply);
1655:   PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);
1656:   PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);
1657:   PCSetOperators(newpc,newmat,newmat);
1658:   PCSetUp(newpc);
1659:   /* return pointers for objects created */
1660:   *fetidp_mat=newmat;
1661:   *fetidp_pc=newpc;
1662:   return(0);
1663: }

1667: /*@
1668:  PCBDDCCreateFETIDPOperators - Create operators for FETIDP

1670:    Collective

1672:    Input Parameters:
1673: +  pc - the BDDC preconditioning context already setup

1675:    Output Parameters:
1676: .  fetidp_mat - shell FETIDP matrix object
1677: .  fetidp_pc  - shell Dirichlet preconditioner for FETIDP matrix

1679:    Options Database Keys:
1680: -    -fetidp_fullyredundant: use or not a fully redundant set of Lagrange multipliers

1682:    Level: developer

1684:    Notes:
1685:      Currently the only operation provided for FETIDP matrix is MatMult

1687: .seealso: PCBDDC
1688: @*/
1689: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
1690: {

1695:   if (pc->setupcalled) {
1696:     PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));
1697:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n");
1698:   return(0);
1699: }
1700: /* -------------------------------------------------------------------------- */
1701: /*MC
1702:    PCBDDC - Balancing Domain Decomposition by Constraints.

1704:    An implementation of the BDDC preconditioner based on

1706: .vb
1707:    [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007
1708:    [2] A. Klawonn and O. B. Widlund. "Dual-Primal FETI Methods for Linear Elasticity", http://cs.nyu.edu/csweb/Research/TechReports/TR2004-855/TR2004-855.pdf
1709:    [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977
1710: .ve

1712:    The matrix to be preconditioned (Pmat) must be of type MATIS.

1714:    Currently works with MATIS matrices with local Neumann matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers.

1716:    It also works with unsymmetric and indefinite problems.

1718:    Unlike 'conventional' interface preconditioners, PCBDDC iterates over all degrees of freedom, not just those on the interface. This allows the use of approximate solvers on the subdomains.

1720:    Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace

1722:    Boundary nodes are split in vertices, edges and faces using information from the local to global mapping of dofs and the local connectivity graph of nodes. The latter can be customized by using PCBDDCSetLocalAdjacencyGraph()

1724:    Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace().

1726:    Change of basis is performed similarly to [2] when requested. When more the one constraint is present on a single connected component (i.e. an edge or a face), a robust method based on local QR factorizations is used.

1728:    The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using MatPartitioning object.

1730:    Options Database Keys:

1732: .    -pc_bddc_use_vertices <1> - use or not vertices in primal space
1733: .    -pc_bddc_use_edges <1> - use or not edges in primal space
1734: .    -pc_bddc_use_faces <0> - use or not faces in primal space
1735: .    -pc_bddc_use_change_of_basis <0> - use change of basis approach (on edges only)
1736: .    -pc_bddc_use_change_on_faces <0> - use change of basis approach on faces if change of basis has been requested
1737: .    -pc_bddc_switch_static <0> - switches from M_2 to M_3 operator (see reference article [1])
1738: .    -pc_bddc_levels <0> - maximum number of levels for multilevel
1739: .    -pc_bddc_coarsening_ratio - H/h ratio at the coarser level
1740: -    -pc_bddc_check_level <0> - set verbosity level of debugging output

1742:    Options for Dirichlet, Neumann or coarse solver can be set with
1743: .vb
1744:       -pc_bddc_dirichlet_
1745:       -pc_bddc_neumann_
1746:       -pc_bddc_coarse_
1747: .ve
1748:    e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg

1750:    When using a multilevel approach, solvers' options at the N-th level can be specified as
1751: .vb
1752:       -pc_bddc_dirichlet_lN_
1753:       -pc_bddc_neumann_lN_
1754:       -pc_bddc_coarse_lN_
1755: .ve
1756:    Note that level number ranges from the finest 0 to the coarsest N.

1758:    Level: intermediate

1760:    Developer notes:

1762:    New deluxe scaling operator will be available soon.

1764:    Contributed by Stefano Zampini

1766: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,  MATIS
1767: M*/

1771: PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
1772: {
1773:   PetscErrorCode      ierr;
1774:   PC_BDDC             *pcbddc;

1777:   /* Creates the private data structure for this preconditioner and attach it to the PC object. */
1778:   PetscNewLog(pc,&pcbddc);
1779:   pc->data  = (void*)pcbddc;

1781:   /* create PCIS data structure */
1782:   PCISCreate(pc);

1784:   /* BDDC customization */
1785:   pcbddc->use_local_adj       = PETSC_TRUE;
1786:   pcbddc->use_vertices        = PETSC_TRUE;
1787:   pcbddc->use_edges           = PETSC_TRUE;
1788:   pcbddc->use_faces           = PETSC_FALSE;
1789:   pcbddc->use_change_of_basis = PETSC_FALSE;
1790:   pcbddc->use_change_on_faces = PETSC_FALSE;
1791:   pcbddc->switch_static       = PETSC_FALSE;
1792:   pcbddc->use_nnsp_true       = PETSC_FALSE; /* not yet exposed */
1793:   pcbddc->dbg_flag            = 0;
1794:   /* private */
1795:   pcbddc->issym                      = PETSC_FALSE;
1796:   pcbddc->BtoNmap                    = 0;
1797:   pcbddc->local_primal_size          = 0;
1798:   pcbddc->n_vertices                 = 0;
1799:   pcbddc->n_actual_vertices          = 0;
1800:   pcbddc->n_constraints              = 0;
1801:   pcbddc->primal_indices_local_idxs  = 0;
1802:   pcbddc->recompute_topography       = PETSC_FALSE;
1803:   pcbddc->coarse_size                = -1;
1804:   pcbddc->new_primal_space           = PETSC_FALSE;
1805:   pcbddc->new_primal_space_local     = PETSC_FALSE;
1806:   pcbddc->global_primal_indices      = 0;
1807:   pcbddc->onearnullspace             = 0;
1808:   pcbddc->onearnullvecs_state        = 0;
1809:   pcbddc->user_primal_vertices       = 0;
1810:   pcbddc->NullSpace                  = 0;
1811:   pcbddc->temp_solution              = 0;
1812:   pcbddc->original_rhs               = 0;
1813:   pcbddc->local_mat                  = 0;
1814:   pcbddc->ChangeOfBasisMatrix        = 0;
1815:   pcbddc->user_ChangeOfBasisMatrix   = 0;
1816:   pcbddc->new_global_mat             = 0;
1817:   pcbddc->coarse_vec                 = 0;
1818:   pcbddc->coarse_rhs                 = 0;
1819:   pcbddc->coarse_ksp                 = 0;
1820:   pcbddc->coarse_phi_B               = 0;
1821:   pcbddc->coarse_phi_D               = 0;
1822:   pcbddc->coarse_psi_B               = 0;
1823:   pcbddc->coarse_psi_D               = 0;
1824:   pcbddc->vec1_P                     = 0;
1825:   pcbddc->vec1_R                     = 0;
1826:   pcbddc->vec2_R                     = 0;
1827:   pcbddc->local_auxmat1              = 0;
1828:   pcbddc->local_auxmat2              = 0;
1829:   pcbddc->R_to_B                     = 0;
1830:   pcbddc->R_to_D                     = 0;
1831:   pcbddc->ksp_D                      = 0;
1832:   pcbddc->ksp_R                      = 0;
1833:   pcbddc->NeumannBoundaries          = 0;
1834:   pcbddc->NeumannBoundariesLocal     = 0;
1835:   pcbddc->DirichletBoundaries        = 0;
1836:   pcbddc->DirichletBoundariesLocal   = 0;
1837:   pcbddc->user_provided_isfordofs    = PETSC_FALSE;
1838:   pcbddc->n_ISForDofs                = 0;
1839:   pcbddc->n_ISForDofsLocal           = 0;
1840:   pcbddc->ISForDofs                  = 0;
1841:   pcbddc->ISForDofsLocal             = 0;
1842:   pcbddc->ConstraintMatrix           = 0;
1843:   pcbddc->use_exact_dirichlet_trick  = PETSC_TRUE;
1844:   pcbddc->coarse_loc_to_glob         = 0;
1845:   pcbddc->coarsening_ratio           = 8;
1846:   pcbddc->current_level              = 0;
1847:   pcbddc->max_levels                 = 0;
1848:   pcbddc->use_coarse_estimates       = PETSC_FALSE;
1849:   pcbddc->coarse_subassembling       = 0;
1850:   pcbddc->coarse_subassembling_init  = 0;

1852:   /* create local graph structure */
1853:   PCBDDCGraphCreate(&pcbddc->mat_graph);

1855:   /* scaling */
1856:   pcbddc->work_scaling          = 0;
1857:   pcbddc->use_deluxe_scaling    = PETSC_FALSE;
1858:   pcbddc->deluxe_threshold      = 1;
1859:   pcbddc->deluxe_rebuild        = PETSC_FALSE;
1860:   pcbddc->deluxe_layers         = -1;
1861:   pcbddc->deluxe_use_useradj    = PETSC_FALSE;
1862:   pcbddc->deluxe_compute_rowadj = PETSC_TRUE;

1864:   /* function pointers */
1865:   pc->ops->apply               = PCApply_BDDC;
1866:   pc->ops->applytranspose      = PCApplyTranspose_BDDC;
1867:   pc->ops->setup               = PCSetUp_BDDC;
1868:   pc->ops->destroy             = PCDestroy_BDDC;
1869:   pc->ops->setfromoptions      = PCSetFromOptions_BDDC;
1870:   pc->ops->view                = 0;
1871:   pc->ops->applyrichardson     = 0;
1872:   pc->ops->applysymmetricleft  = 0;
1873:   pc->ops->applysymmetricright = 0;
1874:   pc->ops->presolve            = PCPreSolve_BDDC;
1875:   pc->ops->postsolve           = PCPostSolve_BDDC;

1877:   /* composing function */
1878:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);
1879:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);
1880:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);
1881:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);
1882:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);
1883:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);
1884:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);
1885:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);
1886:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);
1887:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);
1888:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);
1889:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);
1890:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);
1891:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);
1892:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);
1893:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);
1894:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);
1895:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);
1896:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);
1897:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);
1898:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);
1899:   return(0);
1900: }