Actual source code: bddc.c
petsc-dev 2013-06-18
1: /* TODOLIST
2: DofSplitting and DM attached to pc?
3: Change SetNeumannBoundaries to SetNeumannBoundariesLocal and provide new SetNeumannBoundaries (same Dirichlet)
4: change how to deal with the coarse problem (PCBDDCSetCoarseEnvironment):
5: - simplify coarse problem structure -> PCBDDC or PCREDUDANT, nothing else -> same comm for all levels?
6: - remove coarse enums and allow use of PCBDDCGetCoarseKSP
7: - remove metis dependency -> use MatPartitioning for multilevel -> Assemble serial adjacency in PCBDDCAnalyzeInterface?
8: code refactoring:
9: - pick up better names for static functions
10: change options structure:
11: - insert BDDC into MG framework?
12: provide other ops? Ask to developers
13: remove all unused printf
14: man pages
15: */
17: /* ----------------------------------------------------------------------------------------------------------------------------------------------
18: Implementation of BDDC preconditioner based on:
19: C. Dohrmann "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007
20: ---------------------------------------------------------------------------------------------------------------------------------------------- */
22: #include "bddc.h" /*I "petscpc.h" I*/ /* includes for fortran wrappers */
23: #include bddcprivate.h
24: #include <petscblaslapack.h>
26: /* prototypes for static functions contained in bddc.c */
27: static PetscErrorCode PCBDDCSetUseExactDirichlet(PC,PetscBool);
28: static PetscErrorCode PCBDDCSetLevel(PC,PetscInt);
29: static PetscErrorCode PCBDDCCoarseSetUp(PC);
30: static PetscErrorCode PCBDDCSetUpCoarseEnvironment(PC,PetscScalar*);
32: /* -------------------------------------------------------------------------- */
35: PetscErrorCode PCSetFromOptions_BDDC(PC pc)
36: {
37: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
41: PetscOptionsHead("BDDC options");
42: /* Verbose debugging of main data structures */
43: PetscOptionsInt("-pc_bddc_check_level" ,"Verbose (debugging) output for PCBDDC" ,"none",pcbddc->dbg_flag ,&pcbddc->dbg_flag ,NULL);
44: /* Some customization for default primal space */
45: PetscOptionsBool("-pc_bddc_vertices_only" ,"Use only vertices in coarse space (i.e. discard constraints)","none",pcbddc->vertices_flag ,&pcbddc->vertices_flag ,NULL);
46: PetscOptionsBool("-pc_bddc_constraints_only","Use only constraints in coarse space (i.e. discard vertices)","none",pcbddc->constraints_flag,&pcbddc->constraints_flag,NULL);
47: PetscOptionsBool("-pc_bddc_faces_only" ,"Use only faces among constraints of coarse space (i.e. discard edges)" ,"none",pcbddc->faces_flag ,&pcbddc->faces_flag ,NULL);
48: PetscOptionsBool("-pc_bddc_edges_only" ,"Use only edges among constraints of coarse space (i.e. discard faces)" ,"none",pcbddc->edges_flag ,&pcbddc->edges_flag ,NULL);
49: /* Coarse solver context */
50: static const char * const avail_coarse_problems[] = {"sequential","replicated","parallel","multilevel","CoarseProblemType","PC_BDDC_",0}; /*order of choiches depends on ENUM defined in bddc.h */
51: PetscOptionsEnum("-pc_bddc_coarse_problem_type","Set coarse problem type","none",avail_coarse_problems,(PetscEnum)pcbddc->coarse_problem_type,(PetscEnum*)&pcbddc->coarse_problem_type,NULL);
52: /* Two different application of BDDC to the whole set of dofs, internal and interface */
53: PetscOptionsBool("-pc_bddc_switch_preconditioning_type","Switch between M_2 (default) and M_3 preconditioners (as defined by Dohrmann)","none",pcbddc->inexact_prec_type,&pcbddc->inexact_prec_type,NULL);
54: PetscOptionsBool("-pc_bddc_use_change_of_basis","Use change of basis approach for primal space","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);
55: PetscOptionsBool("-pc_bddc_use_change_on_faces","Use change of basis approach for face constraints","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);
56: if (!pcbddc->use_change_of_basis) {
57: pcbddc->use_change_on_faces = PETSC_FALSE;
58: }
59: PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);
60: PetscOptionsInt("-pc_bddc_max_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);
61: PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);
62: PetscOptionsTail();
63: return(0);
64: }
65: /* -------------------------------------------------------------------------- */
68: static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
69: {
70: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
74: ISDestroy(&pcbddc->user_primal_vertices);
75: PetscObjectReference((PetscObject)PrimalVertices);
76: pcbddc->user_primal_vertices = PrimalVertices;
77: return(0);
78: }
81: /*@
82: PCBDDCSetPrimalVerticesLocalIS - Set user defined primal vertices in PCBDDC.
84: Not collective
86: Input Parameters:
87: + pc - the preconditioning context
88: - PrimalVertices - index sets of primal vertices in local numbering
90: Level: intermediate
92: Notes:
94: .seealso: PCBDDC
95: @*/
96: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
97: {
103: PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));
104: return(0);
105: }
106: /* -------------------------------------------------------------------------- */
109: static PetscErrorCode PCBDDCSetCoarseProblemType_BDDC(PC pc, CoarseProblemType CPT)
110: {
111: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
114: pcbddc->coarse_problem_type = CPT;
115: return(0);
116: }
120: /*@
121: PCBDDCSetCoarseProblemType - Set coarse problem type in PCBDDC.
123: Not collective
125: Input Parameters:
126: + pc - the preconditioning context
127: - CoarseProblemType - pick a better name and explain what this is
129: Level: intermediate
131: Notes:
132: Not collective but all procs must call with same arguments.
134: .seealso: PCBDDC
135: @*/
136: PetscErrorCode PCBDDCSetCoarseProblemType(PC pc, CoarseProblemType CPT)
137: {
142: PetscTryMethod(pc,"PCBDDCSetCoarseProblemType_C",(PC,CoarseProblemType),(pc,CPT));
143: return(0);
144: }
145: /* -------------------------------------------------------------------------- */
148: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k)
149: {
150: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
153: pcbddc->coarsening_ratio=k;
154: return(0);
155: }
159: /*@
160: PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel coarsening
162: Logically collective on PC
164: Input Parameters:
165: + pc - the preconditioning context
166: - k - coarsening ratio
168: Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level.
170: Level: intermediate
172: Notes:
174: .seealso: PCBDDC
175: @*/
176: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k)
177: {
182: PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));
183: return(0);
184: }
185: /* -------------------------------------------------------------------------- */
189: static PetscErrorCode PCBDDCSetMaxLevels_BDDC(PC pc,PetscInt max_levels)
190: {
191: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
194: pcbddc->max_levels=max_levels;
195: return(0);
196: }
200: /*@
201: PCBDDCSetMaxLevels - Sets the maximum number of levels within the multilevel approach.
203: Logically collective on PC
205: Input Parameters:
206: + pc - the preconditioning context
207: - max_levels - the maximum number of levels
209: Default value is 1, i.e. coarse problem will be solved inexactly with one application
210: of PCBDDC preconditioner if the multilevel approach is requested.
212: Level: intermediate
214: Notes:
216: .seealso: PCBDDC
217: @*/
218: PetscErrorCode PCBDDCSetMaxLevels(PC pc,PetscInt max_levels)
219: {
224: PetscTryMethod(pc,"PCBDDCSetMaxLevels_C",(PC,PetscInt),(pc,max_levels));
225: return(0);
226: }
227: /* -------------------------------------------------------------------------- */
231: static PetscErrorCode PCBDDCSetNullSpace_BDDC(PC pc,MatNullSpace NullSpace)
232: {
233: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
237: PetscObjectReference((PetscObject)NullSpace);
238: MatNullSpaceDestroy(&pcbddc->NullSpace);
239: pcbddc->NullSpace=NullSpace;
240: return(0);
241: }
245: /*@
246: PCBDDCSetNullSpace - Set NullSpace of global operator of BDDC preconditioned mat.
248: Logically collective on PC and MatNullSpace
250: Input Parameters:
251: + pc - the preconditioning context
252: - NullSpace - Null space of the linear operator to be preconditioned.
254: Level: intermediate
256: Notes:
258: .seealso: PCBDDC
259: @*/
260: PetscErrorCode PCBDDCSetNullSpace(PC pc,MatNullSpace NullSpace)
261: {
267: PetscTryMethod(pc,"PCBDDCSetNullSpace_C",(PC,MatNullSpace),(pc,NullSpace));
268: return(0);
269: }
270: /* -------------------------------------------------------------------------- */
274: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries)
275: {
276: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
280: ISDestroy(&pcbddc->DirichletBoundaries);
281: PetscObjectReference((PetscObject)DirichletBoundaries);
282: pcbddc->DirichletBoundaries=DirichletBoundaries;
283: return(0);
284: }
288: /*@
289: PCBDDCSetDirichletBoundaries - Set index set defining subdomain part (in local ordering)
290: of Dirichlet boundaries for the global problem.
292: Not collective
294: Input Parameters:
295: + pc - the preconditioning context
296: - DirichletBoundaries - sequential index set defining the subdomain part of Dirichlet boundaries (can be NULL)
298: Level: intermediate
300: Notes:
302: .seealso: PCBDDC
303: @*/
304: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries)
305: {
311: PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));
312: return(0);
313: }
314: /* -------------------------------------------------------------------------- */
318: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries)
319: {
320: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
324: ISDestroy(&pcbddc->NeumannBoundaries);
325: PetscObjectReference((PetscObject)NeumannBoundaries);
326: pcbddc->NeumannBoundaries=NeumannBoundaries;
327: return(0);
328: }
332: /*@
333: PCBDDCSetNeumannBoundaries - Set index set defining subdomain part (in local ordering)
334: of Neumann boundaries for the global problem.
336: Not collective
338: Input Parameters:
339: + pc - the preconditioning context
340: - NeumannBoundaries - sequential index set defining the subdomain part of Neumann boundaries (can be NULL)
342: Level: intermediate
344: Notes:
346: .seealso: PCBDDC
347: @*/
348: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries)
349: {
355: PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));
356: return(0);
357: }
358: /* -------------------------------------------------------------------------- */
362: static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries)
363: {
364: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
367: *DirichletBoundaries = pcbddc->DirichletBoundaries;
368: return(0);
369: }
373: /*@
374: PCBDDCGetDirichletBoundaries - Get index set defining subdomain part (in local ordering)
375: of Dirichlet boundaries for the global problem.
377: Not collective
379: Input Parameters:
380: + pc - the preconditioning context
382: Output Parameters:
383: + DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries
385: Level: intermediate
387: Notes:
389: .seealso: PCBDDC
390: @*/
391: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries)
392: {
397: PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));
398: return(0);
399: }
400: /* -------------------------------------------------------------------------- */
404: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries)
405: {
406: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
409: *NeumannBoundaries = pcbddc->NeumannBoundaries;
410: return(0);
411: }
415: /*@
416: PCBDDCGetNeumannBoundaries - Get index set defining subdomain part (in local ordering)
417: of Neumann boundaries for the global problem.
419: Not collective
421: Input Parameters:
422: + pc - the preconditioning context
424: Output Parameters:
425: + NeumannBoundaries - index set defining the subdomain part of Neumann boundaries
427: Level: intermediate
429: Notes:
431: .seealso: PCBDDC
432: @*/
433: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries)
434: {
439: PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));
440: return(0);
441: }
442: /* -------------------------------------------------------------------------- */
446: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
447: {
448: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
449: PCBDDCGraph mat_graph = pcbddc->mat_graph;
453: /* free old CSR */
454: PCBDDCGraphResetCSR(mat_graph);
455: /* get CSR into graph structure */
456: if (copymode == PETSC_COPY_VALUES) {
457: PetscMalloc((nvtxs+1)*sizeof(PetscInt),&mat_graph->xadj);
458: PetscMalloc(xadj[nvtxs]*sizeof(PetscInt),&mat_graph->adjncy);
459: PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));
460: PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));
461: } else if (copymode == PETSC_OWN_POINTER) {
462: mat_graph->xadj = (PetscInt*)xadj;
463: mat_graph->adjncy = (PetscInt*)adjncy;
464: }
465: mat_graph->nvtxs_csr = nvtxs;
466: return(0);
467: }
471: /*@
472: PCBDDCSetLocalAdjacencyGraph - Set CSR graph of local matrix for use of PCBDDC.
474: Not collective
476: Input Parameters:
477: + pc - the preconditioning context
478: - nvtxs - number of local vertices of the graph
479: - xadj, adjncy - the CSR graph
480: - copymode - either PETSC_COPY_VALUES or PETSC_OWN_POINTER. In the former case the user must free the array passed in;
481: in the latter case, memory must be obtained with PetscMalloc.
483: Level: intermediate
485: Notes:
487: .seealso: PCBDDC
488: @*/
489: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
490: {
491: void (*f)(void) = 0;
498: if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) {
499: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d in %s\n",copymode,__FUNCT__);
500: }
501: PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));
502: /* free arrays if PCBDDC is not the PC type */
503: PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);
504: if (!f && copymode == PETSC_OWN_POINTER) {
505: PetscFree(xadj);
506: PetscFree(adjncy);
507: }
508: return(0);
509: }
510: /* -------------------------------------------------------------------------- */
514: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[])
515: {
516: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
517: PetscInt i;
521: /* Destroy ISes if they were already set */
522: for (i=0;i<pcbddc->n_ISForDofs;i++) {
523: ISDestroy(&pcbddc->ISForDofs[i]);
524: }
525: PetscFree(pcbddc->ISForDofs);
526: /* allocate space then set */
527: PetscMalloc(n_is*sizeof(IS),&pcbddc->ISForDofs);
528: for (i=0;i<n_is;i++) {
529: PetscObjectReference((PetscObject)ISForDofs[i]);
530: pcbddc->ISForDofs[i]=ISForDofs[i];
531: }
532: pcbddc->n_ISForDofs=n_is;
533: return(0);
534: }
538: /*@
539: PCBDDCSetDofsSplitting - Set index sets defining fields of local mat.
541: Not collective
543: Input Parameters:
544: + pc - the preconditioning context
545: - n - number of index sets defining the fields
546: - IS[] - array of IS describing the fields
548: Level: intermediate
550: Notes:
552: .seealso: PCBDDC
553: @*/
554: PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[])
555: {
560: PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));
561: return(0);
562: }
563: /* -------------------------------------------------------------------------- */
566: /* -------------------------------------------------------------------------- */
567: /*
568: PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial
569: guess if a transformation of basis approach has been selected.
571: Input Parameter:
572: + pc - the preconditioner contex
574: Application Interface Routine: PCPreSolve()
576: Notes:
577: The interface routine PCPreSolve() is not usually called directly by
578: the user, but instead is called by KSPSolve().
579: */
580: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
581: {
583: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
584: PC_IS *pcis = (PC_IS*)(pc->data);
585: Mat_IS *matis = (Mat_IS*)pc->pmat->data;
586: Mat temp_mat;
587: IS dirIS;
588: PetscInt dirsize,i,*is_indices;
589: PetscScalar *array_x,*array_diagonal;
590: Vec used_vec;
591: PetscBool guess_nonzero;
594: if (x) {
595: PetscObjectReference((PetscObject)x);
596: used_vec = x;
597: } else {
598: PetscObjectReference((PetscObject)pcbddc->temp_solution);
599: used_vec = pcbddc->temp_solution;
600: VecSet(used_vec,0.0);
601: }
602: /* hack into ksp data structure PCPreSolve comes earlier in src/ksp/ksp/interface/itfunc.c */
603: if (ksp) {
604: KSPGetInitialGuessNonzero(ksp,&guess_nonzero);
605: if ( !guess_nonzero ) {
606: VecSet(used_vec,0.0);
607: }
608: }
609: /* store the original rhs */
610: VecCopy(rhs,pcbddc->original_rhs);
612: /* Take into account zeroed rows -> change rhs and store solution removed */
613: MatGetDiagonal(pc->pmat,pcis->vec1_global);
614: VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);
615: VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);
616: VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);
617: VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
618: VecScatterEnd(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
619: PCBDDCGetDirichletBoundaries(pc,&dirIS);
620: if (dirIS) {
621: ISGetSize(dirIS,&dirsize);
622: VecGetArray(pcis->vec1_N,&array_x);
623: VecGetArray(pcis->vec2_N,&array_diagonal);
624: ISGetIndices(dirIS,(const PetscInt**)&is_indices);
625: for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
626: ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);
627: VecRestoreArray(pcis->vec2_N,&array_diagonal);
628: VecRestoreArray(pcis->vec1_N,&array_x);
629: }
630: VecScatterBegin(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);
631: VecScatterEnd(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);
633: /* remove the computed solution from the rhs */
634: VecScale(used_vec,-1.0);
635: MatMultAdd(pc->pmat,used_vec,rhs,rhs);
636: VecScale(used_vec,-1.0);
638: /* store partially computed solution and set initial guess */
639: if (x) {
640: VecCopy(used_vec,pcbddc->temp_solution);
641: VecSet(used_vec,0.0);
642: if (pcbddc->use_exact_dirichlet && !pcbddc->coarse_psi_B) {
643: VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
644: VecScatterEnd (pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
645: KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);
646: VecScatterBegin(pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);
647: VecScatterEnd (pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);
648: if (ksp) {
649: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
650: }
651: }
652: }
654: /* rhs change of basis */
655: if (pcbddc->use_change_of_basis) {
656: /* swap pointers for local matrices */
657: temp_mat = matis->A;
658: matis->A = pcbddc->local_mat;
659: pcbddc->local_mat = temp_mat;
660: /* Get local rhs and apply transformation of basis */
661: VecScatterBegin(pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
662: VecScatterEnd (pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
663: /* from original basis to modified basis */
664: MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);
665: /* put back modified values into the global vec using INSERT_VALUES copy mode */
666: VecScatterBegin(pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);
667: VecScatterEnd (pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);
668: }
669: if (ksp && pcbddc->NullSpace) {
670: MatNullSpaceRemove(pcbddc->NullSpace,used_vec);
671: MatNullSpaceRemove(pcbddc->NullSpace,rhs);
672: }
673: VecDestroy(&used_vec);
674: return(0);
675: }
676: /* -------------------------------------------------------------------------- */
679: /* -------------------------------------------------------------------------- */
680: /*
681: PCPostSolve_BDDC - Changes the computed solution if a transformation of basis
682: approach has been selected. Also, restores rhs to its original state.
684: Input Parameter:
685: + pc - the preconditioner contex
687: Application Interface Routine: PCPostSolve()
689: Notes:
690: The interface routine PCPostSolve() is not usually called directly by
691: the user, but instead is called by KSPSolve().
692: */
693: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
694: {
696: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
697: PC_IS *pcis = (PC_IS*)(pc->data);
698: Mat_IS *matis = (Mat_IS*)pc->pmat->data;
699: Mat temp_mat;
702: if (pcbddc->use_change_of_basis) {
703: /* swap pointers for local matrices */
704: temp_mat = matis->A;
705: matis->A = pcbddc->local_mat;
706: pcbddc->local_mat = temp_mat;
707: /* restore rhs to its original state */
708: if (rhs) {
709: VecCopy(pcbddc->original_rhs,rhs);
710: }
711: /* Get Local boundary and apply transformation of basis to solution vector */
712: VecScatterBegin(pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
713: VecScatterEnd (pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
714: /* from modified basis to original basis */
715: MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);
716: /* put back modified values into the global vec using INSERT_VALUES copy mode */
717: VecScatterBegin(pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);
718: VecScatterEnd (pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);
719: }
720: /* add solution removed in presolve */
721: if (x) {
722: VecAXPY(x,1.0,pcbddc->temp_solution);
723: }
724: return(0);
725: }
726: /* -------------------------------------------------------------------------- */
729: /* -------------------------------------------------------------------------- */
730: /*
731: PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner
732: by setting data structures and options.
734: Input Parameter:
735: + pc - the preconditioner context
737: Application Interface Routine: PCSetUp()
739: Notes:
740: The interface routine PCSetUp() is not usually called directly by
741: the user, but instead is called by PCApply() if necessary.
742: */
743: PetscErrorCode PCSetUp_BDDC(PC pc)
744: {
746: PC_BDDC* pcbddc = (PC_BDDC*)pc->data;
747: MatStructure flag;
748: PetscBool computeis,computetopography,computesolvers;
751: /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other nonoverlapping preconditioners */
752: /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
753: So, we set to pcnone the Neumann problem of pcis in order to avoid unneeded computation
754: Also, we decide to directly build the (same) Dirichlet problem */
755: PetscOptionsSetValue("-is_localN_pc_type","none");
756: PetscOptionsSetValue("-is_localD_pc_type","none");
757: /* Get stdout for dbg */
758: if (pcbddc->dbg_flag && !pcbddc->dbg_viewer) {
759: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&pcbddc->dbg_viewer);
760: PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);
761: }
762: /* first attempt to split work */
763: if (pc->setupcalled) {
764: computeis = PETSC_FALSE;
765: PCGetOperators(pc,NULL,NULL,&flag);
766: if (flag == SAME_PRECONDITIONER) {
767: computetopography = PETSC_FALSE;
768: computesolvers = PETSC_FALSE;
769: } else if (flag == SAME_NONZERO_PATTERN) {
770: computetopography = PETSC_FALSE;
771: computesolvers = PETSC_TRUE;
772: } else { /* DIFFERENT_NONZERO_PATTERN */
773: computetopography = PETSC_TRUE;
774: computesolvers = PETSC_TRUE;
775: }
776: } else {
777: computeis = PETSC_TRUE;
778: computetopography = PETSC_TRUE;
779: computesolvers = PETSC_TRUE;
780: }
781: /* Set up all the "iterative substructuring" common block */
782: if (computeis) {
783: PCISSetUp(pc);
784: }
785: /* Analyze interface and set up local constraint and change of basis matrices */
786: if (computetopography) {
787: /* reset data */
788: PCBDDCResetTopography(pc);
789: PCBDDCAnalyzeInterface(pc);
790: PCBDDCConstraintsSetUp(pc);
791: }
792: if (computesolvers) {
793: /* reset data */
794: PCBDDCResetSolvers(pc);
795: PCBDDCScalingDestroy(pc);
796: /* Create coarse and local stuffs used for evaluating action of preconditioner */
797: PCBDDCCoarseSetUp(pc);
798: PCBDDCScalingSetUp(pc);
799: }
800: return(0);
801: }
803: /* -------------------------------------------------------------------------- */
804: /*
805: PCApply_BDDC - Applies the BDDC preconditioner to a vector.
807: Input Parameters:
808: . pc - the preconditioner context
809: . r - input vector (global)
811: Output Parameter:
812: . z - output vector (global)
814: Application Interface Routine: PCApply()
815: */
818: PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z)
819: {
820: PC_IS *pcis = (PC_IS*)(pc->data);
821: PC_BDDC *pcbddc = (PC_BDDC*)(pc->data);
822: PetscErrorCode ierr;
823: const PetscScalar one = 1.0;
824: const PetscScalar m_one = -1.0;
825: const PetscScalar zero = 0.0;
827: /* This code is similar to that provided in nn.c for PCNN
828: NN interface preconditioner changed to BDDC
829: Added support for M_3 preconditioner in the reference article (code is active if pcbddc->inexact_prec_type = PETSC_TRUE) */
832: if (!pcbddc->use_exact_dirichlet || pcbddc->coarse_psi_B) {
833: /* First Dirichlet solve */
834: VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
835: VecScatterEnd (pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
836: KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);
837: /*
838: Assembling right hand side for BDDC operator
839: - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE)
840: - pcis->vec1_B the interface part of the global vector z
841: */
842: VecScale(pcis->vec2_D,m_one);
843: MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);
844: if (pcbddc->inexact_prec_type) { MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D); }
845: VecScale(pcis->vec2_D,m_one);
846: VecCopy(r,z);
847: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);
848: VecScatterEnd (pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);
849: PCBDDCScalingRestriction(pc,z,pcis->vec1_B);
850: } else {
851: VecSet(pcis->vec1_D,zero);
852: VecSet(pcis->vec2_D,zero);
853: PCBDDCScalingRestriction(pc,r,pcis->vec1_B);
854: }
856: /* Apply interface preconditioner
857: input/output vecs: pcis->vec1_B and pcis->vec1_D */
858: PCBDDCApplyInterfacePreconditioner(pc);
860: /* Apply transpose of partition of unity operator */
861: PCBDDCScalingExtension(pc,pcis->vec1_B,z);
863: /* Second Dirichlet solve and assembling of output */
864: VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
865: VecScatterEnd (pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
866: MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);
867: if (pcbddc->inexact_prec_type) { MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D); }
868: KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcbddc->vec4_D);
869: VecScale(pcbddc->vec4_D,m_one);
870: if (pcbddc->inexact_prec_type) { VecAXPY (pcbddc->vec4_D,one,pcis->vec1_D); }
871: VecAXPY (pcis->vec2_D,one,pcbddc->vec4_D);
872: VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
873: VecScatterEnd (pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
874: return(0);
875: }
876: /* -------------------------------------------------------------------------- */
880: PetscErrorCode PCDestroy_BDDC(PC pc)
881: {
882: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
886: /* free data created by PCIS */
887: PCISDestroy(pc);
888: /* free BDDC custom data */
889: PCBDDCResetCustomization(pc);
890: /* destroy objects related to topography */
891: PCBDDCResetTopography(pc);
892: /* free allocated graph structure */
893: PetscFree(pcbddc->mat_graph);
894: /* free data for scaling operator */
895: PCBDDCScalingDestroy(pc);
896: /* free solvers stuff */
897: PCBDDCResetSolvers(pc);
898: /* remove functions */
899: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);
900: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);
901: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetMaxLevels_C",NULL);
902: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);
903: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);
904: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);
905: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);
906: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);
907: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseProblemType_C",NULL);
908: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);
909: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);
910: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);
911: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);
912: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);
913: /* Free the private data structure */
914: PetscFree(pc->data);
915: return(0);
916: }
917: /* -------------------------------------------------------------------------- */
921: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
922: {
923: FETIDPMat_ctx mat_ctx;
924: PC_IS* pcis;
925: PC_BDDC* pcbddc;
929: MatShellGetContext(fetidp_mat,&mat_ctx);
930: pcis = (PC_IS*)mat_ctx->pc->data;
931: pcbddc = (PC_BDDC*)mat_ctx->pc->data;
933: /* change of basis for physical rhs if needed
934: It also changes the rhs in case of dirichlet boundaries */
935: (*mat_ctx->pc->ops->presolve)(mat_ctx->pc,NULL,standard_rhs,NULL);
936: /* store vectors for computation of fetidp final solution */
937: VecScatterBegin(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);
938: VecScatterEnd(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);
939: /* scale rhs since it should be unassembled : TODO use counter scaling? (also below) */
940: VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
941: VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
942: /* Apply partition of unity */
943: VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);
944: /* PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B); */
945: if (!pcbddc->inexact_prec_type) {
946: /* compute partially subassembled Schur complement right-hand side */
947: KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);
948: MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);
949: VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);
950: VecSet(standard_rhs,0.0);
951: VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);
952: VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);
953: /* PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B); */
954: VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
955: VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
956: VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);
957: }
958: /* BDDC rhs */
959: VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);
960: if (pcbddc->inexact_prec_type) {
961: VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);
962: }
963: /* apply BDDC */
964: PCBDDCApplyInterfacePreconditioner(mat_ctx->pc);
965: /* Application of B_delta and assembling of rhs for fetidp fluxes */
966: VecSet(fetidp_flux_rhs,0.0);
967: MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);
968: VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);
969: VecScatterEnd (mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);
970: /* restore original rhs */
971: VecCopy(pcbddc->original_rhs,standard_rhs);
972: return(0);
973: }
977: /*@
978: PCBDDCMatFETIDPGetRHS - Get rhs for FETIDP linear system.
980: Collective
982: Input Parameters:
983: + fetidp_mat - the FETIDP mat obtained by a call to PCBDDCCreateFETIDPOperators
984: + standard_rhs - the rhs of your linear system
986: Output Parameters:
987: + fetidp_flux_rhs - the rhs of the FETIDP linear system
989: Level: developer
991: Notes:
993: .seealso: PCBDDC
994: @*/
995: PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
996: {
997: FETIDPMat_ctx mat_ctx;
1001: MatShellGetContext(fetidp_mat,&mat_ctx);
1002: PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));
1003: return(0);
1004: }
1005: /* -------------------------------------------------------------------------- */
1009: static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
1010: {
1011: FETIDPMat_ctx mat_ctx;
1012: PC_IS* pcis;
1013: PC_BDDC* pcbddc;
1017: MatShellGetContext(fetidp_mat,&mat_ctx);
1018: pcis = (PC_IS*)mat_ctx->pc->data;
1019: pcbddc = (PC_BDDC*)mat_ctx->pc->data;
1021: /* apply B_delta^T */
1022: VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
1023: VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
1024: MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);
1025: /* compute rhs for BDDC application */
1026: VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);
1027: if (pcbddc->inexact_prec_type) {
1028: VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);
1029: }
1030: /* apply BDDC */
1031: PCBDDCApplyInterfacePreconditioner(mat_ctx->pc);
1032: /* put values into standard global vector */
1033: VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
1034: VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
1035: if (!pcbddc->inexact_prec_type) {
1036: /* compute values into the interior if solved for the partially subassembled Schur complement */
1037: MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);
1038: VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);
1039: KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);
1040: }
1041: VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
1042: VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
1043: /* final change of basis if needed
1044: Is also sums the dirichlet part removed during RHS assembling */
1045: (*mat_ctx->pc->ops->postsolve)(mat_ctx->pc,NULL,NULL,standard_sol);
1046: return(0);
1048: }
1052: /*@
1053: PCBDDCMatFETIDPGetSolution - Get Solution for FETIDP linear system.
1055: Collective
1057: Input Parameters:
1058: + fetidp_mat - the FETIDP mat obtained by a call to PCBDDCCreateFETIDPOperators
1059: + fetidp_flux_sol - the solution of the FETIDP linear system
1061: Output Parameters:
1062: + standard_sol - the solution on the global domain
1064: Level: developer
1066: Notes:
1068: .seealso: PCBDDC
1069: @*/
1070: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
1071: {
1072: FETIDPMat_ctx mat_ctx;
1076: MatShellGetContext(fetidp_mat,&mat_ctx);
1077: PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));
1078: return(0);
1079: }
1080: /* -------------------------------------------------------------------------- */
1082: extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec);
1083: extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat);
1084: extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec);
1085: extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC);
1089: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
1090: {
1092: FETIDPMat_ctx fetidpmat_ctx;
1093: Mat newmat;
1094: FETIDPPC_ctx fetidppc_ctx;
1095: PC newpc;
1096: MPI_Comm comm;
1100: PetscObjectGetComm((PetscObject)pc,&comm);
1101: /* FETIDP linear matrix */
1102: PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);
1103: PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);
1104: MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);
1105: MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);
1106: MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);
1107: MatSetUp(newmat);
1108: /* FETIDP preconditioner */
1109: PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);
1110: PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);
1111: PCCreate(comm,&newpc);
1112: PCSetType(newpc,PCSHELL);
1113: PCShellSetContext(newpc,fetidppc_ctx);
1114: PCShellSetApply(newpc,FETIDPPCApply);
1115: PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);
1116: PCSetOperators(newpc,newmat,newmat,SAME_PRECONDITIONER);
1117: PCSetUp(newpc);
1118: /* return pointers for objects created */
1119: *fetidp_mat=newmat;
1120: *fetidp_pc=newpc;
1121: return(0);
1122: }
1126: /*@
1127: PCBDDCCreateFETIDPOperators - Create operators for FETIDP.
1129: Collective
1131: Input Parameters:
1132: + pc - the BDDC preconditioning context (setup must be already called)
1134: Level: developer
1136: Notes:
1138: .seealso: PCBDDC
1139: @*/
1140: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
1141: {
1146: if (pc->setupcalled) {
1147: PetscTryMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));
1148: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n");
1149: return(0);
1150: }
1151: /* -------------------------------------------------------------------------- */
1152: /*MC
1153: PCBDDC - Balancing Domain Decomposition by Constraints.
1155: Options Database Keys:
1156: . -pcbddc ??? -
1158: Level: intermediate
1160: Notes: The matrix used with this preconditioner must be of type MATIS
1162: Unlike more 'conventional' interface preconditioners, this iterates over ALL the
1163: degrees of freedom, NOT just those on the interface (this allows the use of approximate solvers
1164: on the subdomains).
1166: Options for the coarse grid preconditioner can be set with -
1167: Options for the Dirichlet subproblem can be set with -
1168: Options for the Neumann subproblem can be set with -
1170: Contributed by Stefano Zampini
1172: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS
1173: M*/
1177: PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
1178: {
1179: PetscErrorCode ierr;
1180: PC_BDDC *pcbddc;
1183: /* Creates the private data structure for this preconditioner and attach it to the PC object. */
1184: PetscNewLog(pc,PC_BDDC,&pcbddc);
1185: pc->data = (void*)pcbddc;
1187: /* create PCIS data structure */
1188: PCISCreate(pc);
1190: /* BDDC specific */
1191: pcbddc->user_primal_vertices = 0;
1192: pcbddc->NullSpace = 0;
1193: pcbddc->temp_solution = 0;
1194: pcbddc->original_rhs = 0;
1195: pcbddc->local_mat = 0;
1196: pcbddc->ChangeOfBasisMatrix = 0;
1197: pcbddc->use_change_of_basis = PETSC_TRUE;
1198: pcbddc->use_change_on_faces = PETSC_FALSE;
1199: pcbddc->coarse_vec = 0;
1200: pcbddc->coarse_rhs = 0;
1201: pcbddc->coarse_ksp = 0;
1202: pcbddc->coarse_phi_B = 0;
1203: pcbddc->coarse_phi_D = 0;
1204: pcbddc->coarse_psi_B = 0;
1205: pcbddc->coarse_psi_D = 0;
1206: pcbddc->vec1_P = 0;
1207: pcbddc->vec1_R = 0;
1208: pcbddc->vec2_R = 0;
1209: pcbddc->local_auxmat1 = 0;
1210: pcbddc->local_auxmat2 = 0;
1211: pcbddc->R_to_B = 0;
1212: pcbddc->R_to_D = 0;
1213: pcbddc->ksp_D = 0;
1214: pcbddc->ksp_R = 0;
1215: pcbddc->local_primal_indices = 0;
1216: pcbddc->inexact_prec_type = PETSC_FALSE;
1217: pcbddc->NeumannBoundaries = 0;
1218: pcbddc->ISForDofs = 0;
1219: pcbddc->ConstraintMatrix = 0;
1220: pcbddc->use_nnsp_true = PETSC_FALSE;
1221: pcbddc->local_primal_sizes = 0;
1222: pcbddc->local_primal_displacements = 0;
1223: pcbddc->coarse_loc_to_glob = 0;
1224: pcbddc->dbg_flag = 0;
1225: pcbddc->coarsening_ratio = 8;
1226: pcbddc->use_exact_dirichlet = PETSC_TRUE;
1227: pcbddc->current_level = 0;
1228: pcbddc->max_levels = 1;
1229: pcbddc->replicated_local_primal_indices = 0;
1230: pcbddc->replicated_local_primal_values = 0;
1232: /* create local graph structure */
1233: PCBDDCGraphCreate(&pcbddc->mat_graph);
1235: /* scaling */
1236: pcbddc->use_deluxe_scaling = PETSC_FALSE;
1237: pcbddc->work_scaling = 0;
1239: /* function pointers */
1240: pc->ops->apply = PCApply_BDDC;
1241: pc->ops->applytranspose = 0;
1242: pc->ops->setup = PCSetUp_BDDC;
1243: pc->ops->destroy = PCDestroy_BDDC;
1244: pc->ops->setfromoptions = PCSetFromOptions_BDDC;
1245: pc->ops->view = 0;
1246: pc->ops->applyrichardson = 0;
1247: pc->ops->applysymmetricleft = 0;
1248: pc->ops->applysymmetricright = 0;
1249: pc->ops->presolve = PCPreSolve_BDDC;
1250: pc->ops->postsolve = PCPostSolve_BDDC;
1252: /* composing function */
1253: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);
1254: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);
1255: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetMaxLevels_C",PCBDDCSetMaxLevels_BDDC);
1256: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);
1257: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);
1258: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);
1259: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);
1260: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);
1261: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseProblemType_C",PCBDDCSetCoarseProblemType_BDDC);
1262: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);
1263: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);
1264: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);
1265: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);
1266: PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);
1267: return(0);
1268: }
1270: /* -------------------------------------------------------------------------- */
1271: /* All static functions from now on */
1272: /* -------------------------------------------------------------------------- */
1276: static PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool use)
1277: {
1278: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
1281: pcbddc->use_exact_dirichlet=use;
1282: return(0);
1283: }
1287: static PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level)
1288: {
1289: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
1292: pcbddc->current_level=level;
1293: return(0);
1294: }
1296: /* -------------------------------------------------------------------------- */
1299: static PetscErrorCode PCBDDCCoarseSetUp(PC pc)
1300: {
1301: PetscErrorCode ierr;
1303: PC_IS* pcis = (PC_IS*)(pc->data);
1304: PC_BDDC* pcbddc = (PC_BDDC*)pc->data;
1305: Mat_IS *matis = (Mat_IS*)pc->pmat->data;
1306: IS is_R_local;
1307: VecType impVecType;
1308: MatType impMatType;
1309: PetscInt n_R=0;
1310: PetscInt n_D=0;
1311: PetscInt n_B=0;
1312: PetscScalar zero=0.0;
1313: PetscScalar one=1.0;
1314: PetscScalar m_one=-1.0;
1315: PetscScalar* array;
1316: PetscScalar *coarse_submat_vals;
1317: PetscInt *idx_R_local;
1318: PetscReal *coarsefunctions_errors,*constraints_errors;
1319: /* auxiliary indices */
1320: PetscInt i,j,k;
1321: /* for verbose output of bddc */
1322: PetscViewer viewer=pcbddc->dbg_viewer;
1323: PetscInt dbg_flag=pcbddc->dbg_flag;
1324: /* for counting coarse dofs */
1325: PetscInt n_vertices,n_constraints;
1326: PetscInt size_of_constraint;
1327: PetscInt *row_cmat_indices;
1328: PetscScalar *row_cmat_values;
1329: PetscInt *vertices;
1332: /* Set Non-overlapping dimensions */
1333: n_B = pcis->n_B; n_D = pcis->n - n_B;
1335: /* transform local matrices if needed */
1336: if (pcbddc->use_change_of_basis) {
1337: Mat change_mat_all;
1338: PetscInt *nnz,*is_indices,*temp_indices;
1340: PetscMalloc(pcis->n*sizeof(PetscInt),&nnz);
1341: ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);
1342: for (i=0;i<n_D;i++) nnz[is_indices[i]] = 1;
1343: ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);
1344: ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);
1345: k=1;
1346: for (i=0;i<n_B;i++) {
1347: MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);
1348: nnz[is_indices[i]]=j;
1349: if (k < j) k = j;
1350: MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);
1351: }
1352: ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);
1353: /* assemble change of basis matrix on the whole set of local dofs */
1354: PetscMalloc(k*sizeof(PetscInt),&temp_indices);
1355: MatCreate(PETSC_COMM_SELF,&change_mat_all);
1356: MatSetSizes(change_mat_all,pcis->n,pcis->n,pcis->n,pcis->n);
1357: MatSetType(change_mat_all,MATSEQAIJ);
1358: MatSeqAIJSetPreallocation(change_mat_all,0,nnz);
1359: ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);
1360: for (i=0;i<n_D;i++) {
1361: MatSetValue(change_mat_all,is_indices[i],is_indices[i],1.0,INSERT_VALUES);
1362: }
1363: ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);
1364: ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);
1365: for (i=0;i<n_B;i++) {
1366: MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);
1367: for (k=0; k<j; k++) temp_indices[k]=is_indices[row_cmat_indices[k]];
1368: MatSetValues(change_mat_all,1,&is_indices[i],j,temp_indices,row_cmat_values,INSERT_VALUES);
1369: MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);
1370: }
1371: MatAssemblyBegin(change_mat_all,MAT_FINAL_ASSEMBLY);
1372: MatAssemblyEnd(change_mat_all,MAT_FINAL_ASSEMBLY);
1373: /* TODO: HOW TO WORK WITH BAIJ? PtAP not provided */
1374: MatGetBlockSize(matis->A,&i);
1375: if (i==1) {
1376: MatPtAP(matis->A,change_mat_all,MAT_INITIAL_MATRIX,1.0,&pcbddc->local_mat);
1377: } else {
1378: Mat work_mat;
1379: MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&work_mat);
1380: MatPtAP(work_mat,change_mat_all,MAT_INITIAL_MATRIX,1.0,&pcbddc->local_mat);
1381: MatDestroy(&work_mat);
1382: }
1383: MatDestroy(&change_mat_all);
1384: PetscFree(nnz);
1385: PetscFree(temp_indices);
1386: } else {
1387: /* without change of basis, the local matrix is unchanged */
1388: PetscObjectReference((PetscObject)matis->A);
1389: pcbddc->local_mat = matis->A;
1390: }
1391: /* need to rebuild PCIS matrices during SNES or TS -> TODO move this to PCIS code */
1392: MatDestroy(&pcis->A_IB);
1393: MatDestroy(&pcis->A_BI);
1394: MatDestroy(&pcis->A_BB);
1395: MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);
1396: MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);
1397: MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);
1398: /* Change global null space passed in by the user if change of basis has been requested */
1399: if (pcbddc->NullSpace && pcbddc->use_change_of_basis) {
1400: PCBDDCNullSpaceAdaptGlobal(pc);
1401: }
1403: /* Set types for local objects needed by BDDC precondtioner */
1404: impMatType = MATSEQDENSE;
1405: impVecType = VECSEQ;
1406: /* get vertex indices from constraint matrix */
1407: PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);
1408: /* Set number of constraints */
1409: n_constraints = pcbddc->local_primal_size-n_vertices;
1410: /* Dohrmann's notation: dofs splitted in R (Remaining: all dofs but the vertices) and V (Vertices) */
1411: VecSet(pcis->vec1_N,one);
1412: VecGetArray(pcis->vec1_N,&array);
1413: for (i=0;i<n_vertices;i++) array[vertices[i]] = zero;
1414: PetscMalloc((pcis->n-n_vertices)*sizeof(PetscInt),&idx_R_local);
1415: for (i=0, n_R=0; i<pcis->n; i++) {
1416: if (array[i] == one) {
1417: idx_R_local[n_R] = i;
1418: n_R++;
1419: }
1420: }
1421: VecRestoreArray(pcis->vec1_N,&array);
1422: PetscFree(vertices);
1423: if (dbg_flag) {
1424: PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");
1425: PetscViewerFlush(viewer);
1426: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d local dimensions\n",PetscGlobalRank);
1427: PetscViewerASCIISynchronizedPrintf(viewer,"local_size = %d, dirichlet_size = %d, boundary_size = %d\n",pcis->n,n_D,n_B);
1428: PetscViewerASCIISynchronizedPrintf(viewer,"r_size = %d, v_size = %d, constraints = %d, local_primal_size = %d\n",n_R,n_vertices,n_constraints,pcbddc->local_primal_size);
1429: PetscViewerASCIISynchronizedPrintf(viewer,"pcbddc->n_vertices = %d, pcbddc->n_constraints = %d\n",pcbddc->n_vertices,pcbddc->n_constraints);
1430: PetscViewerFlush(viewer);
1431: }
1433: /* Allocate needed vectors */
1434: VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);
1435: VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);
1436: VecDuplicate(pcis->vec1_D,&pcbddc->vec4_D);
1437: VecCreate(PETSC_COMM_SELF,&pcbddc->vec1_R);
1438: VecSetSizes(pcbddc->vec1_R,n_R,n_R);
1439: VecSetType(pcbddc->vec1_R,impVecType);
1440: VecDuplicate(pcbddc->vec1_R,&pcbddc->vec2_R);
1441: VecCreate(PETSC_COMM_SELF,&pcbddc->vec1_P);
1442: VecSetSizes(pcbddc->vec1_P,pcbddc->local_primal_size,pcbddc->local_primal_size);
1443: VecSetType(pcbddc->vec1_P,impVecType);
1445: /* Creating some index sets needed */
1446: /* For submatrices */
1447: ISCreateGeneral(PETSC_COMM_SELF,n_R,idx_R_local,PETSC_OWN_POINTER,&is_R_local);
1449: /* For VecScatters pcbddc->R_to_B and (optionally) pcbddc->R_to_D */
1450: {
1451: IS is_aux1,is_aux2;
1452: PetscInt *aux_array1;
1453: PetscInt *aux_array2;
1454: PetscInt *idx_I_local;
1456: PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array1);
1457: PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array2);
1459: ISGetIndices(pcis->is_I_local,(const PetscInt**)&idx_I_local);
1460: VecGetArray(pcis->vec1_N,&array);
1461: for (i=0; i<n_D; i++) array[idx_I_local[i]] = 0;
1462: ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&idx_I_local);
1463: for (i=0, j=0; i<n_R; i++) {
1464: if (array[idx_R_local[i]] == one) {
1465: aux_array1[j] = i;
1466: j++;
1467: }
1468: }
1469: VecRestoreArray(pcis->vec1_N,&array);
1470: ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_COPY_VALUES,&is_aux1);
1471: VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1472: VecScatterEnd (pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1473: VecGetArray(pcis->vec1_B,&array);
1474: for (i=0, j=0; i<n_B; i++) {
1475: if (array[i] == one) {
1476: aux_array2[j] = i; j++;
1477: }
1478: }
1479: VecRestoreArray(pcis->vec1_B,&array);
1480: ISCreateGeneral(PETSC_COMM_SELF,j,aux_array2,PETSC_COPY_VALUES,&is_aux2);
1481: VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_B,is_aux2,&pcbddc->R_to_B);
1482: PetscFree(aux_array1);
1483: PetscFree(aux_array2);
1484: ISDestroy(&is_aux1);
1485: ISDestroy(&is_aux2);
1487: if (pcbddc->inexact_prec_type || dbg_flag ) {
1488: PetscMalloc(n_D*sizeof(PetscInt),&aux_array1);
1489: VecGetArray(pcis->vec1_N,&array);
1490: for (i=0, j=0; i<n_R; i++) {
1491: if (array[idx_R_local[i]] == zero) {
1492: aux_array1[j] = i;
1493: j++;
1494: }
1495: }
1496: VecRestoreArray(pcis->vec1_N,&array);
1497: ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_COPY_VALUES,&is_aux1);
1498: VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_D,(IS)0,&pcbddc->R_to_D);
1499: PetscFree(aux_array1);
1500: ISDestroy(&is_aux1);
1501: }
1502: }
1504: /* Creating PC contexts for local Dirichlet and Neumann problems */
1505: {
1506: Mat A_RR;
1507: PC pc_temp;
1508: MatStructure matstruct;
1509: /* Matrix for Dirichlet problem is A_II */
1510: /* HACK (TODO) A_II can be changed between nonlinear iterations */
1511: if (pc->setupcalled) { /* we dont need to rebuild dirichlet problem the first time we build BDDC */
1512: PCGetOperators(pc,NULL,NULL,&matstruct);
1513: if (matstruct == SAME_NONZERO_PATTERN) {
1514: MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,MAT_REUSE_MATRIX,&pcis->A_II);
1515: } else {
1516: MatDestroy(&pcis->A_II);
1517: MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_II);
1518: }
1519: }
1520: KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_D);
1521: PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_D,(PetscObject)pc,1);
1522: KSPSetOperators(pcbddc->ksp_D,pcis->A_II,pcis->A_II,SAME_PRECONDITIONER);
1523: KSPSetType(pcbddc->ksp_D,KSPPREONLY);
1524: KSPSetOptionsPrefix(pcbddc->ksp_D,"dirichlet_");
1525: /* default */
1526: KSPGetPC(pcbddc->ksp_D,&pc_temp);
1527: PCSetType(pc_temp,PCLU);
1528: /* Allow user's customization */
1529: KSPSetFromOptions(pcbddc->ksp_D);
1530: /* umfpack interface has a bug when matrix dimension is zero */
1531: if (!n_D) {
1532: PCSetType(pc_temp,PCNONE);
1533: }
1534: /* Set Up KSP for Dirichlet problem of BDDC */
1535: KSPSetUp(pcbddc->ksp_D);
1536: /* set ksp_D into pcis data */
1537: KSPDestroy(&pcis->ksp_D);
1538: PetscObjectReference((PetscObject)pcbddc->ksp_D);
1539: pcis->ksp_D = pcbddc->ksp_D;
1540: /* Matrix for Neumann problem is A_RR -> we need to create it */
1541: MatGetSubMatrix(pcbddc->local_mat,is_R_local,is_R_local,MAT_INITIAL_MATRIX,&A_RR);
1542: KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_R);
1543: PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_R,(PetscObject)pc,1);
1544: KSPSetOperators(pcbddc->ksp_R,A_RR,A_RR,SAME_PRECONDITIONER);
1545: KSPSetType(pcbddc->ksp_R,KSPPREONLY);
1546: KSPSetOptionsPrefix(pcbddc->ksp_R,"neumann_");
1547: /* default */
1548: KSPGetPC(pcbddc->ksp_R,&pc_temp);
1549: PCSetType(pc_temp,PCLU);
1550: /* Allow user's customization */
1551: KSPSetFromOptions(pcbddc->ksp_R);
1552: /* umfpack interface has a bug when matrix dimension is zero */
1553: if (!n_R) {
1554: PCSetType(pc_temp,PCNONE);
1555: }
1556: /* Set Up KSP for Neumann problem of BDDC */
1557: KSPSetUp(pcbddc->ksp_R);
1558: /* check Dirichlet and Neumann solvers and adapt them if a nullspace correction is needed */
1559: {
1560: Vec temp_vec;
1561: PetscReal value;
1562: PetscInt use_exact,use_exact_reduced;
1564: VecDuplicate(pcis->vec1_D,&temp_vec);
1565: VecSetRandom(pcis->vec1_D,NULL);
1566: MatMult(pcis->A_II,pcis->vec1_D,pcis->vec2_D);
1567: KSPSolve(pcbddc->ksp_D,pcis->vec2_D,temp_vec);
1568: VecAXPY(temp_vec,m_one,pcis->vec1_D);
1569: VecNorm(temp_vec,NORM_INFINITY,&value);
1570: VecDestroy(&temp_vec);
1571: use_exact = 1;
1572: if (PetscAbsReal(value) > 1.e-4) use_exact = 0;
1573: MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));
1574: PCBDDCSetUseExactDirichlet(pc,(PetscBool)use_exact_reduced);
1575: if (dbg_flag) {
1576: PetscViewerFlush(viewer);
1577: PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");
1578: PetscViewerASCIIPrintf(viewer,"Checking solution of Dirichlet and Neumann problems\n");
1579: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d infinity error for Dirichlet solve = % 1.14e \n",PetscGlobalRank,value);
1580: PetscViewerFlush(viewer);
1581: }
1582: if (n_D && pcbddc->NullSpace && !use_exact_reduced && !pcbddc->inexact_prec_type) {
1583: PCBDDCNullSpaceAssembleCorrection(pc,pcis->is_I_local);
1584: }
1585: VecDuplicate(pcbddc->vec1_R,&temp_vec);
1586: VecSetRandom(pcbddc->vec1_R,NULL);
1587: MatMult(A_RR,pcbddc->vec1_R,pcbddc->vec2_R);
1588: KSPSolve(pcbddc->ksp_R,pcbddc->vec2_R,temp_vec);
1589: VecAXPY(temp_vec,m_one,pcbddc->vec1_R);
1590: VecNorm(temp_vec,NORM_INFINITY,&value);
1591: VecDestroy(&temp_vec);
1592: use_exact = 1;
1593: if (PetscAbsReal(value) > 1.e-4) use_exact = 0;
1594: MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));
1595: if (dbg_flag) {
1596: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d infinity error for Neumann solve = % 1.14e \n",PetscGlobalRank,value);
1597: PetscViewerFlush(viewer);
1598: }
1599: if (n_R && pcbddc->NullSpace && !use_exact_reduced) { /* is it the right logic? */
1600: PCBDDCNullSpaceAssembleCorrection(pc,is_R_local);
1601: }
1602: }
1603: /* free Neumann problem's matrix */
1604: MatDestroy(&A_RR);
1605: }
1607: /* Assemble all remaining stuff needed to apply BDDC */
1608: {
1609: Mat A_RV,A_VR,A_VV;
1610: Mat M1;
1611: Mat C_CR;
1612: Mat AUXMAT;
1613: Vec vec1_C;
1614: Vec vec2_C;
1615: Vec vec1_V;
1616: Vec vec2_V;
1617: IS is_C_local,is_V_local,is_aux1;
1618: ISLocalToGlobalMapping BtoNmap;
1619: PetscInt *nnz;
1620: PetscInt *idx_V_B;
1621: PetscInt *auxindices;
1622: PetscInt index;
1623: PetscScalar* array2;
1624: MatFactorInfo matinfo;
1625: PetscBool setsym=PETSC_FALSE,issym=PETSC_FALSE;
1627: /* Allocating some extra storage just to be safe */
1628: PetscMalloc (pcis->n*sizeof(PetscInt),&nnz);
1629: PetscMalloc (pcis->n*sizeof(PetscInt),&auxindices);
1630: for (i=0;i<pcis->n;i++) auxindices[i]=i;
1632: PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);
1633: /* vertices in boundary numbering */
1634: PetscMalloc(n_vertices*sizeof(PetscInt),&idx_V_B);
1635: ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&BtoNmap);
1636: ISGlobalToLocalMappingApply(BtoNmap,IS_GTOLM_DROP,n_vertices,vertices,&i,idx_V_B);
1637: ISLocalToGlobalMappingDestroy(&BtoNmap);
1638: if (i != n_vertices) {
1639: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Error in boundary numbering for BDDC vertices! %d != %d\n",n_vertices,i);
1640: }
1642: /* some work vectors on vertices and/or constraints */
1643: if (n_vertices) {
1644: VecCreate(PETSC_COMM_SELF,&vec1_V);
1645: VecSetSizes(vec1_V,n_vertices,n_vertices);
1646: VecSetType(vec1_V,impVecType);
1647: VecDuplicate(vec1_V,&vec2_V);
1648: }
1649: if (n_constraints) {
1650: VecCreate(PETSC_COMM_SELF,&vec1_C);
1651: VecSetSizes(vec1_C,n_constraints,n_constraints);
1652: VecSetType(vec1_C,impVecType);
1653: VecDuplicate(vec1_C,&vec2_C);
1654: VecDuplicate(vec1_C,&pcbddc->vec1_C);
1655: }
1656: /* Precompute stuffs needed for preprocessing and application of BDDC*/
1657: if (n_constraints) {
1658: MatCreate(PETSC_COMM_SELF,&pcbddc->local_auxmat2);
1659: MatSetSizes(pcbddc->local_auxmat2,n_R,n_constraints,n_R,n_constraints);
1660: MatSetType(pcbddc->local_auxmat2,impMatType);
1661: MatSeqDenseSetPreallocation(pcbddc->local_auxmat2,NULL);
1663: /* Create Constraint matrix on R nodes: C_{CR} */
1664: ISCreateStride(PETSC_COMM_SELF,n_constraints,n_vertices,1,&is_C_local);
1665: MatGetSubMatrix(pcbddc->ConstraintMatrix,is_C_local,is_R_local,MAT_INITIAL_MATRIX,&C_CR);
1666: ISDestroy(&is_C_local);
1668: /* Assemble local_auxmat2 = - A_{RR}^{-1} C^T_{CR} needed by BDDC application */
1669: for (i=0;i<n_constraints;i++) {
1670: VecSet(pcbddc->vec1_R,zero);
1671: /* Get row of constraint matrix in R numbering */
1672: VecGetArray(pcbddc->vec1_R,&array);
1673: MatGetRow(C_CR,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);
1674: for (j=0;j<size_of_constraint;j++) array[row_cmat_indices[j]] = -row_cmat_values[j];
1675: MatRestoreRow(C_CR,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);
1676: VecRestoreArray(pcbddc->vec1_R,&array);
1678: /* Solve for row of constraint matrix in R numbering */
1679: KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);
1681: /* Set values */
1682: VecGetArray(pcbddc->vec2_R,&array);
1683: MatSetValues(pcbddc->local_auxmat2,n_R,auxindices,1,&i,array,INSERT_VALUES);
1684: VecRestoreArray(pcbddc->vec2_R,&array);
1685: }
1686: MatAssemblyBegin(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);
1687: MatAssemblyEnd(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);
1689: /* Assemble AUXMAT = ( LUFactor )( -C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} */
1690: MatMatMult(C_CR,pcbddc->local_auxmat2,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AUXMAT);
1691: MatFactorInfoInitialize(&matinfo);
1692: ISCreateStride(PETSC_COMM_SELF,n_constraints,0,1,&is_aux1);
1693: MatLUFactor(AUXMAT,is_aux1,is_aux1,&matinfo);
1694: ISDestroy(&is_aux1);
1696: /* Assemble explicitly M1 = ( C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} needed in preproc */
1697: MatCreate(PETSC_COMM_SELF,&M1);
1698: MatSetSizes(M1,n_constraints,n_constraints,n_constraints,n_constraints);
1699: MatSetType(M1,impMatType);
1700: MatSeqDenseSetPreallocation(M1,NULL);
1701: for (i=0;i<n_constraints;i++) {
1702: VecSet(vec1_C,zero);
1703: VecSetValue(vec1_C,i,one,INSERT_VALUES);
1704: VecAssemblyBegin(vec1_C);
1705: VecAssemblyEnd(vec1_C);
1706: MatSolve(AUXMAT,vec1_C,vec2_C);
1707: VecScale(vec2_C,m_one);
1708: VecGetArray(vec2_C,&array);
1709: MatSetValues(M1,n_constraints,auxindices,1,&i,array,INSERT_VALUES);
1710: VecRestoreArray(vec2_C,&array);
1711: }
1712: MatAssemblyBegin(M1,MAT_FINAL_ASSEMBLY);
1713: MatAssemblyEnd(M1,MAT_FINAL_ASSEMBLY);
1714: MatDestroy(&AUXMAT);
1715: /* Assemble local_auxmat1 = M1*C_{CR} needed by BDDC application in KSP and in preproc */
1716: MatMatMult(M1,C_CR,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pcbddc->local_auxmat1);
1718: }
1720: /* Get submatrices from subdomain matrix */
1721: if (n_vertices) {
1722: ISCreateGeneral(PETSC_COMM_SELF,n_vertices,vertices,PETSC_COPY_VALUES,&is_V_local);
1723: MatGetSubMatrix(pcbddc->local_mat,is_R_local,is_V_local,MAT_INITIAL_MATRIX,&A_RV);
1724: MatGetSubMatrix(pcbddc->local_mat,is_V_local,is_R_local,MAT_INITIAL_MATRIX,&A_VR);
1725: MatGetSubMatrix(pcbddc->local_mat,is_V_local,is_V_local,MAT_INITIAL_MATRIX,&A_VV);
1726: ISDestroy(&is_V_local);
1727: }
1729: /* Matrix of coarse basis functions (local) */
1730: MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_B);
1731: MatSetSizes(pcbddc->coarse_phi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);
1732: MatSetType(pcbddc->coarse_phi_B,impMatType);
1733: MatSeqDenseSetPreallocation(pcbddc->coarse_phi_B,NULL);
1734: if (pcbddc->inexact_prec_type || dbg_flag ) {
1735: MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_D);
1736: MatSetSizes(pcbddc->coarse_phi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);
1737: MatSetType(pcbddc->coarse_phi_D,impMatType);
1738: MatSeqDenseSetPreallocation(pcbddc->coarse_phi_D,NULL);
1739: }
1741: if (dbg_flag) {
1742: PetscMalloc(2*pcbddc->local_primal_size*sizeof(*coarsefunctions_errors),&coarsefunctions_errors);
1743: PetscMalloc(2*pcbddc->local_primal_size*sizeof(*constraints_errors),&constraints_errors);
1744: }
1745: /* Subdomain contribution (Non-overlapping) to coarse matrix */
1746: PetscMalloc ((pcbddc->local_primal_size)*(pcbddc->local_primal_size)*sizeof(PetscScalar),&coarse_submat_vals);
1748: /* We are now ready to evaluate coarse basis functions and subdomain contribution to coarse problem */
1749: for (i=0;i<n_vertices;i++){
1750: VecSet(vec1_V,zero);
1751: VecSetValue(vec1_V,i,one,INSERT_VALUES);
1752: VecAssemblyBegin(vec1_V);
1753: VecAssemblyEnd(vec1_V);
1754: /* solution of saddle point problem */
1755: MatMult(A_RV,vec1_V,pcbddc->vec1_R);
1756: KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec1_R);
1757: VecScale(pcbddc->vec1_R,m_one);
1758: if (n_constraints) {
1759: MatMult(pcbddc->local_auxmat1,pcbddc->vec1_R,vec1_C);
1760: MatMultAdd(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);
1761: VecScale(vec1_C,m_one);
1762: }
1763: MatMult(A_VR,pcbddc->vec1_R,vec2_V);
1764: MatMultAdd(A_VV,vec1_V,vec2_V,vec2_V);
1766: /* Set values in coarse basis function and subdomain part of coarse_mat */
1767: /* coarse basis functions */
1768: VecSet(pcis->vec1_B,zero);
1769: VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1770: VecScatterEnd (pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1771: VecGetArray(pcis->vec1_B,&array);
1772: MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&i,array,INSERT_VALUES);
1773: VecRestoreArray(pcis->vec1_B,&array);
1774: MatSetValue(pcbddc->coarse_phi_B,idx_V_B[i],i,one,INSERT_VALUES);
1775: if ( pcbddc->inexact_prec_type || dbg_flag ) {
1776: VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1777: VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1778: VecGetArray(pcis->vec1_D,&array);
1779: MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&i,array,INSERT_VALUES);
1780: VecRestoreArray(pcis->vec1_D,&array);
1781: }
1782: /* subdomain contribution to coarse matrix */
1783: VecGetArray(vec2_V,&array);
1784: for (j=0; j<n_vertices; j++) coarse_submat_vals[i*pcbddc->local_primal_size+j] = array[j]; /* WARNING -> column major ordering */
1785: VecRestoreArray(vec2_V,&array);
1786: if (n_constraints) {
1787: VecGetArray(vec1_C,&array);
1788: for (j=0; j<n_constraints; j++) coarse_submat_vals[i*pcbddc->local_primal_size+j+n_vertices] = array[j]; /* WARNING -> column major ordering */
1789: VecRestoreArray(vec1_C,&array);
1790: }
1792: if ( dbg_flag ) {
1793: /* assemble subdomain vector on nodes */
1794: VecSet(pcis->vec1_N,zero);
1795: VecGetArray(pcis->vec1_N,&array);
1796: VecGetArray(pcbddc->vec1_R,&array2);
1797: for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j];
1798: array[ vertices[i] ] = one;
1799: VecRestoreArray(pcbddc->vec1_R,&array2);
1800: VecRestoreArray(pcis->vec1_N,&array);
1801: /* assemble subdomain vector of lagrange multipliers (i.e. primal nodes) */
1802: VecSet(pcbddc->vec1_P,zero);
1803: VecGetArray(pcbddc->vec1_P,&array2);
1804: VecGetArray(vec2_V,&array);
1805: for (j=0;j<n_vertices;j++) array2[j]=array[j];
1806: VecRestoreArray(vec2_V,&array);
1807: if (n_constraints) {
1808: VecGetArray(vec1_C,&array);
1809: for (j=0;j<n_constraints;j++) array2[j+n_vertices]=array[j];
1810: VecRestoreArray(vec1_C,&array);
1811: }
1812: VecRestoreArray(pcbddc->vec1_P,&array2);
1813: VecScale(pcbddc->vec1_P,m_one);
1814: /* check saddle point solution */
1815: MatMult(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);
1816: MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);
1817: VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[i]);
1818: MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);
1819: VecGetArray(pcbddc->vec1_P,&array);
1820: array[i]=array[i]+m_one; /* shift by the identity matrix */
1821: VecRestoreArray(pcbddc->vec1_P,&array);
1822: VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[i]);
1823: }
1824: }
1826: for (i=0;i<n_constraints;i++){
1827: VecSet(vec2_C,zero);
1828: VecSetValue(vec2_C,i,m_one,INSERT_VALUES);
1829: VecAssemblyBegin(vec2_C);
1830: VecAssemblyEnd(vec2_C);
1831: /* solution of saddle point problem */
1832: MatMult(M1,vec2_C,vec1_C);
1833: MatMult(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R);
1834: VecScale(vec1_C,m_one);
1835: if (n_vertices) { MatMult(A_VR,pcbddc->vec1_R,vec2_V); }
1836: /* Set values in coarse basis function and subdomain part of coarse_mat */
1837: /* coarse basis functions */
1838: index=i+n_vertices;
1839: VecSet(pcis->vec1_B,zero);
1840: VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1841: VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1842: VecGetArray(pcis->vec1_B,&array);
1843: MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&index,array,INSERT_VALUES);
1844: VecRestoreArray(pcis->vec1_B,&array);
1845: if ( pcbddc->inexact_prec_type || dbg_flag ) {
1846: VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1847: VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1848: VecGetArray(pcis->vec1_D,&array);
1849: MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&index,array,INSERT_VALUES);
1850: VecRestoreArray(pcis->vec1_D,&array);
1851: }
1852: /* subdomain contribution to coarse matrix */
1853: if (n_vertices) {
1854: VecGetArray(vec2_V,&array);
1855: for (j=0; j<n_vertices; j++) coarse_submat_vals[index*pcbddc->local_primal_size+j]=array[j]; /* WARNING -> column major ordering */
1856: VecRestoreArray(vec2_V,&array);
1857: }
1858: VecGetArray(vec1_C,&array);
1859: for (j=0; j<n_constraints; j++) coarse_submat_vals[index*pcbddc->local_primal_size+j+n_vertices]=array[j]; /* WARNING -> column major ordering */
1860: VecRestoreArray(vec1_C,&array);
1862: if ( dbg_flag ) {
1863: /* assemble subdomain vector on nodes */
1864: VecSet(pcis->vec1_N,zero);
1865: VecGetArray(pcis->vec1_N,&array);
1866: VecGetArray(pcbddc->vec1_R,&array2);
1867: for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j];
1868: VecRestoreArray(pcbddc->vec1_R,&array2);
1869: VecRestoreArray(pcis->vec1_N,&array);
1870: /* assemble subdomain vector of lagrange multipliers */
1871: VecSet(pcbddc->vec1_P,zero);
1872: VecGetArray(pcbddc->vec1_P,&array2);
1873: if ( n_vertices) {
1874: VecGetArray(vec2_V,&array);
1875: for (j=0;j<n_vertices;j++) array2[j]=-array[j];
1876: VecRestoreArray(vec2_V,&array);
1877: }
1878: VecGetArray(vec1_C,&array);
1879: for (j=0;j<n_constraints;j++) {array2[j+n_vertices]=-array[j];}
1880: VecRestoreArray(vec1_C,&array);
1881: VecRestoreArray(pcbddc->vec1_P,&array2);
1882: /* check saddle point solution */
1883: MatMult(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);
1884: MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);
1885: VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[index]);
1886: MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);
1887: VecGetArray(pcbddc->vec1_P,&array);
1888: array[index]=array[index]+m_one; /* shift by the identity matrix */
1889: VecRestoreArray(pcbddc->vec1_P,&array);
1890: VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[index]);
1891: }
1892: }
1893: MatAssemblyBegin(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);
1894: MatAssemblyEnd(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);
1895: if ( pcbddc->inexact_prec_type || dbg_flag ) {
1896: MatAssemblyBegin(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);
1897: MatAssemblyEnd(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);
1898: }
1899: /* compute other basis functions for non-symmetric problems */
1900: MatIsSymmetricKnown(pc->pmat,&setsym,&issym);
1901: if ( !setsym || (setsym && !issym) ) {
1902: MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_psi_B);
1903: MatSetSizes(pcbddc->coarse_psi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);
1904: MatSetType(pcbddc->coarse_psi_B,impMatType);
1905: MatSeqDenseSetPreallocation(pcbddc->coarse_psi_B,NULL);
1906: if (pcbddc->inexact_prec_type || dbg_flag ) {
1907: MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_psi_D);
1908: MatSetSizes(pcbddc->coarse_psi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);
1909: MatSetType(pcbddc->coarse_psi_D,impMatType);
1910: MatSeqDenseSetPreallocation(pcbddc->coarse_psi_D,NULL);
1911: }
1912: for (i=0;i<pcbddc->local_primal_size;i++) {
1913: if (n_constraints) {
1914: VecSet(vec1_C,zero);
1915: VecGetArray(vec1_C,&array);
1916: for (j=0;j<n_constraints;j++) {
1917: array[j]=coarse_submat_vals[(j+n_vertices)*pcbddc->local_primal_size+i];
1918: }
1919: VecRestoreArray(vec1_C,&array);
1920: }
1921: if (i<n_vertices) {
1922: VecSet(vec1_V,zero);
1923: VecSetValue(vec1_V,i,m_one,INSERT_VALUES);
1924: VecAssemblyBegin(vec1_V);
1925: VecAssemblyEnd(vec1_V);
1926: MatMultTranspose(A_VR,vec1_V,pcbddc->vec1_R);
1927: if (n_constraints) {
1928: MatMultTransposeAdd(C_CR,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);
1929: }
1930: } else {
1931: MatMultTranspose(C_CR,vec1_C,pcbddc->vec1_R);
1932: }
1933: KSPSolveTranspose(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec1_R);
1934: VecSet(pcis->vec1_B,zero);
1935: VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1936: VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
1937: VecGetArray(pcis->vec1_B,&array);
1938: MatSetValues(pcbddc->coarse_psi_B,n_B,auxindices,1,&i,array,INSERT_VALUES);
1939: VecRestoreArray(pcis->vec1_B,&array);
1940: if (i<n_vertices) {
1941: MatSetValue(pcbddc->coarse_psi_B,idx_V_B[i],i,one,INSERT_VALUES);
1942: }
1943: if ( pcbddc->inexact_prec_type || dbg_flag ) {
1944: VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1945: VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
1946: VecGetArray(pcis->vec1_D,&array);
1947: MatSetValues(pcbddc->coarse_psi_D,n_D,auxindices,1,&i,array,INSERT_VALUES);
1948: VecRestoreArray(pcis->vec1_D,&array);
1949: }
1951: if ( dbg_flag ) {
1952: /* assemble subdomain vector on nodes */
1953: VecSet(pcis->vec1_N,zero);
1954: VecGetArray(pcis->vec1_N,&array);
1955: VecGetArray(pcbddc->vec1_R,&array2);
1956: for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j];
1957: if (i<n_vertices) array[vertices[i]] = one;
1958: VecRestoreArray(pcbddc->vec1_R,&array2);
1959: VecRestoreArray(pcis->vec1_N,&array);
1960: /* assemble subdomain vector of lagrange multipliers */
1961: VecGetArray(pcbddc->vec1_P,&array);
1962: for (j=0;j<pcbddc->local_primal_size;j++) {
1963: array[j]=-coarse_submat_vals[j*pcbddc->local_primal_size+i];
1964: }
1965: VecRestoreArray(pcbddc->vec1_P,&array);
1966: /* check saddle point solution */
1967: MatMultTranspose(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);
1968: MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);
1969: VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[i+pcbddc->local_primal_size]);
1970: MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);
1971: VecGetArray(pcbddc->vec1_P,&array);
1972: array[i]=array[i]+m_one; /* shift by the identity matrix */
1973: VecRestoreArray(pcbddc->vec1_P,&array);
1974: VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[i+pcbddc->local_primal_size]);
1975: }
1976: }
1977: MatAssemblyBegin(pcbddc->coarse_psi_B,MAT_FINAL_ASSEMBLY);
1978: MatAssemblyEnd(pcbddc->coarse_psi_B,MAT_FINAL_ASSEMBLY);
1979: if ( pcbddc->inexact_prec_type || dbg_flag ) {
1980: MatAssemblyBegin(pcbddc->coarse_psi_D,MAT_FINAL_ASSEMBLY);
1981: MatAssemblyEnd(pcbddc->coarse_psi_D,MAT_FINAL_ASSEMBLY);
1982: }
1983: }
1984: PetscFree(idx_V_B);
1985: /* Checking coarse_sub_mat and coarse basis functios */
1986: /* Symmetric case : It should be \Phi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */
1987: /* Non-symmetric case : It should be \Psi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */
1988: if (dbg_flag) {
1989: Mat coarse_sub_mat;
1990: Mat TM1,TM2,TM3,TM4;
1991: Mat coarse_phi_D,coarse_phi_B;
1992: Mat coarse_psi_D,coarse_psi_B;
1993: Mat A_II,A_BB,A_IB,A_BI;
1994: MatType checkmattype=MATSEQAIJ;
1995: PetscReal real_value;
1997: MatConvert(pcis->A_II,checkmattype,MAT_INITIAL_MATRIX,&A_II);
1998: MatConvert(pcis->A_IB,checkmattype,MAT_INITIAL_MATRIX,&A_IB);
1999: MatConvert(pcis->A_BI,checkmattype,MAT_INITIAL_MATRIX,&A_BI);
2000: MatConvert(pcis->A_BB,checkmattype,MAT_INITIAL_MATRIX,&A_BB);
2001: MatConvert(pcbddc->coarse_phi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_D);
2002: MatConvert(pcbddc->coarse_phi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_B);
2003: if (pcbddc->coarse_psi_B) {
2004: MatConvert(pcbddc->coarse_psi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_psi_D);
2005: MatConvert(pcbddc->coarse_psi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_psi_B);
2006: }
2007: MatCreateSeqDense(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_size,coarse_submat_vals,&coarse_sub_mat);
2009: PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");
2010: PetscViewerASCIIPrintf(viewer,"Check coarse sub mat and local basis functions\n");
2011: PetscViewerFlush(viewer);
2012: if (pcbddc->coarse_psi_B) {
2013: MatMatMult(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);
2014: MatTransposeMatMult(coarse_psi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM1);
2015: MatDestroy(&AUXMAT);
2016: MatMatMult(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);
2017: MatTransposeMatMult(coarse_psi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM2);
2018: MatDestroy(&AUXMAT);
2019: MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);
2020: MatTransposeMatMult(coarse_psi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);
2021: MatDestroy(&AUXMAT);
2022: MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);
2023: MatTransposeMatMult(coarse_psi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);
2024: MatDestroy(&AUXMAT);
2025: } else {
2026: MatPtAP(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&TM1);
2027: MatPtAP(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&TM2);
2028: MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);
2029: MatTransposeMatMult(coarse_phi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);
2030: MatDestroy(&AUXMAT);
2031: MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);
2032: MatTransposeMatMult(coarse_phi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);
2033: MatDestroy(&AUXMAT);
2034: }
2035: MatAXPY(TM1,one,TM2,DIFFERENT_NONZERO_PATTERN);
2036: MatAXPY(TM1,one,TM3,DIFFERENT_NONZERO_PATTERN);
2037: MatAXPY(TM1,one,TM4,DIFFERENT_NONZERO_PATTERN);
2038: MatConvert(TM1,MATSEQDENSE,MAT_REUSE_MATRIX,&TM1);
2039: MatAXPY(TM1,m_one,coarse_sub_mat,SAME_NONZERO_PATTERN);
2040: MatNorm(TM1,NORM_INFINITY,&real_value);
2041: PetscViewerASCIISynchronizedPrintf(viewer,"----------------------------------\n");
2042: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d \n",PetscGlobalRank);
2043: PetscViewerASCIISynchronizedPrintf(viewer,"matrix error = % 1.14e\n",real_value);
2044: PetscViewerASCIISynchronizedPrintf(viewer,"coarse functions (phi) errors\n");
2045: for (i=0;i<pcbddc->local_primal_size;i++) {
2046: PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i,coarsefunctions_errors[i]);
2047: }
2048: PetscViewerASCIISynchronizedPrintf(viewer,"constraints (phi) errors\n");
2049: for (i=0;i<pcbddc->local_primal_size;i++) {
2050: PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i,constraints_errors[i]);
2051: }
2052: if (pcbddc->coarse_psi_B) {
2053: PetscViewerASCIISynchronizedPrintf(viewer,"coarse functions (psi) errors\n");
2054: for (i=pcbddc->local_primal_size;i<2*pcbddc->local_primal_size;i++) {
2055: PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i-pcbddc->local_primal_size,coarsefunctions_errors[i]);
2056: }
2057: PetscViewerASCIISynchronizedPrintf(viewer,"constraints (psi) errors\n");
2058: for (i=pcbddc->local_primal_size;i<2*pcbddc->local_primal_size;i++) {
2059: PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i-pcbddc->local_primal_size,constraints_errors[i]);
2060: }
2061: }
2062: PetscViewerFlush(viewer);
2063: MatDestroy(&A_II);
2064: MatDestroy(&A_BB);
2065: MatDestroy(&A_IB);
2066: MatDestroy(&A_BI);
2067: MatDestroy(&TM1);
2068: MatDestroy(&TM2);
2069: MatDestroy(&TM3);
2070: MatDestroy(&TM4);
2071: MatDestroy(&coarse_phi_D);
2072: MatDestroy(&coarse_phi_B);
2073: if (pcbddc->coarse_psi_B) {
2074: MatDestroy(&coarse_psi_D);
2075: MatDestroy(&coarse_psi_B);
2076: }
2077: MatDestroy(&coarse_sub_mat);
2078: PetscFree(coarsefunctions_errors);
2079: PetscFree(constraints_errors);
2080: }
2081: /* free memory */
2082: if (n_vertices) {
2083: PetscFree(vertices);
2084: VecDestroy(&vec1_V);
2085: VecDestroy(&vec2_V);
2086: MatDestroy(&A_RV);
2087: MatDestroy(&A_VR);
2088: MatDestroy(&A_VV);
2089: }
2090: if (n_constraints) {
2091: VecDestroy(&vec1_C);
2092: VecDestroy(&vec2_C);
2093: MatDestroy(&M1);
2094: MatDestroy(&C_CR);
2095: }
2096: PetscFree(auxindices);
2097: PetscFree(nnz);
2098: /* create coarse matrix and data structures for message passing associated actual choice of coarse problem type */
2099: PCBDDCSetUpCoarseEnvironment(pc,coarse_submat_vals);
2100: PetscFree(coarse_submat_vals);
2101: }
2102: /* free memory */
2103: ISDestroy(&is_R_local);
2105: return(0);
2106: }
2108: /* -------------------------------------------------------------------------- */
2110: /* BDDC requires metis 5.0.1 for multilevel */
2111: #if defined(PETSC_HAVE_METIS)
2112: #include "metis.h"
2113: #define MetisInt idx_t
2114: #define MetisScalar real_t
2115: #endif
2119: static PetscErrorCode PCBDDCSetUpCoarseEnvironment(PC pc,PetscScalar* coarse_submat_vals)
2120: {
2123: Mat_IS *matis = (Mat_IS*)pc->pmat->data;
2124: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
2125: PC_IS *pcis = (PC_IS*)pc->data;
2126: MPI_Comm prec_comm;
2127: MPI_Comm coarse_comm;
2129: MatNullSpace CoarseNullSpace;
2131: /* common to all choiches */
2132: PetscScalar *temp_coarse_mat_vals;
2133: PetscScalar *ins_coarse_mat_vals;
2134: PetscInt *ins_local_primal_indices;
2135: PetscMPIInt *localsizes2,*localdispl2;
2136: PetscMPIInt size_prec_comm;
2137: PetscMPIInt rank_prec_comm;
2138: PetscMPIInt active_rank=MPI_PROC_NULL;
2139: PetscMPIInt master_proc=0;
2140: PetscInt ins_local_primal_size;
2141: /* specific to MULTILEVEL_BDDC */
2142: PetscMPIInt *ranks_recv=0;
2143: PetscMPIInt count_recv=0;
2144: PetscMPIInt rank_coarse_proc_send_to=-1;
2145: PetscMPIInt coarse_color = MPI_UNDEFINED;
2146: ISLocalToGlobalMapping coarse_ISLG;
2147: /* some other variables */
2149: MatType coarse_mat_type;
2150: PCType coarse_pc_type;
2151: KSPType coarse_ksp_type;
2152: PC pc_temp;
2153: PetscInt i,j,k;
2154: PetscInt max_it_coarse_ksp=1; /* don't increase this value */
2155: /* verbose output viewer */
2156: PetscViewer viewer=pcbddc->dbg_viewer;
2157: PetscInt dbg_flag=pcbddc->dbg_flag;
2159: PetscInt offset,offset2;
2160: PetscMPIInt im_active,active_procs;
2161: PetscInt *dnz,*onz;
2163: PetscBool setsym,issym=PETSC_FALSE;
2166: PetscObjectGetComm((PetscObject)pc,&prec_comm);
2167: ins_local_primal_indices = 0;
2168: ins_coarse_mat_vals = 0;
2169: localsizes2 = 0;
2170: localdispl2 = 0;
2171: temp_coarse_mat_vals = 0;
2172: coarse_ISLG = 0;
2174: MPI_Comm_size(prec_comm,&size_prec_comm);
2175: MPI_Comm_rank(prec_comm,&rank_prec_comm);
2176: MatIsSymmetricKnown(pc->pmat,&setsym,&issym);
2178: /* Assign global numbering to coarse dofs */
2179: {
2180: PetscInt *auxlocal_primal,*aux_idx;
2181: PetscMPIInt mpi_local_primal_size;
2182: PetscScalar coarsesum,*array;
2184: mpi_local_primal_size = (PetscMPIInt)pcbddc->local_primal_size;
2186: /* Construct needed data structures for message passing */
2187: j = 0;
2188: if (rank_prec_comm == 0 || pcbddc->coarse_problem_type == REPLICATED_BDDC || pcbddc->coarse_problem_type == MULTILEVEL_BDDC) {
2189: j = size_prec_comm;
2190: }
2191: PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->local_primal_sizes);
2192: PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->local_primal_displacements);
2193: /* Gather local_primal_size information for all processes */
2194: if (pcbddc->coarse_problem_type == REPLICATED_BDDC || pcbddc->coarse_problem_type == MULTILEVEL_BDDC) {
2195: MPI_Allgather(&mpi_local_primal_size,1,MPIU_INT,&pcbddc->local_primal_sizes[0],1,MPIU_INT,prec_comm);
2196: } else {
2197: MPI_Gather(&mpi_local_primal_size,1,MPIU_INT,&pcbddc->local_primal_sizes[0],1,MPIU_INT,0,prec_comm);
2198: }
2199: pcbddc->replicated_primal_size = 0;
2200: for (i=0; i<j; i++) {
2201: pcbddc->local_primal_displacements[i] = pcbddc->replicated_primal_size ;
2202: pcbddc->replicated_primal_size += pcbddc->local_primal_sizes[i];
2203: }
2205: /* First let's count coarse dofs.
2206: This code fragment assumes that the number of local constraints per connected component
2207: is not greater than the number of nodes defined for the connected component
2208: (otherwise we will surely have linear dependence between constraints and thus a singular coarse problem) */
2209: PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&auxlocal_primal);
2210: PCBDDCGetPrimalVerticesLocalIdx(pc,&i,&aux_idx);
2211: PetscMemcpy(auxlocal_primal,aux_idx,i*sizeof(PetscInt));
2212: PetscFree(aux_idx);
2213: PCBDDCGetPrimalConstraintsLocalIdx(pc,&j,&aux_idx);
2214: PetscMemcpy(&auxlocal_primal[i],aux_idx,j*sizeof(PetscInt));
2215: PetscFree(aux_idx);
2216: /* Compute number of coarse dofs */
2217: PCBDDCSubsetNumbering(prec_comm,matis->mapping,pcbddc->local_primal_size,auxlocal_primal,NULL,&pcbddc->coarse_size,&pcbddc->local_primal_indices);
2219: if (dbg_flag) {
2220: PetscViewerFlush(viewer);
2221: PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");
2222: PetscViewerASCIIPrintf(viewer,"Check coarse indices\n");
2223: VecSet(pcis->vec1_N,0.0);
2224: VecGetArray(pcis->vec1_N,&array);
2225: for (i=0;i<pcbddc->local_primal_size;i++) array[auxlocal_primal[i]]=1.0;
2226: VecRestoreArray(pcis->vec1_N,&array);
2227: VecSet(pcis->vec1_global,0.0);
2228: VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
2229: VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
2230: VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
2231: VecScatterEnd (matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
2232: VecGetArray(pcis->vec1_N,&array);
2233: for (i=0;i<pcis->n;i++) {
2234: if (array[i] == 1.0) {
2235: ISLocalToGlobalMappingApply(matis->mapping,1,&i,&j);
2236: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d: WRONG COARSE INDEX %d (local %d)\n",PetscGlobalRank,j,i);
2237: }
2238: }
2239: PetscViewerFlush(viewer);
2240: for (i=0;i<pcis->n;i++) {
2241: if (PetscRealPart(array[i]) > 0.0) array[i] = 1.0/PetscRealPart(array[i]);
2242: }
2243: VecRestoreArray(pcis->vec1_N,&array);
2244: VecSet(pcis->vec1_global,0.0);
2245: VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
2246: VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
2247: VecSum(pcis->vec1_global,&coarsesum);
2248: PetscViewerASCIIPrintf(viewer,"Size of coarse problem SHOULD be %lf\n",coarsesum);
2249: PetscViewerFlush(viewer);
2250: }
2251: PetscFree(auxlocal_primal);
2252: }
2254: if (dbg_flag) {
2255: PetscViewerASCIIPrintf(viewer,"Size of coarse problem is %d\n",pcbddc->coarse_size);
2256: if (dbg_flag > 1) {
2257: PetscViewerASCIIPrintf(viewer,"Distribution of local primal indices\n");
2258: PetscViewerFlush(viewer);
2259: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d\n",PetscGlobalRank);
2260: for (i=0;i<pcbddc->local_primal_size;i++) {
2261: PetscViewerASCIISynchronizedPrintf(viewer,"local_primal_indices[%d]=%d \n",i,pcbddc->local_primal_indices[i]);
2262: }
2263: PetscViewerFlush(viewer);
2264: }
2265: }
2267: im_active = 0;
2268: if (pcis->n) im_active = 1;
2269: MPI_Allreduce(&im_active,&active_procs,1,MPIU_INT,MPI_SUM,prec_comm);
2271: /* adapt coarse problem type */
2272: #if defined(PETSC_HAVE_METIS)
2273: if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) {
2274: if (pcbddc->current_level < pcbddc->max_levels) {
2275: if ( (active_procs/pcbddc->coarsening_ratio) < 2 ) {
2276: if (dbg_flag) {
2277: PetscViewerASCIIPrintf(viewer,"Not enough active processes on level %d (active %d,ratio %d). Parallel direct solve for coarse problem\n",pcbddc->current_level,active_procs,pcbddc->coarsening_ratio);
2278: PetscViewerFlush(viewer);
2279: }
2280: pcbddc->coarse_problem_type = PARALLEL_BDDC;
2281: }
2282: } else {
2283: if (dbg_flag) {
2284: PetscViewerASCIIPrintf(viewer,"Max number of levels reached. Using parallel direct solve for coarse problem\n",pcbddc->max_levels,active_procs,pcbddc->coarsening_ratio);
2285: PetscViewerFlush(viewer);
2286: }
2287: pcbddc->coarse_problem_type = PARALLEL_BDDC;
2288: }
2289: }
2290: #else
2291: pcbddc->coarse_problem_type = PARALLEL_BDDC;
2292: #endif
2294: switch(pcbddc->coarse_problem_type){
2296: case(MULTILEVEL_BDDC): /* we define a coarse mesh where subdomains are elements */
2297: {
2298: #if defined(PETSC_HAVE_METIS)
2299: /* we need additional variables */
2300: MetisInt n_subdomains,n_parts,objval,ncon,faces_nvtxs;
2301: MetisInt *metis_coarse_subdivision;
2302: MetisInt options[METIS_NOPTIONS];
2303: PetscMPIInt size_coarse_comm,rank_coarse_comm;
2304: PetscMPIInt procs_jumps_coarse_comm;
2305: PetscMPIInt *coarse_subdivision;
2306: PetscMPIInt *total_count_recv;
2307: PetscMPIInt *total_ranks_recv;
2308: PetscMPIInt *displacements_recv;
2309: PetscMPIInt *my_faces_connectivity;
2310: PetscMPIInt *petsc_faces_adjncy;
2311: MetisInt *faces_adjncy;
2312: MetisInt *faces_xadj;
2313: PetscMPIInt *number_of_faces;
2314: PetscMPIInt *faces_displacements;
2315: PetscInt *array_int;
2316: PetscMPIInt my_faces=0;
2317: PetscMPIInt total_faces=0;
2318: PetscInt ranks_stretching_ratio;
2320: /* define some quantities */
2321: pcbddc->coarse_communications_type = SCATTERS_BDDC;
2322: coarse_mat_type = MATIS;
2323: coarse_pc_type = PCBDDC;
2324: coarse_ksp_type = KSPRICHARDSON;
2326: /* details of coarse decomposition */
2327: n_subdomains = active_procs;
2328: n_parts = n_subdomains/pcbddc->coarsening_ratio;
2329: ranks_stretching_ratio = size_prec_comm/active_procs;
2330: procs_jumps_coarse_comm = pcbddc->coarsening_ratio*ranks_stretching_ratio;
2332: #if 0
2333: PetscMPIInt *old_ranks;
2334: PetscInt *new_ranks,*jj,*ii;
2335: MatPartitioning mat_part;
2336: IS coarse_new_decomposition,is_numbering;
2337: PetscViewer viewer_test;
2338: MPI_Comm test_coarse_comm;
2339: PetscMPIInt test_coarse_color;
2340: Mat mat_adj;
2341: /* Create new communicator for coarse problem splitting the old one */
2342: /* procs with coarse_color = MPI_UNDEFINED will have coarse_comm = MPI_COMM_NULL (from mpi standards)
2343: key = rank_prec_comm -> keep same ordering of ranks from the old to the new communicator */
2344: test_coarse_color = ( im_active ? 0 : MPI_UNDEFINED );
2345: test_coarse_comm = MPI_COMM_NULL;
2346: MPI_Comm_split(prec_comm,test_coarse_color,rank_prec_comm,&test_coarse_comm);
2347: if (im_active) {
2348: PetscMalloc(n_subdomains*sizeof(PetscMPIInt),&old_ranks);
2349: PetscMalloc(size_prec_comm*sizeof(PetscInt),&new_ranks);
2350: MPI_Comm_rank(test_coarse_comm,&rank_coarse_comm);
2351: MPI_Comm_size(test_coarse_comm,&j);
2352: MPI_Allgather(&rank_prec_comm,1,MPIU_INT,old_ranks,1,MPIU_INT,test_coarse_comm);
2353: for (i=0; i<size_prec_comm; i++) new_ranks[i] = -1;
2354: for (i=0; i<n_subdomains; i++) new_ranks[old_ranks[i]] = i;
2355: PetscViewerASCIIOpen(test_coarse_comm,"test_mat_part.out",&viewer_test);
2356: k = pcis->n_neigh-1;
2357: PetscMalloc(2*sizeof(PetscInt),&ii);
2358: ii[0]=0;
2359: ii[1]=k;
2360: PetscMalloc(k*sizeof(PetscInt),&jj);
2361: for (i=0; i<k; i++) jj[i]=new_ranks[pcis->neigh[i+1]];
2362: PetscSortInt(k,jj);
2363: MatCreateMPIAdj(test_coarse_comm,1,n_subdomains,ii,jj,NULL,&mat_adj);
2364: MatView(mat_adj,viewer_test);
2365: MatPartitioningCreate(test_coarse_comm,&mat_part);
2366: MatPartitioningSetAdjacency(mat_part,mat_adj);
2367: MatPartitioningSetFromOptions(mat_part);
2368: printf("Setting Nparts %d\n",n_parts);
2369: MatPartitioningSetNParts(mat_part,n_parts);
2370: MatPartitioningView(mat_part,viewer_test);
2371: MatPartitioningApply(mat_part,&coarse_new_decomposition);
2372: ISView(coarse_new_decomposition,viewer_test);
2373: ISPartitioningToNumbering(coarse_new_decomposition,&is_numbering);
2374: ISView(is_numbering,viewer_test);
2375: PetscViewerDestroy(&viewer_test);
2376: ISDestroy(&coarse_new_decomposition);
2377: ISDestroy(&is_numbering);
2378: MatPartitioningDestroy(&mat_part);
2379: PetscFree(old_ranks);
2380: PetscFree(new_ranks);
2381: MPI_Comm_free(&test_coarse_comm);
2382: }
2383: #endif
2385: /* build CSR graph of subdomains' connectivity */
2386: PetscMalloc (pcis->n*sizeof(PetscInt),&array_int);
2387: PetscMemzero(array_int,pcis->n*sizeof(PetscInt));
2388: for (i=1;i<pcis->n_neigh;i++){/* i=1 so I don't count myself -> faces nodes counts to 1 */
2389: for (j=0;j<pcis->n_shared[i];j++){
2390: array_int[ pcis->shared[i][j] ]+=1;
2391: }
2392: }
2393: for (i=1;i<pcis->n_neigh;i++){
2394: for (j=0;j<pcis->n_shared[i];j++){
2395: if (array_int[ pcis->shared[i][j] ] > 0 ){
2396: my_faces++;
2397: break;
2398: }
2399: }
2400: }
2402: MPI_Reduce(&my_faces,&total_faces,1,MPIU_INT,MPI_SUM,master_proc,prec_comm);
2403: PetscMalloc (my_faces*sizeof(PetscInt),&my_faces_connectivity);
2404: my_faces=0;
2405: for (i=1;i<pcis->n_neigh;i++){
2406: for (j=0;j<pcis->n_shared[i];j++){
2407: if (array_int[ pcis->shared[i][j] ] > 0 ){
2408: my_faces_connectivity[my_faces]=pcis->neigh[i];
2409: my_faces++;
2410: break;
2411: }
2412: }
2413: }
2414: if (rank_prec_comm == master_proc) {
2415: PetscMalloc (total_faces*sizeof(PetscMPIInt),&petsc_faces_adjncy);
2416: PetscMalloc (size_prec_comm*sizeof(PetscMPIInt),&number_of_faces);
2417: PetscMalloc (total_faces*sizeof(MetisInt),&faces_adjncy);
2418: PetscMalloc ((n_subdomains+1)*sizeof(MetisInt),&faces_xadj);
2419: PetscMalloc ((size_prec_comm+1)*sizeof(PetscMPIInt),&faces_displacements);
2420: }
2421: MPI_Gather(&my_faces,1,MPIU_INT,&number_of_faces[0],1,MPIU_INT,master_proc,prec_comm);
2422: if (rank_prec_comm == master_proc) {
2423: faces_xadj[0]=0;
2424: faces_displacements[0]=0;
2425: j=0;
2426: for (i=1;i<size_prec_comm+1;i++) {
2427: faces_displacements[i]=faces_displacements[i-1]+number_of_faces[i-1];
2428: if (number_of_faces[i-1]) {
2429: j++;
2430: faces_xadj[j]=faces_xadj[j-1]+number_of_faces[i-1];
2431: }
2432: }
2433: }
2434: MPI_Gatherv(&my_faces_connectivity[0],my_faces,MPIU_INT,&petsc_faces_adjncy[0],number_of_faces,faces_displacements,MPIU_INT,master_proc,prec_comm);
2435: PetscFree(my_faces_connectivity);
2436: PetscFree(array_int);
2437: if (rank_prec_comm == master_proc) {
2438: for (i=0;i<total_faces;i++) faces_adjncy[i]=(MetisInt)(petsc_faces_adjncy[i]/ranks_stretching_ratio); /* cast to MetisInt */
2439: PetscFree(faces_displacements);
2440: PetscFree(number_of_faces);
2441: PetscFree(petsc_faces_adjncy);
2442: }
2444: if ( rank_prec_comm == master_proc ) {
2446: PetscInt heuristic_for_metis=3;
2448: ncon=1;
2449: faces_nvtxs=n_subdomains;
2450: /* partition graoh induced by face connectivity */
2451: PetscMalloc (n_subdomains*sizeof(MetisInt),&metis_coarse_subdivision);
2452: METIS_SetDefaultOptions(options);
2453: /* we need a contiguous partition of the coarse mesh */
2454: options[METIS_OPTION_CONTIG]=1;
2455: options[METIS_OPTION_NITER]=30;
2456: if (pcbddc->coarsening_ratio > 1) {
2457: if (n_subdomains>n_parts*heuristic_for_metis) {
2458: options[METIS_OPTION_IPTYPE]=METIS_IPTYPE_EDGE;
2459: options[METIS_OPTION_OBJTYPE]=METIS_OBJTYPE_CUT;
2460: METIS_PartGraphKway(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision);
2461: if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in METIS_PartGraphKway (metis error code %D) called from PCBDDCSetUpCoarseEnvironment\n",ierr);
2462: } else {
2463: METIS_PartGraphRecursive(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision);
2464: if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in METIS_PartGraphRecursive (metis error code %D) called from PCBDDCSetUpCoarseEnvironment\n",ierr);
2465: }
2466: } else {
2467: for (i=0;i<n_subdomains;i++) metis_coarse_subdivision[i]=i;
2468: }
2469: PetscFree(faces_xadj);
2470: PetscFree(faces_adjncy);
2471: PetscMalloc(size_prec_comm*sizeof(PetscMPIInt),&coarse_subdivision);
2473: /* copy/cast values avoiding possible type conflicts between PETSc, MPI and METIS */
2474: for (i=0;i<size_prec_comm;i++) coarse_subdivision[i]=MPI_PROC_NULL;
2475: for (i=0;i<n_subdomains;i++) coarse_subdivision[ranks_stretching_ratio*i]=(PetscInt)(metis_coarse_subdivision[i]);
2476: PetscFree(metis_coarse_subdivision);
2477: }
2479: /* Create new communicator for coarse problem splitting the old one */
2480: if ( !(rank_prec_comm%procs_jumps_coarse_comm) && rank_prec_comm < procs_jumps_coarse_comm*n_parts ){
2481: coarse_color=0; /* for communicator splitting */
2482: active_rank=rank_prec_comm; /* for insertion of matrix values */
2483: }
2484: /* procs with coarse_color = MPI_UNDEFINED will have coarse_comm = MPI_COMM_NULL (from mpi standards)
2485: key = rank_prec_comm -> keep same ordering of ranks from the old to the new communicator */
2486: MPI_Comm_split(prec_comm,coarse_color,rank_prec_comm,&coarse_comm);
2488: if ( coarse_color == 0 ) {
2489: MPI_Comm_size(coarse_comm,&size_coarse_comm);
2490: MPI_Comm_rank(coarse_comm,&rank_coarse_comm);
2491: } else {
2492: rank_coarse_comm = MPI_PROC_NULL;
2493: }
2495: /* master proc take care of arranging and distributing coarse information */
2496: if (rank_coarse_comm == master_proc) {
2497: PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&displacements_recv);
2498: PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&total_count_recv);
2499: PetscMalloc (n_subdomains*sizeof(PetscMPIInt),&total_ranks_recv);
2500: /* some initializations */
2501: displacements_recv[0]=0;
2502: PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt));
2503: /* count from how many processes the j-th process of the coarse decomposition will receive data */
2504: for (j=0;j<size_coarse_comm;j++) {
2505: for (i=0;i<size_prec_comm;i++) {
2506: if (coarse_subdivision[i]==j) total_count_recv[j]++;
2507: }
2508: }
2509: /* displacements needed for scatterv of total_ranks_recv */
2510: for (i=1; i<size_coarse_comm; i++) displacements_recv[i]=displacements_recv[i-1]+total_count_recv[i-1];
2512: /* Now fill properly total_ranks_recv -> each coarse process will receive the ranks (in prec_comm communicator) of its friend (sending) processes */
2513: PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt));
2514: for (j=0;j<size_coarse_comm;j++) {
2515: for (i=0;i<size_prec_comm;i++) {
2516: if (coarse_subdivision[i]==j) {
2517: total_ranks_recv[displacements_recv[j]+total_count_recv[j]]=i;
2518: total_count_recv[j]+=1;
2519: }
2520: }
2521: }
2522: /*for (j=0;j<size_coarse_comm;j++) {
2523: printf("process %d in new rank will receive from %d processes (original ranks follows)\n",j,total_count_recv[j]);
2524: for (i=0;i<total_count_recv[j];i++) {
2525: printf("%d ",total_ranks_recv[displacements_recv[j]+i]);
2526: }
2527: printf("\n");
2528: }*/
2530: /* identify new decomposition in terms of ranks in the old communicator */
2531: for (i=0;i<n_subdomains;i++) {
2532: coarse_subdivision[ranks_stretching_ratio*i]=coarse_subdivision[ranks_stretching_ratio*i]*procs_jumps_coarse_comm;
2533: }
2534: /*printf("coarse_subdivision in old end new ranks\n");
2535: for (i=0;i<size_prec_comm;i++)
2536: if (coarse_subdivision[i]!=MPI_PROC_NULL) {
2537: printf("%d=(%d %d), ",i,coarse_subdivision[i],coarse_subdivision[i]/procs_jumps_coarse_comm);
2538: } else {
2539: printf("%d=(%d %d), ",i,coarse_subdivision[i],coarse_subdivision[i]);
2540: }
2541: printf("\n");*/
2542: }
2544: /* Scatter new decomposition for send details */
2545: MPI_Scatter(&coarse_subdivision[0],1,MPIU_INT,&rank_coarse_proc_send_to,1,MPIU_INT,master_proc,prec_comm);
2546: /* Scatter receiving details to members of coarse decomposition */
2547: if ( coarse_color == 0) {
2548: MPI_Scatter(&total_count_recv[0],1,MPIU_INT,&count_recv,1,MPIU_INT,master_proc,coarse_comm);
2549: PetscMalloc (count_recv*sizeof(PetscMPIInt),&ranks_recv);
2550: MPI_Scatterv(&total_ranks_recv[0],total_count_recv,displacements_recv,MPIU_INT,&ranks_recv[0],count_recv,MPIU_INT,master_proc,coarse_comm);
2551: }
2553: /*printf("I will send my matrix data to proc %d\n",rank_coarse_proc_send_to);
2554: if (coarse_color == 0) {
2555: printf("I will receive some matrix data from %d processes (ranks follows)\n",count_recv);
2556: for (i=0;i<count_recv;i++)
2557: printf("%d ",ranks_recv[i]);
2558: printf("\n");
2559: }*/
2561: if (rank_prec_comm == master_proc) {
2562: PetscFree(coarse_subdivision);
2563: PetscFree(total_count_recv);
2564: PetscFree(total_ranks_recv);
2565: PetscFree(displacements_recv);
2566: }
2567: #endif
2568: break;
2569: }
2571: case(REPLICATED_BDDC):
2573: pcbddc->coarse_communications_type = GATHERS_BDDC;
2574: coarse_mat_type = MATSEQAIJ;
2575: coarse_pc_type = PCLU;
2576: coarse_ksp_type = KSPPREONLY;
2577: coarse_comm = PETSC_COMM_SELF;
2578: active_rank = rank_prec_comm;
2579: break;
2581: case(PARALLEL_BDDC):
2583: pcbddc->coarse_communications_type = SCATTERS_BDDC;
2584: coarse_mat_type = MATAIJ;
2585: coarse_pc_type = PCREDUNDANT;
2586: coarse_ksp_type = KSPPREONLY;
2587: coarse_comm = prec_comm;
2588: active_rank = rank_prec_comm;
2589: break;
2591: case(SEQUENTIAL_BDDC):
2592: pcbddc->coarse_communications_type = GATHERS_BDDC;
2593: coarse_mat_type = MATAIJ;
2594: coarse_pc_type = PCLU;
2595: coarse_ksp_type = KSPPREONLY;
2596: coarse_comm = PETSC_COMM_SELF;
2597: active_rank = master_proc;
2598: break;
2599: }
2601: switch(pcbddc->coarse_communications_type){
2603: case(SCATTERS_BDDC):
2604: {
2605: if (pcbddc->coarse_problem_type==MULTILEVEL_BDDC) {
2607: IS coarse_IS;
2609: if(pcbddc->coarsening_ratio == 1) {
2610: ins_local_primal_size = pcbddc->local_primal_size;
2611: ins_local_primal_indices = pcbddc->local_primal_indices;
2612: if (coarse_color == 0) { PetscFree(ranks_recv); }
2613: /* nonzeros */
2614: PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&dnz);
2615: PetscMemzero(dnz,ins_local_primal_size*sizeof(PetscInt));
2616: for (i=0;i<ins_local_primal_size;i++) {
2617: dnz[i] = ins_local_primal_size;
2618: }
2619: } else {
2620: PetscMPIInt send_size;
2621: PetscMPIInt *send_buffer;
2622: PetscInt *aux_ins_indices;
2623: PetscInt ii,jj;
2624: MPI_Request *requests;
2626: PetscMalloc(count_recv*sizeof(PetscMPIInt),&localdispl2);
2627: /* reusing pcbddc->local_primal_displacements and pcbddc->replicated_primal_size */
2628: PetscFree(pcbddc->local_primal_displacements);
2629: PetscMalloc((count_recv+1)*sizeof(PetscMPIInt),&pcbddc->local_primal_displacements);
2630: pcbddc->replicated_primal_size = count_recv;
2631: j = 0;
2632: for (i=0;i<count_recv;i++) {
2633: pcbddc->local_primal_displacements[i] = j;
2634: j += pcbddc->local_primal_sizes[ranks_recv[i]];
2635: }
2636: pcbddc->local_primal_displacements[count_recv] = j;
2637: PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);
2638: /* allocate auxiliary space */
2639: PetscMalloc(count_recv*sizeof(PetscMPIInt),&localsizes2);
2640: PetscMalloc(pcbddc->coarse_size*sizeof(PetscInt),&aux_ins_indices);
2641: PetscMemzero(aux_ins_indices,pcbddc->coarse_size*sizeof(PetscInt));
2642: /* allocate stuffs for message massing */
2643: PetscMalloc((count_recv+1)*sizeof(MPI_Request),&requests);
2644: for (i=0;i<count_recv+1;i++) { requests[i]=MPI_REQUEST_NULL; }
2645: /* send indices to be inserted */
2646: for (i=0;i<count_recv;i++) {
2647: send_size = pcbddc->local_primal_sizes[ranks_recv[i]];
2648: MPI_Irecv(&pcbddc->replicated_local_primal_indices[pcbddc->local_primal_displacements[i]],send_size,MPIU_INT,ranks_recv[i],999,prec_comm,&requests[i]);
2649: }
2650: if (rank_coarse_proc_send_to != MPI_PROC_NULL ) {
2651: send_size = pcbddc->local_primal_size;
2652: PetscMalloc(send_size*sizeof(PetscMPIInt),&send_buffer);
2653: for (i=0;i<send_size;i++) {
2654: send_buffer[i]=(PetscMPIInt)pcbddc->local_primal_indices[i];
2655: }
2656: MPI_Isend(send_buffer,send_size,MPIU_INT,rank_coarse_proc_send_to,999,prec_comm,&requests[count_recv]);
2657: }
2658: MPI_Waitall(count_recv+1,requests,MPI_STATUSES_IGNORE);
2659: if (rank_coarse_proc_send_to != MPI_PROC_NULL ) {
2660: PetscFree(send_buffer);
2661: }
2662: j = 0;
2663: for (i=0;i<count_recv;i++) {
2664: ii = pcbddc->local_primal_displacements[i+1]-pcbddc->local_primal_displacements[i];
2665: localsizes2[i] = ii*ii;
2666: localdispl2[i] = j;
2667: j += localsizes2[i];
2668: jj = pcbddc->local_primal_displacements[i];
2669: /* it counts the coarse subdomains sharing the coarse node */
2670: for (k=0;k<ii;k++) {
2671: aux_ins_indices[pcbddc->replicated_local_primal_indices[jj+k]] += 1;
2672: }
2673: }
2674: /* temp_coarse_mat_vals used to store matrix values to be received */
2675: PetscMalloc(j*sizeof(PetscScalar),&temp_coarse_mat_vals);
2676: /* evaluate how many values I will insert in coarse mat */
2677: ins_local_primal_size = 0;
2678: for (i=0;i<pcbddc->coarse_size;i++) {
2679: if (aux_ins_indices[i]) {
2680: ins_local_primal_size++;
2681: }
2682: }
2683: /* evaluate indices I will insert in coarse mat */
2684: PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);
2685: j = 0;
2686: for(i=0;i<pcbddc->coarse_size;i++) {
2687: if(aux_ins_indices[i]) {
2688: ins_local_primal_indices[j] = i;
2689: j++;
2690: }
2691: }
2692: /* processes partecipating in coarse problem receive matrix data from their friends */
2693: for (i=0;i<count_recv;i++) {
2694: MPI_Irecv(&temp_coarse_mat_vals[localdispl2[i]],localsizes2[i],MPIU_SCALAR,ranks_recv[i],666,prec_comm,&requests[i]);
2695: }
2696: if (rank_coarse_proc_send_to != MPI_PROC_NULL ) {
2697: send_size = pcbddc->local_primal_size*pcbddc->local_primal_size;
2698: MPI_Isend(&coarse_submat_vals[0],send_size,MPIU_SCALAR,rank_coarse_proc_send_to,666,prec_comm,&requests[count_recv]);
2699: }
2700: MPI_Waitall(count_recv+1,requests,MPI_STATUSES_IGNORE);
2701: /* nonzeros */
2702: PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&dnz);
2703: PetscMemzero(dnz,ins_local_primal_size*sizeof(PetscInt));
2704: /* use aux_ins_indices to realize a global to local mapping */
2705: j=0;
2706: for(i=0;i<pcbddc->coarse_size;i++){
2707: if(aux_ins_indices[i]==0){
2708: aux_ins_indices[i]=-1;
2709: } else {
2710: aux_ins_indices[i]=j;
2711: j++;
2712: }
2713: }
2714: for (i=0;i<count_recv;i++) {
2715: j = pcbddc->local_primal_sizes[ranks_recv[i]];
2716: for (k=0;k<j;k++) {
2717: dnz[aux_ins_indices[pcbddc->replicated_local_primal_indices[pcbddc->local_primal_displacements[i]+k]]] += j;
2718: }
2719: }
2720: /* check */
2721: for (i=0;i<ins_local_primal_size;i++) {
2722: if (dnz[i] > ins_local_primal_size) {
2723: dnz[i] = ins_local_primal_size;
2724: }
2725: }
2726: PetscFree(requests);
2727: PetscFree(aux_ins_indices);
2728: if (coarse_color == 0) { PetscFree(ranks_recv); }
2729: }
2730: /* create local to global mapping needed by coarse MATIS */
2731: if (coarse_comm != MPI_COMM_NULL ) {MPI_Comm_free(&coarse_comm);}
2732: coarse_comm = prec_comm;
2733: active_rank = rank_prec_comm;
2734: ISCreateGeneral(coarse_comm,ins_local_primal_size,ins_local_primal_indices,PETSC_COPY_VALUES,&coarse_IS);
2735: ISLocalToGlobalMappingCreateIS(coarse_IS,&coarse_ISLG);
2736: ISDestroy(&coarse_IS);
2737: } else if (pcbddc->coarse_problem_type==PARALLEL_BDDC) {
2738: /* arrays for values insertion */
2739: ins_local_primal_size = pcbddc->local_primal_size;
2740: PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);
2741: PetscMalloc(ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);
2742: for (j=0;j<ins_local_primal_size;j++){
2743: ins_local_primal_indices[j]=pcbddc->local_primal_indices[j];
2744: for (i=0;i<ins_local_primal_size;i++) {
2745: ins_coarse_mat_vals[j*ins_local_primal_size+i]=coarse_submat_vals[j*ins_local_primal_size+i];
2746: }
2747: }
2748: }
2749: break;
2751: }
2753: case(GATHERS_BDDC):
2754: {
2756: PetscMPIInt mysize,mysize2;
2757: PetscMPIInt *send_buffer;
2759: if (rank_prec_comm==active_rank) {
2760: PetscMalloc ( pcbddc->replicated_primal_size*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);
2761: PetscMalloc ( pcbddc->replicated_primal_size*sizeof(PetscScalar),&pcbddc->replicated_local_primal_values);
2762: PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localsizes2);
2763: PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localdispl2);
2764: /* arrays for values insertion */
2765: for (i=0;i<size_prec_comm;i++) localsizes2[i]=pcbddc->local_primal_sizes[i]*pcbddc->local_primal_sizes[i];
2766: localdispl2[0]=0;
2767: for (i=1;i<size_prec_comm;i++) localdispl2[i]=localsizes2[i-1]+localdispl2[i-1];
2768: j=0;
2769: for (i=0;i<size_prec_comm;i++) j+=localsizes2[i];
2770: PetscMalloc ( j*sizeof(PetscScalar),&temp_coarse_mat_vals);
2771: }
2773: mysize=pcbddc->local_primal_size;
2774: mysize2=pcbddc->local_primal_size*pcbddc->local_primal_size;
2775: PetscMalloc(mysize*sizeof(PetscMPIInt),&send_buffer);
2776: for (i=0; i<mysize; i++) send_buffer[i]=(PetscMPIInt)pcbddc->local_primal_indices[i];
2778: if (pcbddc->coarse_problem_type == SEQUENTIAL_BDDC){
2779: MPI_Gatherv(send_buffer,mysize,MPIU_INT,&pcbddc->replicated_local_primal_indices[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,master_proc,prec_comm);
2780: MPI_Gatherv(&coarse_submat_vals[0],mysize2,MPIU_SCALAR,&temp_coarse_mat_vals[0],localsizes2,localdispl2,MPIU_SCALAR,master_proc,prec_comm);
2781: } else {
2782: MPI_Allgatherv(send_buffer,mysize,MPIU_INT,&pcbddc->replicated_local_primal_indices[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,prec_comm);
2783: MPI_Allgatherv(&coarse_submat_vals[0],mysize2,MPIU_SCALAR,&temp_coarse_mat_vals[0],localsizes2,localdispl2,MPIU_SCALAR,prec_comm);
2784: }
2785: PetscFree(send_buffer);
2786: break;
2787: }/* switch on coarse problem and communications associated with finished */
2788: }
2790: /* Now create and fill up coarse matrix */
2791: if ( rank_prec_comm == active_rank ) {
2793: Mat matis_coarse_local_mat;
2795: if (pcbddc->coarse_problem_type != MULTILEVEL_BDDC) {
2796: MatCreate(coarse_comm,&pcbddc->coarse_mat);
2797: MatSetSizes(pcbddc->coarse_mat,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size);
2798: MatSetType(pcbddc->coarse_mat,coarse_mat_type);
2799: MatSetOptionsPrefix(pcbddc->coarse_mat,"coarse_");
2800: MatSetFromOptions(pcbddc->coarse_mat);
2801: MatSetUp(pcbddc->coarse_mat);
2802: MatSetOption(pcbddc->coarse_mat,MAT_ROW_ORIENTED,PETSC_FALSE); /* local values stored in column major */
2803: MatSetOption(pcbddc->coarse_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
2804: } else {
2805: MatCreateIS(coarse_comm,1,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size,coarse_ISLG,&pcbddc->coarse_mat);
2806: MatSetUp(pcbddc->coarse_mat);
2807: MatISGetLocalMat(pcbddc->coarse_mat,&matis_coarse_local_mat);
2808: MatSetOptionsPrefix(pcbddc->coarse_mat,"coarse_");
2809: MatSetFromOptions(pcbddc->coarse_mat);
2810: MatSetUp(matis_coarse_local_mat);
2811: MatSetOption(matis_coarse_local_mat,MAT_ROW_ORIENTED,PETSC_FALSE); /* local values stored in column major */
2812: MatSetOption(matis_coarse_local_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
2813: }
2814: /* preallocation */
2815: if (pcbddc->coarse_problem_type != MULTILEVEL_BDDC) {
2817: PetscInt lrows,lcols,bs;
2819: MatGetLocalSize(pcbddc->coarse_mat,&lrows,&lcols);
2820: MatPreallocateInitialize(coarse_comm,lrows,lcols,dnz,onz);
2821: MatGetBlockSize(pcbddc->coarse_mat,&bs);
2823: if (pcbddc->coarse_problem_type == PARALLEL_BDDC) {
2825: Vec vec_dnz,vec_onz;
2826: PetscScalar *my_dnz,*my_onz,*array;
2827: PetscInt *mat_ranges,*row_ownership;
2828: PetscInt coarse_index_row,coarse_index_col,owner;
2830: VecCreate(prec_comm,&vec_dnz);
2831: VecSetBlockSize(vec_dnz,bs);
2832: VecSetSizes(vec_dnz,PETSC_DECIDE,pcbddc->coarse_size);
2833: VecSetType(vec_dnz,VECMPI);
2834: VecDuplicate(vec_dnz,&vec_onz);
2836: PetscMalloc(pcbddc->local_primal_size*sizeof(PetscScalar),&my_dnz);
2837: PetscMalloc(pcbddc->local_primal_size*sizeof(PetscScalar),&my_onz);
2838: PetscMemzero(my_dnz,pcbddc->local_primal_size*sizeof(PetscScalar));
2839: PetscMemzero(my_onz,pcbddc->local_primal_size*sizeof(PetscScalar));
2841: PetscMalloc(pcbddc->coarse_size*sizeof(PetscInt),&row_ownership);
2842: MatGetOwnershipRanges(pcbddc->coarse_mat,(const PetscInt**)&mat_ranges);
2843: for (i=0;i<size_prec_comm;i++) {
2844: for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
2845: row_ownership[j]=i;
2846: }
2847: }
2849: for (i=0;i<pcbddc->local_primal_size;i++) {
2850: coarse_index_row = pcbddc->local_primal_indices[i];
2851: owner = row_ownership[coarse_index_row];
2852: for (j=i;j<pcbddc->local_primal_size;j++) {
2853: owner = row_ownership[coarse_index_row];
2854: coarse_index_col = pcbddc->local_primal_indices[j];
2855: if (coarse_index_col > mat_ranges[owner]-1 && coarse_index_col < mat_ranges[owner+1] ) {
2856: my_dnz[i] += 1.0;
2857: } else {
2858: my_onz[i] += 1.0;
2859: }
2860: if (i != j) {
2861: owner = row_ownership[coarse_index_col];
2862: if (coarse_index_row > mat_ranges[owner]-1 && coarse_index_row < mat_ranges[owner+1] ) {
2863: my_dnz[j] += 1.0;
2864: } else {
2865: my_onz[j] += 1.0;
2866: }
2867: }
2868: }
2869: }
2870: VecSet(vec_dnz,0.0);
2871: VecSet(vec_onz,0.0);
2872: if (pcbddc->local_primal_size) {
2873: VecSetValues(vec_dnz,pcbddc->local_primal_size,pcbddc->local_primal_indices,my_dnz,ADD_VALUES);
2874: VecSetValues(vec_onz,pcbddc->local_primal_size,pcbddc->local_primal_indices,my_onz,ADD_VALUES);
2875: }
2876: VecAssemblyBegin(vec_dnz);
2877: VecAssemblyBegin(vec_onz);
2878: VecAssemblyEnd(vec_dnz);
2879: VecAssemblyEnd(vec_onz);
2880: j = mat_ranges[rank_prec_comm+1]-mat_ranges[rank_prec_comm];
2881: VecGetArray(vec_dnz,&array);
2882: for (i=0; i<j; i++) dnz[i] = (PetscInt)PetscRealPart(array[i]);
2884: VecRestoreArray(vec_dnz,&array);
2885: VecGetArray(vec_onz,&array);
2886: for (i=0;i<j;i++) onz[i] = (PetscInt)PetscRealPart(array[i]);
2888: VecRestoreArray(vec_onz,&array);
2889: PetscFree(my_dnz);
2890: PetscFree(my_onz);
2891: PetscFree(row_ownership);
2892: VecDestroy(&vec_dnz);
2893: VecDestroy(&vec_onz);
2894: } else {
2895: for (k=0;k<size_prec_comm;k++){
2896: offset=pcbddc->local_primal_displacements[k];
2897: offset2=localdispl2[k];
2898: ins_local_primal_size = pcbddc->local_primal_sizes[k];
2899: PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);
2900: for (j=0;j<ins_local_primal_size;j++){
2901: ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j];
2902: }
2903: for (j=0;j<ins_local_primal_size;j++) {
2904: MatPreallocateSet(ins_local_primal_indices[j],ins_local_primal_size,ins_local_primal_indices,dnz,onz);
2905: }
2906: PetscFree(ins_local_primal_indices);
2907: }
2908: }
2910: /* check */
2911: for (i=0;i<lrows;i++) {
2912: if (dnz[i]>lcols) dnz[i]=lcols;
2913: if (onz[i]>pcbddc->coarse_size-lcols) onz[i]=pcbddc->coarse_size-lcols;
2914: }
2915: MatSeqAIJSetPreallocation(pcbddc->coarse_mat,0,dnz);
2916: MatMPIAIJSetPreallocation(pcbddc->coarse_mat,0,dnz,0,onz);
2917: MatPreallocateFinalize(dnz,onz);
2918: } else {
2919: MatSeqAIJSetPreallocation(matis_coarse_local_mat,0,dnz);
2920: PetscFree(dnz);
2921: }
2922: /* insert values */
2923: if (pcbddc->coarse_problem_type == PARALLEL_BDDC) {
2924: MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);
2925: } else if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) {
2926: if (pcbddc->coarsening_ratio == 1) {
2927: ins_coarse_mat_vals = coarse_submat_vals;
2928: MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,INSERT_VALUES);
2929: } else {
2930: PetscFree(ins_local_primal_indices);
2931: for (k=0;k<pcbddc->replicated_primal_size;k++) {
2932: offset = pcbddc->local_primal_displacements[k];
2933: offset2 = localdispl2[k];
2934: ins_local_primal_size = pcbddc->local_primal_displacements[k+1]-pcbddc->local_primal_displacements[k];
2935: PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);
2936: for (j=0;j<ins_local_primal_size;j++){
2937: ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j];
2938: }
2939: ins_coarse_mat_vals = &temp_coarse_mat_vals[offset2];
2940: MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);
2941: PetscFree(ins_local_primal_indices);
2942: }
2943: }
2944: ins_local_primal_indices = 0;
2945: ins_coarse_mat_vals = 0;
2946: } else {
2947: for (k=0;k<size_prec_comm;k++){
2948: offset=pcbddc->local_primal_displacements[k];
2949: offset2=localdispl2[k];
2950: ins_local_primal_size = pcbddc->local_primal_sizes[k];
2951: PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);
2952: for (j=0;j<ins_local_primal_size;j++){
2953: ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j];
2954: }
2955: ins_coarse_mat_vals = &temp_coarse_mat_vals[offset2];
2956: MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);
2957: PetscFree(ins_local_primal_indices);
2958: }
2959: ins_local_primal_indices = 0;
2960: ins_coarse_mat_vals = 0;
2961: }
2962: MatAssemblyBegin(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);
2963: MatAssemblyEnd(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);
2964: /* symmetry of coarse matrix */
2965: if (issym) {
2966: MatSetOption(pcbddc->coarse_mat,MAT_SYMMETRIC,PETSC_TRUE);
2967: }
2968: MatGetVecs(pcbddc->coarse_mat,&pcbddc->coarse_vec,&pcbddc->coarse_rhs);
2969: }
2971: /* create loc to glob scatters if needed */
2972: if (pcbddc->coarse_communications_type == SCATTERS_BDDC) {
2973: IS local_IS,global_IS;
2974: ISCreateStride(PETSC_COMM_SELF,pcbddc->local_primal_size,0,1,&local_IS);
2975: ISCreateGeneral(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_indices,PETSC_COPY_VALUES,&global_IS);
2976: VecScatterCreate(pcbddc->vec1_P,local_IS,pcbddc->coarse_vec,global_IS,&pcbddc->coarse_loc_to_glob);
2977: ISDestroy(&local_IS);
2978: ISDestroy(&global_IS);
2979: }
2981: /* free memory no longer needed */
2982: if (coarse_ISLG) { ISLocalToGlobalMappingDestroy(&coarse_ISLG); }
2983: if (ins_local_primal_indices) { PetscFree(ins_local_primal_indices); }
2984: if (ins_coarse_mat_vals) { PetscFree(ins_coarse_mat_vals); }
2985: if (localsizes2) { PetscFree(localsizes2); }
2986: if (localdispl2) { PetscFree(localdispl2); }
2987: if (temp_coarse_mat_vals) { PetscFree(temp_coarse_mat_vals); }
2989: /* Compute coarse null space */
2990: CoarseNullSpace = 0;
2991: if (pcbddc->NullSpace) {
2992: PCBDDCNullSpaceAssembleCoarse(pc,&CoarseNullSpace);
2993: }
2995: /* KSP for coarse problem */
2996: if (rank_prec_comm == active_rank) {
2997: PetscBool isbddc=PETSC_FALSE;
2999: KSPCreate(coarse_comm,&pcbddc->coarse_ksp);
3000: PetscObjectIncrementTabLevel((PetscObject)pcbddc->coarse_ksp,(PetscObject)pc,1);
3001: KSPSetOperators(pcbddc->coarse_ksp,pcbddc->coarse_mat,pcbddc->coarse_mat,SAME_PRECONDITIONER);
3002: KSPSetTolerances(pcbddc->coarse_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,max_it_coarse_ksp);
3003: KSPSetType(pcbddc->coarse_ksp,coarse_ksp_type);
3004: KSPGetPC(pcbddc->coarse_ksp,&pc_temp);
3005: PCSetType(pc_temp,coarse_pc_type);
3006: /* Allow user's customization */
3007: KSPSetOptionsPrefix(pcbddc->coarse_ksp,"coarse_");
3008: /* Set Up PC for coarse problem BDDC */
3009: if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) {
3010: i = pcbddc->current_level+1;
3011: PCBDDCSetLevel(pc_temp,i);
3012: PCBDDCSetCoarseningRatio(pc_temp,pcbddc->coarsening_ratio);
3013: PCBDDCSetMaxLevels(pc_temp,pcbddc->max_levels);
3014: PCBDDCSetCoarseProblemType(pc_temp,MULTILEVEL_BDDC);
3015: if (CoarseNullSpace) {
3016: PCBDDCSetNullSpace(pc_temp,CoarseNullSpace);
3017: }
3018: if (dbg_flag) {
3019: PetscViewerASCIIPrintf(viewer,"----------------Level %d: Setting up level %d---------------\n",pcbddc->current_level,i);
3020: PetscViewerFlush(viewer);
3021: }
3022: } else {
3023: if (CoarseNullSpace) {
3024: KSPSetNullSpace(pcbddc->coarse_ksp,CoarseNullSpace);
3025: }
3026: }
3027: KSPSetFromOptions(pcbddc->coarse_ksp);
3028: KSPSetUp(pcbddc->coarse_ksp);
3030: KSPGetTolerances(pcbddc->coarse_ksp,NULL,NULL,NULL,&j);
3031: KSPGetPC(pcbddc->coarse_ksp,&pc_temp);
3032: PetscObjectTypeCompare((PetscObject)pc_temp,PCBDDC,&isbddc);
3033: if (j == 1) {
3034: KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_NONE);
3035: if (isbddc) {
3036: PCBDDCSetUseExactDirichlet(pc_temp,PETSC_FALSE);
3037: }
3038: }
3039: }
3040: /* Check coarse problem if requested */
3041: if ( dbg_flag && rank_prec_comm == active_rank ) {
3042: KSP check_ksp;
3043: PC check_pc;
3044: Vec check_vec;
3045: PetscReal abs_infty_error,infty_error,lambda_min,lambda_max;
3046: KSPType check_ksp_type;
3048: /* Create ksp object suitable for extreme eigenvalues' estimation */
3049: KSPCreate(coarse_comm,&check_ksp);
3050: KSPSetOperators(check_ksp,pcbddc->coarse_mat,pcbddc->coarse_mat,SAME_PRECONDITIONER);
3051: KSPSetTolerances(check_ksp,1.e-12,1.e-12,PETSC_DEFAULT,pcbddc->coarse_size);
3052: if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) {
3053: if (issym) check_ksp_type = KSPCG;
3054: else check_ksp_type = KSPGMRES;
3055: KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);
3056: } else {
3057: check_ksp_type = KSPPREONLY;
3058: }
3059: KSPSetType(check_ksp,check_ksp_type);
3060: KSPGetPC(pcbddc->coarse_ksp,&check_pc);
3061: KSPSetPC(check_ksp,check_pc);
3062: KSPSetUp(check_ksp);
3063: /* create random vec */
3064: VecDuplicate(pcbddc->coarse_vec,&check_vec);
3065: VecSetRandom(check_vec,NULL);
3066: if (CoarseNullSpace) {
3067: MatNullSpaceRemove(CoarseNullSpace,check_vec);
3068: }
3069: MatMult(pcbddc->coarse_mat,check_vec,pcbddc->coarse_rhs);
3070: /* solve coarse problem */
3071: KSPSolve(check_ksp,pcbddc->coarse_rhs,pcbddc->coarse_vec);
3072: if (CoarseNullSpace) {
3073: MatNullSpaceRemove(CoarseNullSpace,pcbddc->coarse_vec);
3074: }
3075: /* check coarse problem residual error */
3076: VecAXPY(check_vec,-1.0,pcbddc->coarse_vec);
3077: VecNorm(check_vec,NORM_INFINITY,&infty_error);
3078: MatMult(pcbddc->coarse_mat,check_vec,pcbddc->coarse_rhs);
3079: VecNorm(pcbddc->coarse_rhs,NORM_INFINITY,&abs_infty_error);
3080: VecDestroy(&check_vec);
3081: /* get eigenvalue estimation if inexact */
3082: if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) {
3083: KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);
3084: KSPGetIterationNumber(check_ksp,&k);
3085: PetscViewerASCIIPrintf(viewer,"Coarse problem eigenvalues estimated with %d iterations of %s.\n",k,check_ksp_type);
3086: PetscViewerASCIIPrintf(viewer,"Coarse problem eigenvalues: % 1.14e %1.14e\n",lambda_min,lambda_max);
3087: }
3088: PetscViewerASCIIPrintf(viewer,"Coarse problem exact infty_error : %1.14e\n",infty_error);
3089: PetscViewerASCIIPrintf(viewer,"Coarse problem residual infty_error: %1.14e\n",abs_infty_error);
3090: KSPDestroy(&check_ksp);
3091: }
3092: if (dbg_flag) {
3093: PetscViewerFlush(viewer);
3094: }
3095: MatNullSpaceDestroy(&CoarseNullSpace);
3097: return(0);
3098: }