xref: /petsc/src/ksp/pc/impls/bddc/bddcfetidp.c (revision 94ef8dde638caef1d0cd84a7dc8a2db65fcda8b6) 
1 #include <../src/ksp/pc/impls/bddc/bddc.h>
2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3 #include <petscblaslapack.h>
4 
5 static PetscErrorCode MatMult_BDdelta_deluxe_nonred(Mat A, Vec x, Vec y)
6 {
7   BDdelta_DN     ctx;
8   PetscErrorCode ierr;
9 
10   PetscFunctionBegin;
11   ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
12   ierr = MatMultTranspose(ctx->BD,x,ctx->work);CHKERRQ(ierr);
13   ierr = KSPSolveTranspose(ctx->kBD,ctx->work,y);CHKERRQ(ierr);
14   PetscFunctionReturn(0);
15 }
16 
17 static PetscErrorCode MatMultTranspose_BDdelta_deluxe_nonred(Mat A, Vec x, Vec y)
18 {
19   BDdelta_DN     ctx;
20   PetscErrorCode ierr;
21 
22   PetscFunctionBegin;
23   ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
24   ierr = KSPSolve(ctx->kBD,x,ctx->work);CHKERRQ(ierr);
25   ierr = MatMult(ctx->BD,ctx->work,y);CHKERRQ(ierr);
26   PetscFunctionReturn(0);
27 }
28 
29 static PetscErrorCode MatDestroy_BDdelta_deluxe_nonred(Mat A)
30 {
31   BDdelta_DN     ctx;
32   PetscErrorCode ierr;
33 
34   PetscFunctionBegin;
35   ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
36   ierr = MatDestroy(&ctx->BD);CHKERRQ(ierr);
37   ierr = KSPDestroy(&ctx->kBD);CHKERRQ(ierr);
38   ierr = VecDestroy(&ctx->work);CHKERRQ(ierr);
39   ierr = PetscFree(ctx);CHKERRQ(ierr);
40   PetscFunctionReturn(0);
41 }
42 
43 
44 PetscErrorCode PCBDDCCreateFETIDPMatContext(PC pc, FETIDPMat_ctx *fetidpmat_ctx)
45 {
46   FETIDPMat_ctx  newctx;
47   PetscErrorCode ierr;
48 
49   PetscFunctionBegin;
50   ierr = PetscNew(&newctx);CHKERRQ(ierr);
51   /* increase the reference count for BDDC preconditioner */
52   ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
53   newctx->pc              = pc;
54   *fetidpmat_ctx          = newctx;
55   PetscFunctionReturn(0);
56 }
57 
58 PetscErrorCode PCBDDCCreateFETIDPPCContext(PC pc, FETIDPPC_ctx *fetidppc_ctx)
59 {
60   FETIDPPC_ctx   newctx;
61   PetscErrorCode ierr;
62 
63   PetscFunctionBegin;
64   ierr = PetscNew(&newctx);CHKERRQ(ierr);
65   /* increase the reference count for BDDC preconditioner */
66   ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
67   newctx->pc              = pc;
68   *fetidppc_ctx           = newctx;
69   PetscFunctionReturn(0);
70 }
71 
72 PetscErrorCode PCBDDCDestroyFETIDPMat(Mat A)
73 {
74   FETIDPMat_ctx  mat_ctx;
75   PetscErrorCode ierr;
76 
77   PetscFunctionBegin;
78   ierr = MatShellGetContext(A,(void**)&mat_ctx);CHKERRQ(ierr);
79   ierr = VecDestroy(&mat_ctx->lambda_local);CHKERRQ(ierr);
80   ierr = VecDestroy(&mat_ctx->temp_solution_D);CHKERRQ(ierr);
81   ierr = VecDestroy(&mat_ctx->temp_solution_B);CHKERRQ(ierr);
82   ierr = MatDestroy(&mat_ctx->B_delta);CHKERRQ(ierr);
83   ierr = MatDestroy(&mat_ctx->B_Ddelta);CHKERRQ(ierr);
84   ierr = MatDestroy(&mat_ctx->B_BB);CHKERRQ(ierr);
85   ierr = MatDestroy(&mat_ctx->B_BI);CHKERRQ(ierr);
86   ierr = MatDestroy(&mat_ctx->Bt_BB);CHKERRQ(ierr);
87   ierr = MatDestroy(&mat_ctx->Bt_BI);CHKERRQ(ierr);
88   ierr = MatDestroy(&mat_ctx->C);CHKERRQ(ierr);
89   ierr = VecDestroy(&mat_ctx->rhs_flip);CHKERRQ(ierr);
90   ierr = VecDestroy(&mat_ctx->vP);CHKERRQ(ierr);
91   ierr = VecDestroy(&mat_ctx->xPg);CHKERRQ(ierr);
92   ierr = VecDestroy(&mat_ctx->yPg);CHKERRQ(ierr);
93   ierr = VecScatterDestroy(&mat_ctx->l2g_lambda);CHKERRQ(ierr);
94   ierr = VecScatterDestroy(&mat_ctx->l2g_p);CHKERRQ(ierr);
95   ierr = VecScatterDestroy(&mat_ctx->g2g_p);CHKERRQ(ierr);
96   ierr = PCDestroy(&mat_ctx->pc);CHKERRQ(ierr); /* decrease PCBDDC reference count */
97   ierr = PetscFree(mat_ctx);CHKERRQ(ierr);
98   PetscFunctionReturn(0);
99 }
100 
101 PetscErrorCode PCBDDCDestroyFETIDPPC(PC pc)
102 {
103   FETIDPPC_ctx   pc_ctx;
104   PetscErrorCode ierr;
105 
106   PetscFunctionBegin;
107   ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr);
108   ierr = VecDestroy(&pc_ctx->lambda_local);CHKERRQ(ierr);
109   ierr = MatDestroy(&pc_ctx->B_Ddelta);CHKERRQ(ierr);
110   ierr = VecScatterDestroy(&pc_ctx->l2g_lambda);CHKERRQ(ierr);
111   ierr = MatDestroy(&pc_ctx->S_j);CHKERRQ(ierr);
112   ierr = PCDestroy(&pc_ctx->pc);CHKERRQ(ierr); /* decrease PCBDDC reference count */
113   ierr = KSPDestroy(&pc_ctx->kP);CHKERRQ(ierr);
114   ierr = VecDestroy(&pc_ctx->xPg);CHKERRQ(ierr);
115   ierr = VecDestroy(&pc_ctx->yPg);CHKERRQ(ierr);
116   ierr = PetscFree(pc_ctx);CHKERRQ(ierr);
117   PetscFunctionReturn(0);
118 }
119 
120 PetscErrorCode PCBDDCSetupFETIDPMatContext(FETIDPMat_ctx fetidpmat_ctx )
121 {
122   PetscErrorCode ierr;
123   PC_IS          *pcis=(PC_IS*)fetidpmat_ctx->pc->data;
124   PC_BDDC        *pcbddc=(PC_BDDC*)fetidpmat_ctx->pc->data;
125   PCBDDCGraph    mat_graph=pcbddc->mat_graph;
126 #if defined(PETSC_USE_DEBUG)
127   Mat_IS         *matis  = (Mat_IS*)fetidpmat_ctx->pc->pmat->data;
128 #endif
129   MPI_Comm       comm;
130   Mat            ScalingMat,BD1,BD2;
131   Vec            fetidp_global;
132   IS             IS_l2g_lambda;
133   IS             subset,subset_mult,subset_n,isvert;
134   PetscBool      skip_node,fully_redundant;
135   PetscInt       i,j,k,s,n_boundary_dofs,n_global_lambda,n_vertices,partial_sum;
136   PetscInt       cum,n_local_lambda,n_lambda_for_dof,dual_size,n_neg_values,n_pos_values;
137   PetscMPIInt    rank,size,buf_size,neigh;
138   PetscScalar    scalar_value;
139   const PetscInt *vertex_indices;
140   PetscInt       *dual_dofs_boundary_indices,*aux_local_numbering_1;
141   const PetscInt *aux_global_numbering;
142   PetscInt       *aux_sums,*cols_B_delta,*l2g_indices;
143   PetscScalar    *array,*scaling_factors,*vals_B_delta;
144   PetscScalar    **all_factors;
145   PetscInt       *aux_local_numbering_2;
146   PetscLayout    llay;
147 
148   /* saddlepoint */
149   ISLocalToGlobalMapping l2gmap_p;
150   PetscLayout            play;
151   IS                     gP,pP;
152   PetscInt               nPl,nPg,nPgl;
153 
154   PetscFunctionBegin;
155   ierr = PetscObjectGetComm((PetscObject)(fetidpmat_ctx->pc),&comm);CHKERRQ(ierr);
156   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
157   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
158 
159   /* saddlepoint */
160   nPl      = 0;
161   nPg      = 0;
162   nPgl     = 0;
163   gP       = NULL;
164   pP       = NULL;
165   l2gmap_p = NULL;
166   play     = NULL;
167   ierr = PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_pP",(PetscObject*)&pP);CHKERRQ(ierr);
168   if (pP) { /* saddle point */
169     /* subdomain pressures in global numbering */
170     ierr = PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_gP",(PetscObject*)&gP);CHKERRQ(ierr);
171     if (!gP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"gP not present");
172     ierr = ISGetLocalSize(gP,&nPl);CHKERRQ(ierr);
173     ierr = VecCreate(PETSC_COMM_SELF,&fetidpmat_ctx->vP);CHKERRQ(ierr);
174     ierr = VecSetSizes(fetidpmat_ctx->vP,nPl,nPl);CHKERRQ(ierr);
175     ierr = VecSetType(fetidpmat_ctx->vP,VECSTANDARD);CHKERRQ(ierr);
176     ierr = VecSetUp(fetidpmat_ctx->vP);CHKERRQ(ierr);
177 
178     /* interface pressure matrix */
179     ierr = PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_C",(PetscObject*)&fetidpmat_ctx->C);CHKERRQ(ierr);
180     if (!fetidpmat_ctx->C) { /* null pressure block, compute layout and global numbering for interface pressures */
181       IS Pg;
182 
183       ierr = ISRenumber(gP,NULL,&nPg,&Pg);CHKERRQ(ierr);
184       ierr = ISLocalToGlobalMappingCreateIS(Pg,&l2gmap_p);CHKERRQ(ierr);
185       ierr = ISDestroy(&Pg);CHKERRQ(ierr);
186       ierr = PetscLayoutCreate(comm,&play);CHKERRQ(ierr);
187       ierr = PetscLayoutSetBlockSize(play,1);CHKERRQ(ierr);
188       ierr = PetscLayoutSetSize(play,nPg);CHKERRQ(ierr);
189       ierr = ISGetLocalSize(pP,&nPgl);CHKERRQ(ierr);
190       ierr = PetscLayoutSetLocalSize(play,nPgl);CHKERRQ(ierr);
191       ierr = PetscLayoutSetUp(play);CHKERRQ(ierr);
192     } else {
193       ierr = PetscObjectReference((PetscObject)fetidpmat_ctx->C);CHKERRQ(ierr);
194       ierr = MatGetLocalToGlobalMapping(fetidpmat_ctx->C,&l2gmap_p,NULL);CHKERRQ(ierr);
195       ierr = PetscObjectReference((PetscObject)l2gmap_p);CHKERRQ(ierr);
196       ierr = MatGetSize(fetidpmat_ctx->C,&nPg,NULL);CHKERRQ(ierr);
197       ierr = MatGetLocalSize(fetidpmat_ctx->C,NULL,&nPgl);CHKERRQ(ierr);
198       ierr = MatGetLayouts(fetidpmat_ctx->C,NULL,&llay);CHKERRQ(ierr);
199       ierr = PetscLayoutReference(llay,&play);CHKERRQ(ierr);
200     }
201     ierr = VecCreateMPIWithArray(comm,1,nPgl,nPg,NULL,&fetidpmat_ctx->xPg);CHKERRQ(ierr);
202     ierr = VecCreateMPIWithArray(comm,1,nPgl,nPg,NULL,&fetidpmat_ctx->yPg);CHKERRQ(ierr);
203 
204     /* import matrices for interface pressures coupling */
205     ierr = PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_B_BI",(PetscObject*)&fetidpmat_ctx->B_BI);CHKERRQ(ierr);
206     if (!fetidpmat_ctx->B_BI) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B_BI not present");
207     ierr = PetscObjectReference((PetscObject)fetidpmat_ctx->B_BI);CHKERRQ(ierr);
208 
209     ierr = PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_B_BB",(PetscObject*)&fetidpmat_ctx->B_BB);CHKERRQ(ierr);
210     if (!fetidpmat_ctx->B_BB) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B_BB not present");
211     ierr = PetscObjectReference((PetscObject)fetidpmat_ctx->B_BB);CHKERRQ(ierr);
212 
213     ierr = PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_Bt_BI",(PetscObject*)&fetidpmat_ctx->Bt_BI);CHKERRQ(ierr);
214     if (!fetidpmat_ctx->Bt_BI) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bt_BI not present");
215     ierr = PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BI);CHKERRQ(ierr);
216 
217     ierr = PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_Bt_BB",(PetscObject*)&fetidpmat_ctx->Bt_BB);CHKERRQ(ierr);
218     if (!fetidpmat_ctx->Bt_BB) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bt_BB not present");
219     ierr = PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BB);CHKERRQ(ierr);
220 
221     ierr = PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_flip" ,(PetscObject*)&fetidpmat_ctx->rhs_flip);CHKERRQ(ierr);
222     if (fetidpmat_ctx->rhs_flip) {
223       ierr = PetscObjectReference((PetscObject)fetidpmat_ctx->rhs_flip);CHKERRQ(ierr);
224     }
225   }
226 
227   /* Default type of lagrange multipliers is non-redundant */
228   fully_redundant = fetidpmat_ctx->fully_redundant;
229 
230   /* Evaluate local and global number of lagrange multipliers */
231   ierr = VecSet(pcis->vec1_N,0.0);CHKERRQ(ierr);
232   n_local_lambda = 0;
233   partial_sum = 0;
234   n_boundary_dofs = 0;
235   s = 0;
236 
237   /* Get Vertices used to define the BDDC */
238   ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr);
239   ierr = ISGetLocalSize(isvert,&n_vertices);CHKERRQ(ierr);
240   ierr = ISGetIndices(isvert,&vertex_indices);CHKERRQ(ierr);
241 
242   dual_size = pcis->n_B-n_vertices;
243   ierr = PetscMalloc1(dual_size,&dual_dofs_boundary_indices);CHKERRQ(ierr);
244   ierr = PetscMalloc1(dual_size,&aux_local_numbering_1);CHKERRQ(ierr);
245   ierr = PetscMalloc1(dual_size,&aux_local_numbering_2);CHKERRQ(ierr);
246 
247   ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
248   for (i=0;i<pcis->n;i++){
249     j = mat_graph->count[i]; /* RECALL: mat_graph->count[i] does not count myself */
250     if (j > 0) n_boundary_dofs++;
251     skip_node = PETSC_FALSE;
252     if (s < n_vertices && vertex_indices[s] == i) { /* it works for a sorted set of vertices */
253       skip_node = PETSC_TRUE;
254       s++;
255     }
256     if (j < 1) skip_node = PETSC_TRUE;
257     if (mat_graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) skip_node = PETSC_TRUE;
258     if (!skip_node) {
259       if (fully_redundant) {
260         /* fully redundant set of lagrange multipliers */
261         n_lambda_for_dof = (j*(j+1))/2;
262       } else {
263         n_lambda_for_dof = j;
264       }
265       n_local_lambda += j;
266       /* needed to evaluate global number of lagrange multipliers */
267       array[i]=(1.0*n_lambda_for_dof)/(j+1.0); /* already scaled for the next global sum */
268       /* store some data needed */
269       dual_dofs_boundary_indices[partial_sum] = n_boundary_dofs-1;
270       aux_local_numbering_1[partial_sum] = i;
271       aux_local_numbering_2[partial_sum] = n_lambda_for_dof;
272       partial_sum++;
273     }
274   }
275   ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
276   ierr = ISRestoreIndices(isvert,&vertex_indices);CHKERRQ(ierr);
277   ierr = PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr);
278   dual_size = partial_sum;
279 
280   /* compute global ordering of lagrange multipliers and associate l2g map */
281   ierr = ISCreateGeneral(comm,partial_sum,aux_local_numbering_1,PETSC_COPY_VALUES,&subset_n);CHKERRQ(ierr);
282   ierr = ISLocalToGlobalMappingApplyIS(pcis->mapping,subset_n,&subset);CHKERRQ(ierr);
283   ierr = ISDestroy(&subset_n);CHKERRQ(ierr);
284   ierr = ISCreateGeneral(comm,partial_sum,aux_local_numbering_2,PETSC_OWN_POINTER,&subset_mult);CHKERRQ(ierr);
285   ierr = ISRenumber(subset,subset_mult,&fetidpmat_ctx->n_lambda,&subset_n);CHKERRQ(ierr);
286   ierr = ISDestroy(&subset);CHKERRQ(ierr);
287 
288 #if defined(PETSC_USE_DEBUG)
289   ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
290   ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
291   ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
292   ierr = VecSum(pcis->vec1_global,&scalar_value);CHKERRQ(ierr);
293   i = (PetscInt)PetscRealPart(scalar_value);
294   if (i != fetidpmat_ctx->n_lambda) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Global number of multipliers mismatch! (%d!=%d)\n",fetidpmat_ctx->n_lambda,i);
295 #endif
296 
297   /* init data for scaling factors exchange */
298   if (!pcbddc->use_deluxe_scaling) {
299     PetscInt    *ptrs_buffer,neigh_position;
300     PetscScalar *send_buffer,*recv_buffer;
301     MPI_Request *send_reqs,*recv_reqs;
302 
303     partial_sum = 0;
304     ierr = PetscMalloc1(pcis->n_neigh,&ptrs_buffer);CHKERRQ(ierr);
305     ierr = PetscMalloc1(PetscMax(pcis->n_neigh-1,0),&send_reqs);CHKERRQ(ierr);
306     ierr = PetscMalloc1(PetscMax(pcis->n_neigh-1,0),&recv_reqs);CHKERRQ(ierr);
307     ierr = PetscMalloc1(pcis->n+1,&all_factors);CHKERRQ(ierr);
308     if (pcis->n_neigh > 0) ptrs_buffer[0]=0;
309     for (i=1;i<pcis->n_neigh;i++) {
310       partial_sum += pcis->n_shared[i];
311       ptrs_buffer[i] = ptrs_buffer[i-1]+pcis->n_shared[i];
312     }
313     ierr = PetscMalloc1(partial_sum,&send_buffer);CHKERRQ(ierr);
314     ierr = PetscMalloc1(partial_sum,&recv_buffer);CHKERRQ(ierr);
315     ierr = PetscMalloc1(partial_sum,&all_factors[0]);CHKERRQ(ierr);
316     for (i=0;i<pcis->n-1;i++) {
317       j = mat_graph->count[i];
318       all_factors[i+1]=all_factors[i]+j;
319     }
320 
321     /* scatter B scaling to N vec */
322     ierr = VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
323     ierr = VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
324     /* communications */
325     ierr = VecGetArrayRead(pcis->vec1_N,(const PetscScalar**)&array);CHKERRQ(ierr);
326     for (i=1;i<pcis->n_neigh;i++) {
327       for (j=0;j<pcis->n_shared[i];j++) {
328         send_buffer[ptrs_buffer[i-1]+j]=array[pcis->shared[i][j]];
329       }
330       ierr = PetscMPIIntCast(ptrs_buffer[i]-ptrs_buffer[i-1],&buf_size);CHKERRQ(ierr);
331       ierr = PetscMPIIntCast(pcis->neigh[i],&neigh);CHKERRQ(ierr);
332       ierr = MPI_Isend(&send_buffer[ptrs_buffer[i-1]],buf_size,MPIU_SCALAR,neigh,0,comm,&send_reqs[i-1]);CHKERRQ(ierr);
333       ierr = MPI_Irecv(&recv_buffer[ptrs_buffer[i-1]],buf_size,MPIU_SCALAR,neigh,0,comm,&recv_reqs[i-1]);CHKERRQ(ierr);
334     }
335     ierr = VecRestoreArrayRead(pcis->vec1_N,(const PetscScalar**)&array);CHKERRQ(ierr);
336     if (pcis->n_neigh > 0) {
337       ierr = MPI_Waitall(pcis->n_neigh-1,recv_reqs,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
338     }
339     /* put values in correct places */
340     for (i=1;i<pcis->n_neigh;i++) {
341       for (j=0;j<pcis->n_shared[i];j++) {
342         k = pcis->shared[i][j];
343         neigh_position = 0;
344         while(mat_graph->neighbours_set[k][neigh_position] != pcis->neigh[i]) {neigh_position++;}
345         all_factors[k][neigh_position]=recv_buffer[ptrs_buffer[i-1]+j];
346       }
347     }
348     if (pcis->n_neigh > 0) {
349       ierr = MPI_Waitall(pcis->n_neigh-1,send_reqs,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
350     }
351     ierr = PetscFree(send_reqs);CHKERRQ(ierr);
352     ierr = PetscFree(recv_reqs);CHKERRQ(ierr);
353     ierr = PetscFree(send_buffer);CHKERRQ(ierr);
354     ierr = PetscFree(recv_buffer);CHKERRQ(ierr);
355     ierr = PetscFree(ptrs_buffer);CHKERRQ(ierr);
356   }
357 
358   /* Compute B and B_delta (local actions) */
359   ierr = PetscMalloc1(pcis->n_neigh,&aux_sums);CHKERRQ(ierr);
360   ierr = PetscMalloc1(n_local_lambda,&l2g_indices);CHKERRQ(ierr);
361   ierr = PetscMalloc1(n_local_lambda,&vals_B_delta);CHKERRQ(ierr);
362   ierr = PetscMalloc1(n_local_lambda,&cols_B_delta);CHKERRQ(ierr);
363   if (!pcbddc->use_deluxe_scaling) {
364     ierr = PetscMalloc1(n_local_lambda,&scaling_factors);CHKERRQ(ierr);
365   } else {
366     scaling_factors = NULL;
367     all_factors     = NULL;
368   }
369   ierr = ISGetIndices(subset_n,&aux_global_numbering);CHKERRQ(ierr);
370   partial_sum=0;
371   cum = 0;
372   for (i=0;i<dual_size;i++) {
373     n_global_lambda = aux_global_numbering[cum];
374     j = mat_graph->count[aux_local_numbering_1[i]];
375     aux_sums[0]=0;
376     for (s=1;s<j;s++) {
377       aux_sums[s]=aux_sums[s-1]+j-s+1;
378     }
379     if (all_factors) array = all_factors[aux_local_numbering_1[i]];
380     n_neg_values = 0;
381     while(n_neg_values < j && mat_graph->neighbours_set[aux_local_numbering_1[i]][n_neg_values] < rank) {n_neg_values++;}
382     n_pos_values = j - n_neg_values;
383     if (fully_redundant) {
384       for (s=0;s<n_neg_values;s++) {
385         l2g_indices    [partial_sum+s]=aux_sums[s]+n_neg_values-s-1+n_global_lambda;
386         cols_B_delta   [partial_sum+s]=dual_dofs_boundary_indices[i];
387         vals_B_delta   [partial_sum+s]=-1.0;
388         if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum+s]=array[s];
389       }
390       for (s=0;s<n_pos_values;s++) {
391         l2g_indices    [partial_sum+s+n_neg_values]=aux_sums[n_neg_values]+s+n_global_lambda;
392         cols_B_delta   [partial_sum+s+n_neg_values]=dual_dofs_boundary_indices[i];
393         vals_B_delta   [partial_sum+s+n_neg_values]=1.0;
394         if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum+s+n_neg_values]=array[s+n_neg_values];
395       }
396       partial_sum += j;
397     } else {
398       /* l2g_indices and default cols and vals of B_delta */
399       for (s=0;s<j;s++) {
400         l2g_indices    [partial_sum+s]=n_global_lambda+s;
401         cols_B_delta   [partial_sum+s]=dual_dofs_boundary_indices[i];
402         vals_B_delta   [partial_sum+s]=0.0;
403       }
404       /* B_delta */
405       if ( n_neg_values > 0 ) { /* there's a rank next to me to the left */
406         vals_B_delta   [partial_sum+n_neg_values-1]=-1.0;
407       }
408       if ( n_neg_values < j ) { /* there's a rank next to me to the right */
409         vals_B_delta   [partial_sum+n_neg_values]=1.0;
410       }
411       /* scaling as in Klawonn-Widlund 1999 */
412       if (!pcbddc->use_deluxe_scaling) {
413         for (s=0;s<n_neg_values;s++) {
414           scalar_value = 0.0;
415           for (k=0;k<s+1;k++) scalar_value += array[k];
416           scaling_factors[partial_sum+s] = -scalar_value;
417         }
418         for (s=0;s<n_pos_values;s++) {
419           scalar_value = 0.0;
420           for (k=s+n_neg_values;k<j;k++) scalar_value += array[k];
421           scaling_factors[partial_sum+s+n_neg_values] = scalar_value;
422         }
423       }
424       partial_sum += j;
425     }
426     cum += aux_local_numbering_2[i];
427   }
428   ierr = ISRestoreIndices(subset_n,&aux_global_numbering);CHKERRQ(ierr);
429   ierr = ISDestroy(&subset_mult);CHKERRQ(ierr);
430   ierr = ISDestroy(&subset_n);CHKERRQ(ierr);
431   ierr = PetscFree(aux_sums);CHKERRQ(ierr);
432   ierr = PetscFree(aux_local_numbering_1);CHKERRQ(ierr);
433   ierr = PetscFree(dual_dofs_boundary_indices);CHKERRQ(ierr);
434   if (all_factors) {
435     ierr = PetscFree(all_factors[0]);CHKERRQ(ierr);
436     ierr = PetscFree(all_factors);CHKERRQ(ierr);
437   }
438 
439   /* Create local part of B_delta */
440   ierr = MatCreate(PETSC_COMM_SELF,&fetidpmat_ctx->B_delta);CHKERRQ(ierr);
441   ierr = MatSetSizes(fetidpmat_ctx->B_delta,n_local_lambda,pcis->n_B,n_local_lambda,pcis->n_B);CHKERRQ(ierr);
442   ierr = MatSetType(fetidpmat_ctx->B_delta,MATSEQAIJ);CHKERRQ(ierr);
443   ierr = MatSeqAIJSetPreallocation(fetidpmat_ctx->B_delta,1,NULL);CHKERRQ(ierr);
444   ierr = MatSetOption(fetidpmat_ctx->B_delta,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
445   for (i=0;i<n_local_lambda;i++) {
446     ierr = MatSetValue(fetidpmat_ctx->B_delta,i,cols_B_delta[i],vals_B_delta[i],INSERT_VALUES);CHKERRQ(ierr);
447   }
448   ierr = PetscFree(vals_B_delta);CHKERRQ(ierr);
449   ierr = MatAssemblyBegin(fetidpmat_ctx->B_delta,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
450   ierr = MatAssemblyEnd(fetidpmat_ctx->B_delta,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
451 
452   BD1 = NULL;
453   BD2 = NULL;
454   if (fully_redundant) {
455     if (pcbddc->use_deluxe_scaling) SETERRQ(comm,PETSC_ERR_SUP,"Deluxe FETIDP with fully-redundant multipliers to be implemented");
456     ierr = MatCreate(PETSC_COMM_SELF,&ScalingMat);CHKERRQ(ierr);
457     ierr = MatSetSizes(ScalingMat,n_local_lambda,n_local_lambda,n_local_lambda,n_local_lambda);CHKERRQ(ierr);
458     ierr = MatSetType(ScalingMat,MATSEQAIJ);CHKERRQ(ierr);
459     ierr = MatSeqAIJSetPreallocation(ScalingMat,1,NULL);CHKERRQ(ierr);
460     for (i=0;i<n_local_lambda;i++) {
461       ierr = MatSetValue(ScalingMat,i,i,scaling_factors[i],INSERT_VALUES);CHKERRQ(ierr);
462     }
463     ierr = MatAssemblyBegin(ScalingMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
464     ierr = MatAssemblyEnd(ScalingMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
465     ierr = MatMatMult(ScalingMat,fetidpmat_ctx->B_delta,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&fetidpmat_ctx->B_Ddelta);CHKERRQ(ierr);
466     ierr = MatDestroy(&ScalingMat);CHKERRQ(ierr);
467   } else {
468     ierr = MatCreate(PETSC_COMM_SELF,&fetidpmat_ctx->B_Ddelta);CHKERRQ(ierr);
469     ierr = MatSetSizes(fetidpmat_ctx->B_Ddelta,n_local_lambda,pcis->n_B,n_local_lambda,pcis->n_B);CHKERRQ(ierr);
470     if (!pcbddc->use_deluxe_scaling) {
471       ierr = MatSetType(fetidpmat_ctx->B_Ddelta,MATSEQAIJ);CHKERRQ(ierr);
472       ierr = MatSeqAIJSetPreallocation(fetidpmat_ctx->B_Ddelta,1,NULL);CHKERRQ(ierr);
473       for (i=0;i<n_local_lambda;i++) {
474         ierr = MatSetValue(fetidpmat_ctx->B_Ddelta,i,cols_B_delta[i],scaling_factors[i],INSERT_VALUES);CHKERRQ(ierr);
475       }
476       ierr = MatAssemblyBegin(fetidpmat_ctx->B_Ddelta,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
477       ierr = MatAssemblyEnd(fetidpmat_ctx->B_Ddelta,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
478     } else {
479       /* scaling as in Klawonn-Widlund 1999 */
480       PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
481       PCBDDCSubSchurs     sub_schurs = pcbddc->sub_schurs;
482       Mat                 T;
483       PetscScalar         *W,lwork,*Bwork;
484       const PetscInt      *idxs;
485       PetscInt            cum,mss,*nnz;
486       PetscBLASInt        *pivots,B_lwork,B_N,B_ierr;
487 
488       if (!pcbddc->deluxe_singlemat) SETERRQ(comm,PETSC_ERR_USER,"Cannot compute B_Ddelta! rerun with -pc_bddc_deluxe_singlemat");
489       if (deluxe_ctx->change) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot compute B_Ddelta with deluxe scaling with active change context");
490 
491       mss  = 0;
492       ierr = PetscCalloc1(pcis->n_B,&nnz);CHKERRQ(ierr);
493       if (sub_schurs->is_Ej_all) {
494         ierr = ISGetIndices(sub_schurs->is_Ej_all,&idxs);CHKERRQ(ierr);
495         for (i=0,cum=0;i<sub_schurs->n_subs;i++) {
496           PetscInt subset_size;
497 
498           ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
499           for (j=0;j<subset_size;j++) nnz[idxs[j+cum]] = subset_size;
500           mss  = PetscMax(mss,subset_size);
501           cum += subset_size;
502         }
503       }
504       ierr = MatCreate(PETSC_COMM_SELF,&T);CHKERRQ(ierr);
505       ierr = MatSetSizes(T,pcis->n_B,pcis->n_B,pcis->n_B,pcis->n_B);CHKERRQ(ierr);
506       ierr = MatSetType(T,MATSEQAIJ);CHKERRQ(ierr);
507       ierr = MatSeqAIJSetPreallocation(T,0,nnz);CHKERRQ(ierr);
508       ierr = PetscFree(nnz);CHKERRQ(ierr);
509 
510       /* workspace allocation */
511       B_lwork = 0;
512       if (mss) {
513         PetscScalar dummy = 1;
514 
515         B_lwork = -1;
516         ierr = PetscBLASIntCast(mss,&B_N);CHKERRQ(ierr);
517         ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
518         PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,&dummy,&B_N,&B_N,&lwork,&B_lwork,&B_ierr));
519         ierr = PetscFPTrapPop();CHKERRQ(ierr);
520         if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GETRI Lapack routine %d",(int)B_ierr);
521         ierr = PetscBLASIntCast((PetscInt)PetscRealPart(lwork),&B_lwork);CHKERRQ(ierr);
522       }
523       ierr = PetscMalloc3(mss*mss,&W,mss,&pivots,B_lwork,&Bwork);CHKERRQ(ierr);
524 
525       for (i=0,cum=0;i<sub_schurs->n_subs;i++) {
526         PetscScalar  *M;
527         PetscInt     subset_size;
528 
529         ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
530         ierr = PetscBLASIntCast(subset_size,&B_N);CHKERRQ(ierr);
531         ierr = MatDenseGetArray(deluxe_ctx->seq_mat[i],&M);CHKERRQ(ierr);
532         ierr = PetscMemcpy(W,M,subset_size*subset_size*sizeof(PetscScalar));CHKERRQ(ierr);
533         ierr = MatDenseRestoreArray(deluxe_ctx->seq_mat[i],&M);CHKERRQ(ierr);
534         ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
535         PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,W,&B_N,pivots,&B_ierr));
536         if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
537         PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,W,&B_N,pivots,Bwork,&B_lwork,&B_ierr));
538         if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
539         ierr = PetscFPTrapPop();CHKERRQ(ierr);
540         ierr = MatSetValues(T,subset_size,idxs+cum,subset_size,idxs+cum,W,INSERT_VALUES);CHKERRQ(ierr);
541         cum += subset_size;
542       }
543       ierr = MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
544       ierr = MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
545       ierr = MatMatTransposeMult(T,fetidpmat_ctx->B_delta,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&BD1);CHKERRQ(ierr);
546       ierr = MatMatMult(fetidpmat_ctx->B_delta,BD1,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&BD2);CHKERRQ(ierr);
547       ierr = MatDestroy(&T);CHKERRQ(ierr);
548       ierr = PetscFree3(W,pivots,Bwork);CHKERRQ(ierr);
549       if (sub_schurs->is_Ej_all) {
550         ierr = ISRestoreIndices(sub_schurs->is_Ej_all,&idxs);CHKERRQ(ierr);
551       }
552     }
553   }
554   ierr = PetscFree(scaling_factors);CHKERRQ(ierr);
555   ierr = PetscFree(cols_B_delta);CHKERRQ(ierr);
556 
557   /* Layout of multipliers */
558   ierr = PetscLayoutCreate(comm,&llay);CHKERRQ(ierr);
559   ierr = PetscLayoutSetBlockSize(llay,1);CHKERRQ(ierr);
560   ierr = PetscLayoutSetSize(llay,fetidpmat_ctx->n_lambda);CHKERRQ(ierr);
561   ierr = PetscLayoutSetUp(llay);CHKERRQ(ierr);
562   ierr = PetscLayoutGetLocalSize(llay,&fetidpmat_ctx->n);CHKERRQ(ierr);
563 
564   /* Local work vector of multipliers */
565   ierr = VecCreate(PETSC_COMM_SELF,&fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
566   ierr = VecSetSizes(fetidpmat_ctx->lambda_local,n_local_lambda,n_local_lambda);CHKERRQ(ierr);
567   ierr = VecSetType(fetidpmat_ctx->lambda_local,VECSEQ);CHKERRQ(ierr);
568 
569   if (BD2) {
570     ISLocalToGlobalMapping l2g;
571     Mat                    T,TA,*pT;
572     IS                     is;
573     PetscInt               nl,N;
574     BDdelta_DN             ctx;
575 
576     ierr = PetscLayoutGetLocalSize(llay,&nl);CHKERRQ(ierr);
577     ierr = PetscLayoutGetSize(llay,&N);CHKERRQ(ierr);
578     ierr = MatCreate(comm,&T);CHKERRQ(ierr);
579     ierr = MatSetSizes(T,nl,nl,N,N);CHKERRQ(ierr);
580     ierr = MatSetType(T,MATIS);CHKERRQ(ierr);
581     ierr = ISLocalToGlobalMappingCreate(comm,1,n_local_lambda,l2g_indices,PETSC_COPY_VALUES,&l2g);CHKERRQ(ierr);
582     ierr = MatSetLocalToGlobalMapping(T,l2g,l2g);CHKERRQ(ierr);
583     ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr);
584     ierr = MatISSetLocalMat(T,BD2);CHKERRQ(ierr);
585     ierr = MatDestroy(&BD2);CHKERRQ(ierr);
586     ierr = MatISGetMPIXAIJ(T,MAT_INITIAL_MATRIX,&TA);CHKERRQ(ierr);
587     ierr = MatDestroy(&T);CHKERRQ(ierr);
588     ierr = ISCreateGeneral(comm,n_local_lambda,l2g_indices,PETSC_USE_POINTER,&is);CHKERRQ(ierr);
589     ierr = MatCreateSubMatrices(TA,1,&is,&is,MAT_INITIAL_MATRIX,&pT);CHKERRQ(ierr);
590     ierr = MatDestroy(&TA);CHKERRQ(ierr);
591     ierr = ISDestroy(&is);CHKERRQ(ierr);
592     BD2  = pT[0];
593     ierr = PetscFree(pT);CHKERRQ(ierr);
594 
595     /* B_Ddelta for non-redundant multipliers with deluxe scaling */
596     ierr = PetscNew(&ctx);CHKERRQ(ierr);
597     ierr = MatSetType(fetidpmat_ctx->B_Ddelta,MATSHELL);CHKERRQ(ierr);
598     ierr = MatShellSetContext(fetidpmat_ctx->B_Ddelta,(void *)ctx);CHKERRQ(ierr);
599     ierr = MatShellSetOperation(fetidpmat_ctx->B_Ddelta,MATOP_MULT,(void (*)(void))MatMult_BDdelta_deluxe_nonred);CHKERRQ(ierr);
600     ierr = MatShellSetOperation(fetidpmat_ctx->B_Ddelta,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_BDdelta_deluxe_nonred);CHKERRQ(ierr);
601     ierr = MatShellSetOperation(fetidpmat_ctx->B_Ddelta,MATOP_DESTROY,(void (*)(void))MatDestroy_BDdelta_deluxe_nonred);CHKERRQ(ierr);
602     ierr = MatSetUp(fetidpmat_ctx->B_Ddelta);CHKERRQ(ierr);
603 
604     ierr = PetscObjectReference((PetscObject)BD1);CHKERRQ(ierr);
605     ctx->BD = BD1;
606     ierr = KSPCreate(PETSC_COMM_SELF,&ctx->kBD);CHKERRQ(ierr);
607     ierr = KSPSetOperators(ctx->kBD,BD2,BD2);CHKERRQ(ierr);
608     ierr = VecDuplicate(fetidpmat_ctx->lambda_local,&ctx->work);CHKERRQ(ierr);
609     fetidpmat_ctx->deluxe_nonred = PETSC_TRUE;
610   }
611   ierr = MatDestroy(&BD1);CHKERRQ(ierr);
612   ierr = MatDestroy(&BD2);CHKERRQ(ierr);
613 
614   /* fetidpmat sizes */
615   fetidpmat_ctx->n += nPgl;
616   fetidpmat_ctx->N  = fetidpmat_ctx->n_lambda+nPg;
617 
618   /* Global vector for FETI-DP linear system */
619   ierr = VecCreate(comm,&fetidp_global);CHKERRQ(ierr);
620   ierr = VecSetSizes(fetidp_global,fetidpmat_ctx->n,fetidpmat_ctx->N);CHKERRQ(ierr);
621   ierr = VecSetType(fetidp_global,VECMPI);CHKERRQ(ierr);
622   ierr = VecSetUp(fetidp_global);CHKERRQ(ierr);
623 
624   /* Decide layout for fetidp dofs: if it is a saddle point problem
625      pressure is ordered first in the local part of the global vector
626      of the FETI-DP linear system */
627   if (nPg) {
628     IS             IS_l2g_p,ais;
629     PetscLayout    alay;
630     const PetscInt *idxs,*pranges,*aranges,*lranges;
631     PetscInt       *l2g_indices_p,rst;
632 
633     ierr = PetscMalloc1(nPl,&l2g_indices_p);CHKERRQ(ierr);
634     ierr = VecGetLayout(fetidp_global,&alay);CHKERRQ(ierr);
635     ierr = PetscLayoutGetRanges(alay,&aranges);CHKERRQ(ierr);
636     ierr = PetscLayoutGetRanges(play,&pranges);CHKERRQ(ierr);
637     ierr = PetscLayoutGetRanges(llay,&lranges);CHKERRQ(ierr);
638     ierr = ISLocalToGlobalMappingGetIndices(l2gmap_p,&idxs);CHKERRQ(ierr);
639     /* shift local to global indices for pressure */
640     for (i=0;i<nPl;i++) {
641       PetscInt owner;
642 
643       ierr = PetscLayoutFindOwner(play,idxs[i],&owner);CHKERRQ(ierr);
644       l2g_indices_p[i] = idxs[i]-pranges[owner]+aranges[owner];
645     }
646     ierr = ISLocalToGlobalMappingRestoreIndices(l2gmap_p,&idxs);CHKERRQ(ierr);
647     ierr = ISCreateGeneral(comm,nPl,l2g_indices_p,PETSC_OWN_POINTER,&IS_l2g_p);CHKERRQ(ierr);
648 
649     /* local to global scatter for interface pressure */
650     ierr = VecScatterCreate(fetidpmat_ctx->vP,NULL,fetidp_global,IS_l2g_p,&fetidpmat_ctx->l2g_p);CHKERRQ(ierr);
651     ierr = ISDestroy(&IS_l2g_p);CHKERRQ(ierr);
652 
653     /* shift local to global indices for multipliers */
654     for (i=0;i<n_local_lambda;i++) {
655       PetscInt owner,ps;
656 
657       ierr = PetscLayoutFindOwner(llay,l2g_indices[i],&owner);CHKERRQ(ierr);
658       ps = pranges[owner+1]-pranges[owner];
659       l2g_indices[i] = l2g_indices[i]-lranges[owner]+aranges[owner]+ps;
660     }
661 
662     /* scatter from alldofs to interface pressures global fetidp vector */
663     ierr = PetscLayoutGetRange(alay,&rst,NULL);CHKERRQ(ierr);
664     ierr = ISCreateStride(comm,nPgl,rst,1,&ais);CHKERRQ(ierr);
665     ierr = VecScatterCreate(pcis->vec1_global,pP,fetidp_global,ais,&fetidpmat_ctx->g2g_p);CHKERRQ(ierr);
666     ierr = ISDestroy(&ais);CHKERRQ(ierr);
667   }
668   ierr = PetscLayoutDestroy(&llay);CHKERRQ(ierr);
669   ierr = PetscLayoutDestroy(&play);CHKERRQ(ierr);
670   ierr = ISCreateGeneral(comm,n_local_lambda,l2g_indices,PETSC_OWN_POINTER,&IS_l2g_lambda);CHKERRQ(ierr);
671 
672   /* scatter from local to global multipliers */
673   ierr = VecScatterCreate(fetidpmat_ctx->lambda_local,NULL,fetidp_global,IS_l2g_lambda,&fetidpmat_ctx->l2g_lambda);CHKERRQ(ierr);
674   ierr = ISDestroy(&IS_l2g_lambda);CHKERRQ(ierr);
675   ierr = ISLocalToGlobalMappingDestroy(&l2gmap_p);CHKERRQ(ierr);
676   ierr = VecDestroy(&fetidp_global);CHKERRQ(ierr);
677 
678   /* Create some work vectors needed by fetidp */
679   ierr = VecDuplicate(pcis->vec1_B,&fetidpmat_ctx->temp_solution_B);CHKERRQ(ierr);
680   ierr = VecDuplicate(pcis->vec1_D,&fetidpmat_ctx->temp_solution_D);CHKERRQ(ierr);
681   PetscFunctionReturn(0);
682 }
683 
684 
685 PetscErrorCode PCBDDCSetupFETIDPPCContext(Mat fetimat, FETIDPPC_ctx fetidppc_ctx)
686 {
687   FETIDPMat_ctx  mat_ctx;
688   PC_BDDC        *pcbddc = (PC_BDDC*)fetidppc_ctx->pc->data;
689   PC_IS          *pcis = (PC_IS*)fetidppc_ctx->pc->data;
690   PetscBool      lumped = PETSC_FALSE;
691   PetscErrorCode ierr;
692 
693   PetscFunctionBegin;
694   ierr = MatShellGetContext(fetimat,(void**)&mat_ctx);CHKERRQ(ierr);
695   /* get references from objects created when setting up feti mat context */
696   ierr = PetscObjectReference((PetscObject)mat_ctx->lambda_local);CHKERRQ(ierr);
697   fetidppc_ctx->lambda_local = mat_ctx->lambda_local;
698   ierr = PetscObjectReference((PetscObject)mat_ctx->B_Ddelta);CHKERRQ(ierr);
699   fetidppc_ctx->B_Ddelta = mat_ctx->B_Ddelta;
700   if (mat_ctx->deluxe_nonred) {
701     PC               pc,mpc;
702     BDdelta_DN       ctx;
703     MatSolverPackage solver;
704     const char       *prefix;
705 
706     ierr = MatShellGetContext(mat_ctx->B_Ddelta,&ctx);CHKERRQ(ierr);
707     ierr = KSPSetType(ctx->kBD,KSPPREONLY);CHKERRQ(ierr);
708     ierr = KSPGetPC(ctx->kBD,&mpc);CHKERRQ(ierr);
709     ierr = KSPGetPC(pcbddc->ksp_D,&pc);CHKERRQ(ierr);
710     ierr = PCSetType(mpc,PCLU);CHKERRQ(ierr);
711     ierr = PCFactorGetMatSolverPackage(pc,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
712     if (solver) {
713       ierr = PCFactorSetMatSolverPackage(mpc,solver);CHKERRQ(ierr);
714     }
715     ierr = MatGetOptionsPrefix(fetimat,&prefix);CHKERRQ(ierr);
716     ierr = KSPSetOptionsPrefix(ctx->kBD,prefix);CHKERRQ(ierr);
717     ierr = KSPAppendOptionsPrefix(ctx->kBD,"bddelta_");CHKERRQ(ierr);
718     ierr = KSPSetFromOptions(ctx->kBD);CHKERRQ(ierr);
719   }
720 
721   ierr = PetscObjectReference((PetscObject)mat_ctx->l2g_lambda);CHKERRQ(ierr);
722   fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda;
723   /* Dirichlet preconditioner */
724   ierr = PetscOptionsGetBool(NULL,((PetscObject)fetimat)->prefix,"-pc_lumped",&lumped,NULL);CHKERRQ(ierr);
725   if (!lumped) {
726     IS        iP;
727     PetscBool discrete_harmonic = PETSC_FALSE;
728 
729     ierr = PetscObjectQuery((PetscObject)fetidppc_ctx->pc,"__KSPFETIDP_iP",(PetscObject*)&iP);CHKERRQ(ierr);
730     if (iP) {
731       ierr = PetscOptionsGetBool(NULL,((PetscObject)fetimat)->prefix,"-pc_discrete_harmonic",&discrete_harmonic,NULL);CHKERRQ(ierr);
732     }
733     if (discrete_harmonic) {
734       KSP        sksp;
735       PC         pc;
736       Mat        A_II,A_IB,A_BI;
737       IS         aB;
738       PetscInt   nb;
739       PetscBool  isshell;
740       KSPType    ksptype;
741       const char *prefix;
742 
743       /*
744         We constructs a Schur complement for
745 
746         | A_II A_ID |
747         | A_DI A_DD |
748 
749         instead of
750 
751         | A_II  B^t_II A_ID |
752         | B_II -C_II   B_ID |
753         | A_DI  B^t_ID A_DD |
754 
755       */
756       ierr = ISGetLocalSize(pcis->is_B_local,&nb);CHKERRQ(ierr);
757       ierr = ISCreateStride(PetscObjectComm((PetscObject)pcis->A_II),nb,0,1,&aB);CHKERRQ(ierr);
758       ierr = MatCreateSubMatrix(pcis->A_II,iP,iP,MAT_INITIAL_MATRIX,&A_II);CHKERRQ(ierr);
759       ierr = MatCreateSubMatrix(pcis->A_IB,iP,aB,MAT_INITIAL_MATRIX,&A_IB);CHKERRQ(ierr);
760       ierr = MatCreateSubMatrix(pcis->A_BI,aB,iP,MAT_INITIAL_MATRIX,&A_BI);CHKERRQ(ierr);
761       ierr = MatCreateSchurComplement(A_II,A_II,A_IB,A_BI,pcis->A_BB,&fetidppc_ctx->S_j);CHKERRQ(ierr);
762 
763       /* propagate settings of solver */
764       ierr = MatSchurComplementGetKSP(fetidppc_ctx->S_j,&sksp);CHKERRQ(ierr);
765       ierr = KSPGetType(pcis->ksp_D,&ksptype);CHKERRQ(ierr);
766       ierr = KSPSetType(sksp,ksptype);CHKERRQ(ierr);
767       ierr = KSPGetPC(pcis->ksp_D,&pc);CHKERRQ(ierr);
768       ierr = PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&isshell);CHKERRQ(ierr);
769       if (!isshell) {
770         MatSolverPackage solver;
771         PCType           pctype;
772 
773         ierr = PCGetType(pc,&pctype);CHKERRQ(ierr);
774         ierr = PCFactorGetMatSolverPackage(pc,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
775         ierr = KSPGetPC(sksp,&pc);CHKERRQ(ierr);
776         ierr = PCSetType(pc,pctype);CHKERRQ(ierr);
777         if (solver) {
778           ierr = PCFactorSetMatSolverPackage(pc,solver);CHKERRQ(ierr);
779         }
780       } else {
781         ierr = KSPGetPC(sksp,&pc);CHKERRQ(ierr);
782         ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
783       }
784       ierr = MatDestroy(&A_II);CHKERRQ(ierr);
785       ierr = MatDestroy(&A_IB);CHKERRQ(ierr);
786       ierr = MatDestroy(&A_BI);CHKERRQ(ierr);
787       ierr = ISDestroy(&aB);CHKERRQ(ierr);
788       ierr = MatGetOptionsPrefix(fetimat,&prefix);CHKERRQ(ierr);
789       ierr = KSPSetOptionsPrefix(sksp,prefix);CHKERRQ(ierr);
790       ierr = KSPAppendOptionsPrefix(sksp,"harmonic_");CHKERRQ(ierr);
791       ierr = KSPSetFromOptions(sksp);CHKERRQ(ierr);
792     } else { /* default Dirichlet preconditioner is pde-harmonic */
793       ierr = MatCreateSchurComplement(pcis->A_II,pcis->A_II,pcis->A_IB,pcis->A_BI,pcis->A_BB,&fetidppc_ctx->S_j);CHKERRQ(ierr);
794       ierr = MatSchurComplementSetKSP(fetidppc_ctx->S_j,pcis->ksp_D);CHKERRQ(ierr);
795     }
796   } else {
797     ierr = PetscObjectReference((PetscObject)pcis->A_BB);CHKERRQ(ierr);
798     fetidppc_ctx->S_j = pcis->A_BB;
799   }
800   /* saddle-point */
801   if (mat_ctx->xPg) {
802     ierr = PetscObjectReference((PetscObject)mat_ctx->xPg);CHKERRQ(ierr);
803     fetidppc_ctx->xPg = mat_ctx->xPg;
804     ierr = PetscObjectReference((PetscObject)mat_ctx->yPg);CHKERRQ(ierr);
805     fetidppc_ctx->yPg = mat_ctx->yPg;
806     ierr = PetscObjectQuery((PetscObject)fetidppc_ctx->pc,"__KSPFETIDP_PKSP",(PetscObject*)&fetidppc_ctx->kP);CHKERRQ(ierr);
807     ierr = PetscObjectReference((PetscObject)fetidppc_ctx->kP);CHKERRQ(ierr);
808   }
809   PetscFunctionReturn(0);
810 }
811 
812 PetscErrorCode FETIDPMatMult(Mat fetimat, Vec x, Vec y)
813 {
814   FETIDPMat_ctx  mat_ctx;
815   PC_BDDC        *pcbddc;
816   PC_IS          *pcis;
817   PetscErrorCode ierr;
818 
819   PetscFunctionBegin;
820   ierr = MatShellGetContext(fetimat,(void**)&mat_ctx);CHKERRQ(ierr);
821   pcis = (PC_IS*)mat_ctx->pc->data;
822   pcbddc = (PC_BDDC*)mat_ctx->pc->data;
823   /* Application of B_delta^T */
824   ierr = VecSet(pcis->vec1_B,0.);CHKERRQ(ierr);
825   ierr = VecScatterBegin(mat_ctx->l2g_lambda,x,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
826   ierr = VecScatterEnd(mat_ctx->l2g_lambda,x,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
827   ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
828 
829   /* Add contribution from saddle point */
830   if (mat_ctx->l2g_p) {
831     ierr = VecScatterBegin(mat_ctx->l2g_p,x,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
832     ierr = VecScatterEnd(mat_ctx->l2g_p,x,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
833     if (pcbddc->switch_static) {
834       ierr = MatMult(mat_ctx->Bt_BI,mat_ctx->vP,pcis->vec1_D);CHKERRQ(ierr);
835     }
836     ierr = MatMultAdd(mat_ctx->Bt_BB,mat_ctx->vP,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr);
837   } else {
838     if (pcbddc->switch_static) {
839       ierr = VecSet(pcis->vec1_D,0.0);CHKERRQ(ierr);
840     }
841   }
842   /* Application of \widetilde{S}^-1 */
843   ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr);
844   ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr);
845   ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr);
846   ierr = VecSet(y,0.0);CHKERRQ(ierr);
847   /* Application of B_delta */
848   ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr);
849   /* Contribution from boundary pressures */
850   if (mat_ctx->C) {
851     const PetscScalar *lx;
852     PetscScalar       *ly;
853 
854     /* pressures ordered first in x and y */
855     ierr = VecGetArrayRead(x,&lx);CHKERRQ(ierr);
856     ierr = VecGetArray(y,&ly);CHKERRQ(ierr);
857     ierr = VecPlaceArray(mat_ctx->xPg,lx);CHKERRQ(ierr);
858     ierr = VecPlaceArray(mat_ctx->yPg,ly);CHKERRQ(ierr);
859     ierr = MatMult(mat_ctx->C,mat_ctx->xPg,mat_ctx->yPg);CHKERRQ(ierr);
860     ierr = VecResetArray(mat_ctx->xPg);CHKERRQ(ierr);
861     ierr = VecResetArray(mat_ctx->yPg);CHKERRQ(ierr);
862     ierr = VecRestoreArrayRead(x,&lx);CHKERRQ(ierr);
863     ierr = VecRestoreArray(y,&ly);CHKERRQ(ierr);
864   }
865   /* Add contribution from saddle point */
866   if (mat_ctx->l2g_p) {
867     ierr = MatMult(mat_ctx->B_BB,pcis->vec1_B,mat_ctx->vP);CHKERRQ(ierr);
868     if (pcbddc->switch_static) {
869       ierr = MatMultAdd(mat_ctx->B_BI,pcis->vec1_D,mat_ctx->vP,mat_ctx->vP);CHKERRQ(ierr);
870     }
871     ierr = VecScatterBegin(mat_ctx->l2g_p,mat_ctx->vP,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
872     ierr = VecScatterEnd(mat_ctx->l2g_p,mat_ctx->vP,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
873   }
874   ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
875   ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
876   PetscFunctionReturn(0);
877 }
878 
879 PetscErrorCode FETIDPMatMultTranspose(Mat fetimat, Vec x, Vec y)
880 {
881   FETIDPMat_ctx  mat_ctx;
882   PC_IS          *pcis;
883   PetscErrorCode ierr;
884 
885   PetscFunctionBegin;
886   ierr = MatShellGetContext(fetimat,(void**)&mat_ctx);CHKERRQ(ierr);
887   pcis = (PC_IS*)mat_ctx->pc->data;
888   /* Application of B_delta^T */
889   ierr = VecScatterBegin(mat_ctx->l2g_lambda,x,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
890   ierr = VecScatterEnd(mat_ctx->l2g_lambda,x,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
891   ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
892   /* Application of \widetilde{S}^-1 */
893   ierr = VecSet(pcis->vec1_D,0.0);CHKERRQ(ierr);
894   ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_TRUE);CHKERRQ(ierr);
895   /* Application of B_delta */
896   ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr);
897   ierr = VecSet(y,0.0);CHKERRQ(ierr);
898   ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
899   ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
900   PetscFunctionReturn(0);
901 }
902 
903 PetscErrorCode FETIDPPCApply(PC fetipc, Vec x, Vec y)
904 {
905   FETIDPPC_ctx   pc_ctx;
906   PC_IS          *pcis;
907   PetscErrorCode ierr;
908 
909   PetscFunctionBegin;
910   ierr = PCShellGetContext(fetipc,(void**)&pc_ctx);CHKERRQ(ierr);
911   pcis = (PC_IS*)pc_ctx->pc->data;
912   /* Application of B_Ddelta^T */
913   ierr = VecScatterBegin(pc_ctx->l2g_lambda,x,pc_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
914   ierr = VecScatterEnd(pc_ctx->l2g_lambda,x,pc_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
915   ierr = VecSet(pcis->vec2_B,0.0);CHKERRQ(ierr);
916   ierr = MatMultTranspose(pc_ctx->B_Ddelta,pc_ctx->lambda_local,pcis->vec2_B);CHKERRQ(ierr);
917   /* Application of local Schur complement */
918   ierr = MatMult(pc_ctx->S_j,pcis->vec2_B,pcis->vec1_B);CHKERRQ(ierr);
919   /* Application of B_Ddelta */
920   ierr = MatMult(pc_ctx->B_Ddelta,pcis->vec1_B,pc_ctx->lambda_local);CHKERRQ(ierr);
921   ierr = VecSet(y,0.0);CHKERRQ(ierr);
922   ierr = VecScatterBegin(pc_ctx->l2g_lambda,pc_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
923   ierr = VecScatterEnd(pc_ctx->l2g_lambda,pc_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
924   /* interface pressure preconditioner */
925   if (pc_ctx->kP) {
926     const PetscScalar *lx;
927     PetscScalar *ly;
928 
929     /* pressures ordered first in x and y */
930     ierr = VecGetArrayRead(x,&lx);CHKERRQ(ierr);
931     ierr = VecGetArray(y,&ly);CHKERRQ(ierr);
932     ierr = VecPlaceArray(pc_ctx->xPg,lx);CHKERRQ(ierr);
933     ierr = VecPlaceArray(pc_ctx->yPg,ly);CHKERRQ(ierr);
934     ierr = KSPSolve(pc_ctx->kP,pc_ctx->xPg,pc_ctx->yPg);CHKERRQ(ierr);
935     ierr = VecResetArray(pc_ctx->xPg);CHKERRQ(ierr);
936     ierr = VecResetArray(pc_ctx->yPg);CHKERRQ(ierr);
937     ierr = VecRestoreArrayRead(x,&lx);CHKERRQ(ierr);
938     ierr = VecRestoreArray(y,&ly);CHKERRQ(ierr);
939   }
940   PetscFunctionReturn(0);
941 }
942 
943 PetscErrorCode FETIDPPCApplyTranspose(PC fetipc, Vec x, Vec y)
944 {
945   FETIDPPC_ctx   pc_ctx;
946   PC_IS          *pcis;
947   PetscErrorCode ierr;
948 
949   PetscFunctionBegin;
950   ierr = PCShellGetContext(fetipc,(void**)&pc_ctx);CHKERRQ(ierr);
951   pcis = (PC_IS*)pc_ctx->pc->data;
952   /* Application of B_Ddelta^T */
953   ierr = VecScatterBegin(pc_ctx->l2g_lambda,x,pc_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
954   ierr = VecScatterEnd(pc_ctx->l2g_lambda,x,pc_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
955   ierr = VecSet(pcis->vec2_B,0.0);CHKERRQ(ierr);
956   ierr = MatMultTranspose(pc_ctx->B_Ddelta,pc_ctx->lambda_local,pcis->vec2_B);CHKERRQ(ierr);
957   /* Application of local Schur complement */
958   ierr = MatMultTranspose(pc_ctx->S_j,pcis->vec2_B,pcis->vec1_B);CHKERRQ(ierr);
959   /* Application of B_Ddelta */
960   ierr = MatMult(pc_ctx->B_Ddelta,pcis->vec1_B,pc_ctx->lambda_local);CHKERRQ(ierr);
961   ierr = VecSet(y,0.0);CHKERRQ(ierr);
962   ierr = VecScatterBegin(pc_ctx->l2g_lambda,pc_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
963   ierr = VecScatterEnd(pc_ctx->l2g_lambda,pc_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
964   PetscFunctionReturn(0);
965 }
966 
967 PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer)
968 {
969   FETIDPPC_ctx      pc_ctx;
970   PetscBool         iascii;
971   PetscViewer       sviewer;
972   PetscErrorCode    ierr;
973 
974   PetscFunctionBegin;
975   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
976   if (iascii) {
977     PetscMPIInt rank;
978     PetscBool   isschur,isshell;
979 
980     ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr);
981     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRQ(ierr);
982     ierr = PetscObjectTypeCompare((PetscObject)pc_ctx->S_j,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
983     if (isschur) {
984       ierr = PetscViewerASCIIPrintf(viewer,"  FETI-DP multipliers Dirichlet preconditioner (just from rank 0)\n");CHKERRQ(ierr);
985     } else {
986       ierr = PetscViewerASCIIPrintf(viewer,"  FETI-DP multipliers Lumped preconditioner (just from rank 0)\n");CHKERRQ(ierr);
987     }
988     ierr = PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)pc_ctx->S_j),&sviewer);CHKERRQ(ierr);
989     if (!rank) {
990       ierr = PetscViewerPushFormat(sviewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
991       ierr = PetscViewerASCIIAddTab(sviewer,2);CHKERRQ(ierr);
992       ierr = MatView(pc_ctx->S_j,sviewer);CHKERRQ(ierr);
993       ierr = PetscViewerASCIISubtractTab(sviewer,2);CHKERRQ(ierr);
994       ierr = PetscViewerPopFormat(sviewer);CHKERRQ(ierr);
995     }
996     ierr = PetscObjectTypeCompare((PetscObject)pc_ctx->B_Ddelta,MATSHELL,&isshell);CHKERRQ(ierr);
997     if (isshell) {
998       BDdelta_DN ctx;
999       ierr = PetscViewerASCIIPrintf(viewer,"  FETI-DP BDdelta: DB^t * (B D^-1 B^t)^-1 for deluxe scaling (just from rank 0)\n");CHKERRQ(ierr);
1000       ierr = MatShellGetContext(pc_ctx->B_Ddelta,&ctx);CHKERRQ(ierr);
1001       if (!rank) {
1002         ierr = PetscViewerASCIIAddTab(sviewer,2);CHKERRQ(ierr);
1003         ierr = KSPView(ctx->kBD,sviewer);CHKERRQ(ierr);
1004         ierr = PetscViewerPushFormat(sviewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
1005         ierr = MatView(ctx->BD,sviewer);CHKERRQ(ierr);
1006         ierr = PetscViewerASCIISubtractTab(sviewer,2);CHKERRQ(ierr);
1007         ierr = PetscViewerPopFormat(sviewer);CHKERRQ(ierr);
1008       }
1009     }
1010     ierr = PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)pc_ctx->S_j),&sviewer);CHKERRQ(ierr);
1011     ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
1012     if (pc_ctx->kP) {
1013       ierr = PetscViewerASCIIPrintf(viewer,"  FETI-DP pressure preconditioner\n");CHKERRQ(ierr);
1014       ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr);
1015       ierr = KSPView(pc_ctx->kP,viewer);CHKERRQ(ierr);
1016       ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr);
1017     }
1018   }
1019   PetscFunctionReturn(0);
1020 }
1021