1 #include <../src/ksp/pc/impls/bddc/bddc.h>
2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3 #include <../src/mat/impls/dense/seq/dense.h>
4 #include <petscblaslapack.h>
5
6 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt*,PetscInt,PetscBT,PetscInt*,PetscInt*,PetscInt*);
7 static PetscErrorCode PCBDDCComputeExplicitSchur(Mat,PetscBool,MatReuse,Mat*);
8 static PetscErrorCode PCBDDCReuseSolvers_Interior(PC,Vec,Vec);
9 static PetscErrorCode PCBDDCReuseSolvers_Correction(PC,Vec,Vec);
10
11 /* if v2 is not present, correction is done in-place */
PCBDDCReuseSolversBenignAdapt(PCBDDCReuseSolvers ctx,Vec v,Vec v2,PetscBool sol,PetscBool full)12 PetscErrorCode PCBDDCReuseSolversBenignAdapt(PCBDDCReuseSolvers ctx, Vec v, Vec v2, PetscBool sol, PetscBool full)
13 {
14 PetscScalar *array;
15 PetscScalar *array2;
16 PetscErrorCode ierr;
17
18 PetscFunctionBegin;
19 if (!ctx->benign_n) PetscFunctionReturn(0);
20 if (sol && full) {
21 PetscInt n_I,size_schur;
22
23 /* get sizes */
24 ierr = MatGetSize(ctx->benign_csAIB,&size_schur,NULL);CHKERRQ(ierr);
25 ierr = VecGetSize(v,&n_I);CHKERRQ(ierr);
26 n_I = n_I - size_schur;
27 /* get schur sol from array */
28 ierr = VecGetArray(v,&array);CHKERRQ(ierr);
29 ierr = VecPlaceArray(ctx->benign_dummy_schur_vec,array+n_I);CHKERRQ(ierr);
30 ierr = VecRestoreArray(v,&array);CHKERRQ(ierr);
31 /* apply interior sol correction */
32 ierr = MatMultTranspose(ctx->benign_csAIB,ctx->benign_dummy_schur_vec,ctx->benign_corr_work);CHKERRQ(ierr);
33 ierr = VecResetArray(ctx->benign_dummy_schur_vec);CHKERRQ(ierr);
34 ierr = MatMultAdd(ctx->benign_AIIm1ones,ctx->benign_corr_work,v,v);CHKERRQ(ierr);
35 }
36 if (v2) {
37 PetscInt nl;
38
39 ierr = VecGetArrayRead(v,(const PetscScalar**)&array);CHKERRQ(ierr);
40 ierr = VecGetLocalSize(v2,&nl);CHKERRQ(ierr);
41 ierr = VecGetArray(v2,&array2);CHKERRQ(ierr);
42 ierr = PetscArraycpy(array2,array,nl);CHKERRQ(ierr);
43 } else {
44 ierr = VecGetArray(v,&array);CHKERRQ(ierr);
45 array2 = array;
46 }
47 if (!sol) { /* change rhs */
48 PetscInt n;
49 for (n=0;n<ctx->benign_n;n++) {
50 PetscScalar sum = 0.;
51 const PetscInt *cols;
52 PetscInt nz,i;
53
54 ierr = ISGetLocalSize(ctx->benign_zerodiag_subs[n],&nz);CHKERRQ(ierr);
55 ierr = ISGetIndices(ctx->benign_zerodiag_subs[n],&cols);CHKERRQ(ierr);
56 for (i=0;i<nz-1;i++) sum += array[cols[i]];
57 #if defined(PETSC_USE_COMPLEX)
58 sum = -(PetscRealPart(sum)/nz + PETSC_i*(PetscImaginaryPart(sum)/nz));
59 #else
60 sum = -sum/nz;
61 #endif
62 for (i=0;i<nz-1;i++) array2[cols[i]] += sum;
63 ctx->benign_save_vals[n] = array2[cols[nz-1]];
64 array2[cols[nz-1]] = sum;
65 ierr = ISRestoreIndices(ctx->benign_zerodiag_subs[n],&cols);CHKERRQ(ierr);
66 }
67 } else {
68 PetscInt n;
69 for (n=0;n<ctx->benign_n;n++) {
70 PetscScalar sum = 0.;
71 const PetscInt *cols;
72 PetscInt nz,i;
73 ierr = ISGetLocalSize(ctx->benign_zerodiag_subs[n],&nz);CHKERRQ(ierr);
74 ierr = ISGetIndices(ctx->benign_zerodiag_subs[n],&cols);CHKERRQ(ierr);
75 for (i=0;i<nz-1;i++) sum += array[cols[i]];
76 #if defined(PETSC_USE_COMPLEX)
77 sum = -(PetscRealPart(sum)/nz + PETSC_i*(PetscImaginaryPart(sum)/nz));
78 #else
79 sum = -sum/nz;
80 #endif
81 for (i=0;i<nz-1;i++) array2[cols[i]] += sum;
82 array2[cols[nz-1]] = ctx->benign_save_vals[n];
83 ierr = ISRestoreIndices(ctx->benign_zerodiag_subs[n],&cols);CHKERRQ(ierr);
84 }
85 }
86 if (v2) {
87 ierr = VecRestoreArrayRead(v,(const PetscScalar**)&array);CHKERRQ(ierr);
88 ierr = VecRestoreArray(v2,&array2);CHKERRQ(ierr);
89 } else {
90 ierr = VecRestoreArray(v,&array);CHKERRQ(ierr);
91 }
92 if (!sol && full) {
93 Vec usedv;
94 PetscInt n_I,size_schur;
95
96 /* get sizes */
97 ierr = MatGetSize(ctx->benign_csAIB,&size_schur,NULL);CHKERRQ(ierr);
98 ierr = VecGetSize(v,&n_I);CHKERRQ(ierr);
99 n_I = n_I - size_schur;
100 /* compute schur rhs correction */
101 if (v2) {
102 usedv = v2;
103 } else {
104 usedv = v;
105 }
106 /* apply schur rhs correction */
107 ierr = MatMultTranspose(ctx->benign_AIIm1ones,usedv,ctx->benign_corr_work);CHKERRQ(ierr);
108 ierr = VecGetArrayRead(usedv,(const PetscScalar**)&array);CHKERRQ(ierr);
109 ierr = VecPlaceArray(ctx->benign_dummy_schur_vec,array+n_I);CHKERRQ(ierr);
110 ierr = VecRestoreArrayRead(usedv,(const PetscScalar**)&array);CHKERRQ(ierr);
111 ierr = MatMultAdd(ctx->benign_csAIB,ctx->benign_corr_work,ctx->benign_dummy_schur_vec,ctx->benign_dummy_schur_vec);CHKERRQ(ierr);
112 ierr = VecResetArray(ctx->benign_dummy_schur_vec);CHKERRQ(ierr);
113 }
114 PetscFunctionReturn(0);
115 }
116
PCBDDCReuseSolvers_Solve_Private(PC pc,Vec rhs,Vec sol,PetscBool transpose,PetscBool full)117 static PetscErrorCode PCBDDCReuseSolvers_Solve_Private(PC pc, Vec rhs, Vec sol, PetscBool transpose, PetscBool full)
118 {
119 PCBDDCReuseSolvers ctx;
120 PetscBool copy = PETSC_FALSE;
121 PetscErrorCode ierr;
122
123 PetscFunctionBegin;
124 ierr = PCShellGetContext(pc,(void **)&ctx);CHKERRQ(ierr);
125 if (full) {
126 #if defined(PETSC_HAVE_MUMPS)
127 ierr = MatMumpsSetIcntl(ctx->F,26,-1);CHKERRQ(ierr);
128 #endif
129 #if defined(PETSC_HAVE_MKL_PARDISO)
130 ierr = MatMkl_PardisoSetCntl(ctx->F,70,0);CHKERRQ(ierr);
131 #endif
132 copy = ctx->has_vertices;
133 } else { /* interior solver */
134 #if defined(PETSC_HAVE_MUMPS)
135 ierr = MatMumpsSetIcntl(ctx->F,26,0);CHKERRQ(ierr);
136 #endif
137 #if defined(PETSC_HAVE_MKL_PARDISO)
138 ierr = MatMkl_PardisoSetCntl(ctx->F,70,1);CHKERRQ(ierr);
139 #endif
140 copy = PETSC_TRUE;
141 }
142 /* copy rhs into factored matrix workspace */
143 if (copy) {
144 PetscInt n;
145 PetscScalar *array,*array_solver;
146
147 ierr = VecGetLocalSize(rhs,&n);CHKERRQ(ierr);
148 ierr = VecGetArrayRead(rhs,(const PetscScalar**)&array);CHKERRQ(ierr);
149 ierr = VecGetArray(ctx->rhs,&array_solver);CHKERRQ(ierr);
150 ierr = PetscArraycpy(array_solver,array,n);CHKERRQ(ierr);
151 ierr = VecRestoreArray(ctx->rhs,&array_solver);CHKERRQ(ierr);
152 ierr = VecRestoreArrayRead(rhs,(const PetscScalar**)&array);CHKERRQ(ierr);
153
154 ierr = PCBDDCReuseSolversBenignAdapt(ctx,ctx->rhs,NULL,PETSC_FALSE,full);CHKERRQ(ierr);
155 if (transpose) {
156 ierr = MatSolveTranspose(ctx->F,ctx->rhs,ctx->sol);CHKERRQ(ierr);
157 } else {
158 ierr = MatSolve(ctx->F,ctx->rhs,ctx->sol);CHKERRQ(ierr);
159 }
160 ierr = PCBDDCReuseSolversBenignAdapt(ctx,ctx->sol,NULL,PETSC_TRUE,full);CHKERRQ(ierr);
161
162 /* get back data to caller worskpace */
163 ierr = VecGetArrayRead(ctx->sol,(const PetscScalar**)&array_solver);CHKERRQ(ierr);
164 ierr = VecGetArray(sol,&array);CHKERRQ(ierr);
165 ierr = PetscArraycpy(array,array_solver,n);CHKERRQ(ierr);
166 ierr = VecRestoreArray(sol,&array);CHKERRQ(ierr);
167 ierr = VecRestoreArrayRead(ctx->sol,(const PetscScalar**)&array_solver);CHKERRQ(ierr);
168 } else {
169 if (ctx->benign_n) {
170 ierr = PCBDDCReuseSolversBenignAdapt(ctx,rhs,ctx->rhs,PETSC_FALSE,full);CHKERRQ(ierr);
171 if (transpose) {
172 ierr = MatSolveTranspose(ctx->F,ctx->rhs,sol);CHKERRQ(ierr);
173 } else {
174 ierr = MatSolve(ctx->F,ctx->rhs,sol);CHKERRQ(ierr);
175 }
176 ierr = PCBDDCReuseSolversBenignAdapt(ctx,sol,NULL,PETSC_TRUE,full);CHKERRQ(ierr);
177 } else {
178 if (transpose) {
179 ierr = MatSolveTranspose(ctx->F,rhs,sol);CHKERRQ(ierr);
180 } else {
181 ierr = MatSolve(ctx->F,rhs,sol);CHKERRQ(ierr);
182 }
183 }
184 }
185 /* restore defaults */
186 #if defined(PETSC_HAVE_MUMPS)
187 ierr = MatMumpsSetIcntl(ctx->F,26,-1);CHKERRQ(ierr);
188 #endif
189 #if defined(PETSC_HAVE_MKL_PARDISO)
190 ierr = MatMkl_PardisoSetCntl(ctx->F,70,0);CHKERRQ(ierr);
191 #endif
192 PetscFunctionReturn(0);
193 }
194
PCBDDCReuseSolvers_Correction(PC pc,Vec rhs,Vec sol)195 static PetscErrorCode PCBDDCReuseSolvers_Correction(PC pc, Vec rhs, Vec sol)
196 {
197 PetscErrorCode ierr;
198
199 PetscFunctionBegin;
200 ierr = PCBDDCReuseSolvers_Solve_Private(pc,rhs,sol,PETSC_FALSE,PETSC_TRUE);CHKERRQ(ierr);
201 PetscFunctionReturn(0);
202 }
203
PCBDDCReuseSolvers_CorrectionTranspose(PC pc,Vec rhs,Vec sol)204 static PetscErrorCode PCBDDCReuseSolvers_CorrectionTranspose(PC pc, Vec rhs, Vec sol)
205 {
206 PetscErrorCode ierr;
207
208 PetscFunctionBegin;
209 ierr = PCBDDCReuseSolvers_Solve_Private(pc,rhs,sol,PETSC_TRUE,PETSC_TRUE);CHKERRQ(ierr);
210 PetscFunctionReturn(0);
211 }
212
PCBDDCReuseSolvers_Interior(PC pc,Vec rhs,Vec sol)213 static PetscErrorCode PCBDDCReuseSolvers_Interior(PC pc, Vec rhs, Vec sol)
214 {
215 PetscErrorCode ierr;
216
217 PetscFunctionBegin;
218 ierr = PCBDDCReuseSolvers_Solve_Private(pc,rhs,sol,PETSC_FALSE,PETSC_FALSE);CHKERRQ(ierr);
219 PetscFunctionReturn(0);
220 }
221
PCBDDCReuseSolvers_InteriorTranspose(PC pc,Vec rhs,Vec sol)222 static PetscErrorCode PCBDDCReuseSolvers_InteriorTranspose(PC pc, Vec rhs, Vec sol)
223 {
224 PetscErrorCode ierr;
225
226 PetscFunctionBegin;
227 ierr = PCBDDCReuseSolvers_Solve_Private(pc,rhs,sol,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr);
228 PetscFunctionReturn(0);
229 }
230
PCBDDCReuseSolvers_View(PC pc,PetscViewer viewer)231 static PetscErrorCode PCBDDCReuseSolvers_View(PC pc, PetscViewer viewer)
232 {
233 PCBDDCReuseSolvers ctx;
234 PetscBool iascii;
235 PetscErrorCode ierr;
236
237 PetscFunctionBegin;
238 ierr = PCShellGetContext(pc,(void **)&ctx);CHKERRQ(ierr);
239 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
240 if (iascii) {
241 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
242 }
243 ierr = MatView(ctx->F,viewer);CHKERRQ(ierr);
244 if (iascii) {
245 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
246 }
247 PetscFunctionReturn(0);
248 }
249
PCBDDCReuseSolversReset(PCBDDCReuseSolvers reuse)250 static PetscErrorCode PCBDDCReuseSolversReset(PCBDDCReuseSolvers reuse)
251 {
252 PetscInt i;
253 PetscErrorCode ierr;
254
255 PetscFunctionBegin;
256 ierr = MatDestroy(&reuse->F);CHKERRQ(ierr);
257 ierr = VecDestroy(&reuse->sol);CHKERRQ(ierr);
258 ierr = VecDestroy(&reuse->rhs);CHKERRQ(ierr);
259 ierr = PCDestroy(&reuse->interior_solver);CHKERRQ(ierr);
260 ierr = PCDestroy(&reuse->correction_solver);CHKERRQ(ierr);
261 ierr = ISDestroy(&reuse->is_R);CHKERRQ(ierr);
262 ierr = ISDestroy(&reuse->is_B);CHKERRQ(ierr);
263 ierr = VecScatterDestroy(&reuse->correction_scatter_B);CHKERRQ(ierr);
264 ierr = VecDestroy(&reuse->sol_B);CHKERRQ(ierr);
265 ierr = VecDestroy(&reuse->rhs_B);CHKERRQ(ierr);
266 for (i=0;i<reuse->benign_n;i++) {
267 ierr = ISDestroy(&reuse->benign_zerodiag_subs[i]);CHKERRQ(ierr);
268 }
269 ierr = PetscFree(reuse->benign_zerodiag_subs);CHKERRQ(ierr);
270 ierr = PetscFree(reuse->benign_save_vals);CHKERRQ(ierr);
271 ierr = MatDestroy(&reuse->benign_csAIB);CHKERRQ(ierr);
272 ierr = MatDestroy(&reuse->benign_AIIm1ones);CHKERRQ(ierr);
273 ierr = VecDestroy(&reuse->benign_corr_work);CHKERRQ(ierr);
274 ierr = VecDestroy(&reuse->benign_dummy_schur_vec);CHKERRQ(ierr);
275 PetscFunctionReturn(0);
276 }
277
PCBDDCComputeExplicitSchur(Mat M,PetscBool issym,MatReuse reuse,Mat * S)278 static PetscErrorCode PCBDDCComputeExplicitSchur(Mat M, PetscBool issym, MatReuse reuse, Mat *S)
279 {
280 Mat B, C, D, Bd, Cd, AinvBd;
281 KSP ksp;
282 PC pc;
283 PetscBool isLU, isILU, isCHOL, Bdense, Cdense;
284 PetscReal fill = 2.0;
285 PetscInt n_I;
286 PetscMPIInt size;
287 PetscErrorCode ierr;
288
289 PetscFunctionBegin;
290 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)M),&size);CHKERRQ(ierr);
291 if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parallel matrices");
292 if (reuse == MAT_REUSE_MATRIX) {
293 PetscBool Sdense;
294
295 ierr = PetscObjectTypeCompare((PetscObject)*S, MATSEQDENSE, &Sdense);CHKERRQ(ierr);
296 if (!Sdense) SETERRQ(PetscObjectComm((PetscObject)M),PETSC_ERR_SUP,"S should dense");
297 }
298 ierr = MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);CHKERRQ(ierr);
299 ierr = MatSchurComplementGetKSP(M, &ksp);CHKERRQ(ierr);
300 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
301 ierr = PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);CHKERRQ(ierr);
302 ierr = PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);CHKERRQ(ierr);
303 ierr = PetscObjectTypeCompare((PetscObject) pc, PCCHOLESKY, &isCHOL);CHKERRQ(ierr);
304 ierr = PetscObjectTypeCompare((PetscObject) B, MATSEQDENSE, &Bdense);CHKERRQ(ierr);
305 ierr = PetscObjectTypeCompare((PetscObject) C, MATSEQDENSE, &Cdense);CHKERRQ(ierr);
306 ierr = MatGetSize(B,&n_I,NULL);CHKERRQ(ierr);
307 if (n_I) {
308 if (!Bdense) {
309 ierr = MatConvert(B, MATSEQDENSE, MAT_INITIAL_MATRIX, &Bd);CHKERRQ(ierr);
310 } else {
311 Bd = B;
312 }
313
314 if (isLU || isILU || isCHOL) {
315 Mat fact;
316 ierr = KSPSetUp(ksp);CHKERRQ(ierr);
317 ierr = PCFactorGetMatrix(pc, &fact);CHKERRQ(ierr);
318 ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
319 ierr = MatMatSolve(fact, Bd, AinvBd);CHKERRQ(ierr);
320 } else {
321 PetscBool ex = PETSC_TRUE;
322
323 if (ex) {
324 Mat Ainvd;
325
326 ierr = PCComputeOperator(pc, MATDENSE, &Ainvd);CHKERRQ(ierr);
327 ierr = MatMatMult(Ainvd, Bd, MAT_INITIAL_MATRIX, fill, &AinvBd);CHKERRQ(ierr);
328 ierr = MatDestroy(&Ainvd);CHKERRQ(ierr);
329 } else {
330 Vec sol,rhs;
331 PetscScalar *arrayrhs,*arraysol;
332 PetscInt i,nrhs,n;
333
334 ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
335 ierr = MatGetSize(Bd,&n,&nrhs);CHKERRQ(ierr);
336 ierr = MatDenseGetArray(Bd,&arrayrhs);CHKERRQ(ierr);
337 ierr = MatDenseGetArray(AinvBd,&arraysol);CHKERRQ(ierr);
338 ierr = KSPGetSolution(ksp,&sol);CHKERRQ(ierr);
339 ierr = KSPGetRhs(ksp,&rhs);CHKERRQ(ierr);
340 for (i=0;i<nrhs;i++) {
341 ierr = VecPlaceArray(rhs,arrayrhs+i*n);CHKERRQ(ierr);
342 ierr = VecPlaceArray(sol,arraysol+i*n);CHKERRQ(ierr);
343 ierr = KSPSolve(ksp,rhs,sol);CHKERRQ(ierr);
344 ierr = VecResetArray(rhs);CHKERRQ(ierr);
345 ierr = VecResetArray(sol);CHKERRQ(ierr);
346 }
347 ierr = MatDenseRestoreArray(Bd,&arrayrhs);CHKERRQ(ierr);
348 ierr = MatDenseRestoreArray(AinvBd,&arrayrhs);CHKERRQ(ierr);
349 }
350 }
351 if (!Bdense & !issym) {
352 ierr = MatDestroy(&Bd);CHKERRQ(ierr);
353 }
354
355 if (!issym) {
356 if (!Cdense) {
357 ierr = MatConvert(C, MATSEQDENSE, MAT_INITIAL_MATRIX, &Cd);CHKERRQ(ierr);
358 } else {
359 Cd = C;
360 }
361 ierr = MatMatMult(Cd, AinvBd, reuse, fill, S);CHKERRQ(ierr);
362 if (!Cdense) {
363 ierr = MatDestroy(&Cd);CHKERRQ(ierr);
364 }
365 } else {
366 ierr = MatTransposeMatMult(Bd, AinvBd, reuse, fill, S);CHKERRQ(ierr);
367 if (!Bdense) {
368 ierr = MatDestroy(&Bd);CHKERRQ(ierr);
369 }
370 }
371 ierr = MatDestroy(&AinvBd);CHKERRQ(ierr);
372 }
373
374 if (D) {
375 Mat Dd;
376 PetscBool Ddense;
377
378 ierr = PetscObjectTypeCompare((PetscObject)D,MATSEQDENSE,&Ddense);CHKERRQ(ierr);
379 if (!Ddense) {
380 ierr = MatConvert(D, MATSEQDENSE, MAT_INITIAL_MATRIX, &Dd);CHKERRQ(ierr);
381 } else {
382 Dd = D;
383 }
384 if (n_I) {
385 ierr = MatAYPX(*S,-1.0,Dd,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
386 } else {
387 if (reuse == MAT_INITIAL_MATRIX) {
388 ierr = MatDuplicate(Dd,MAT_COPY_VALUES,S);CHKERRQ(ierr);
389 } else {
390 ierr = MatCopy(Dd,*S,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
391 }
392 }
393 if (!Ddense) {
394 ierr = MatDestroy(&Dd);CHKERRQ(ierr);
395 }
396 } else {
397 ierr = MatScale(*S,-1.0);CHKERRQ(ierr);
398 }
399 PetscFunctionReturn(0);
400 }
401
PCBDDCSubSchursSetUp(PCBDDCSubSchurs sub_schurs,Mat Ain,Mat Sin,PetscBool exact_schur,PetscInt xadj[],PetscInt adjncy[],PetscInt nlayers,Vec scaling,PetscBool compute_Stilda,PetscBool reuse_solvers,PetscBool benign_trick,PetscInt benign_n,PetscInt benign_p0_lidx[],IS benign_zerodiag_subs[],Mat change,IS change_primal)402 PetscErrorCode PCBDDCSubSchursSetUp(PCBDDCSubSchurs sub_schurs, Mat Ain, Mat Sin, PetscBool exact_schur, PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers, Vec scaling, PetscBool compute_Stilda, PetscBool reuse_solvers, PetscBool benign_trick, PetscInt benign_n, PetscInt benign_p0_lidx[], IS benign_zerodiag_subs[], Mat change, IS change_primal)
403 {
404 Mat F,A_II,A_IB,A_BI,A_BB,AE_II;
405 Mat S_all;
406 Vec gstash,lstash;
407 VecScatter sstash;
408 IS is_I,is_I_layer;
409 IS all_subsets,all_subsets_mult,all_subsets_n;
410 PetscScalar *stasharray,*Bwork;
411 PetscInt *nnz,*all_local_idx_N;
412 PetscInt *auxnum1,*auxnum2;
413 PetscInt i,subset_size,max_subset_size;
414 PetscInt n_B,extra,local_size,global_size;
415 PetscInt local_stash_size;
416 PetscBLASInt B_N,B_ierr,B_lwork,*pivots;
417 MPI_Comm comm_n;
418 PetscBool deluxe = PETSC_TRUE;
419 PetscBool use_potr = PETSC_FALSE, use_sytr = PETSC_FALSE;
420 PetscViewer matl_dbg_viewer = NULL;
421 PetscErrorCode ierr;
422
423 PetscFunctionBegin;
424 ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr);
425 ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr);
426 if (Ain) {
427 ierr = PetscObjectReference((PetscObject)Ain);CHKERRQ(ierr);
428 sub_schurs->A = Ain;
429 }
430
431 ierr = PetscObjectReference((PetscObject)Sin);CHKERRQ(ierr);
432 sub_schurs->S = Sin;
433 if (sub_schurs->schur_explicit) {
434 sub_schurs->schur_explicit = (PetscBool)(!!sub_schurs->A);
435 }
436
437 /* preliminary checks */
438 if (!sub_schurs->schur_explicit && compute_Stilda) SETERRQ(PetscObjectComm((PetscObject)sub_schurs->l2gmap),PETSC_ERR_SUP,"Adaptive selection of constraints requires MUMPS and/or MKL_PARDISO");
439
440 if (benign_trick) sub_schurs->is_posdef = PETSC_FALSE;
441
442 /* debug (MATLAB) */
443 if (sub_schurs->debug) {
444 PetscMPIInt size,rank;
445 PetscInt nr,*print_schurs_ranks,print_schurs = PETSC_FALSE;
446 PetscBool flg;
447
448 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)sub_schurs->l2gmap),&size);CHKERRQ(ierr);
449 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)sub_schurs->l2gmap),&rank);CHKERRQ(ierr);
450 nr = size;
451 ierr = PetscMalloc1(nr,&print_schurs_ranks);CHKERRQ(ierr);
452 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)sub_schurs->l2gmap),sub_schurs->prefix,"BDDC sub_schurs options","PC");CHKERRQ(ierr);
453 ierr = PetscOptionsIntArray("-sub_schurs_debug_ranks","Ranks to debug (all if the option is not used)",NULL,print_schurs_ranks,&nr,&flg);CHKERRQ(ierr);
454 if (!flg) print_schurs = PETSC_TRUE;
455 else {
456 print_schurs = PETSC_FALSE;
457 for (i=0;i<nr;i++) if (print_schurs_ranks[i] == (PetscInt)rank) { print_schurs = PETSC_TRUE; break; }
458 }
459 ierr = PetscOptionsEnd();CHKERRQ(ierr);
460 ierr = PetscFree(print_schurs_ranks);CHKERRQ(ierr);
461 if (print_schurs) {
462 char filename[256];
463
464 ierr = PetscSNPrintf(filename,sizeof(filename),"sub_schurs_Schur_r%d.m",PetscGlobalRank);CHKERRQ(ierr);
465 ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF,filename,&matl_dbg_viewer);CHKERRQ(ierr);
466 ierr = PetscViewerPushFormat(matl_dbg_viewer,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
467 }
468 }
469
470
471 /* restrict work on active processes */
472 if (sub_schurs->restrict_comm) {
473 PetscSubcomm subcomm;
474 PetscMPIInt color,rank;
475
476 color = 0;
477 if (!sub_schurs->n_subs) color = 1; /* this can happen if we are in a multilevel case or if the subdomain is disconnected */
478 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)sub_schurs->l2gmap),&rank);CHKERRQ(ierr);
479 ierr = PetscSubcommCreate(PetscObjectComm((PetscObject)sub_schurs->l2gmap),&subcomm);CHKERRQ(ierr);
480 ierr = PetscSubcommSetNumber(subcomm,2);CHKERRQ(ierr);
481 ierr = PetscSubcommSetTypeGeneral(subcomm,color,rank);CHKERRQ(ierr);
482 ierr = PetscCommDuplicate(PetscSubcommChild(subcomm),&comm_n,NULL);CHKERRQ(ierr);
483 ierr = PetscSubcommDestroy(&subcomm);CHKERRQ(ierr);
484 if (!sub_schurs->n_subs) {
485 ierr = PetscCommDestroy(&comm_n);CHKERRQ(ierr);
486 PetscFunctionReturn(0);
487 }
488 } else {
489 ierr = PetscCommDuplicate(PetscObjectComm((PetscObject)sub_schurs->l2gmap),&comm_n,NULL);CHKERRQ(ierr);
490 }
491
492 /* get Schur complement matrices */
493 if (!sub_schurs->schur_explicit) {
494 Mat tA_IB,tA_BI,tA_BB;
495 PetscBool isseqsbaij;
496 ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,&tA_IB,&tA_BI,&tA_BB);CHKERRQ(ierr);
497 ierr = PetscObjectTypeCompare((PetscObject)tA_BB,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr);
498 if (isseqsbaij) {
499 ierr = MatConvert(tA_BB,MATSEQAIJ,MAT_INITIAL_MATRIX,&A_BB);CHKERRQ(ierr);
500 ierr = MatConvert(tA_IB,MATSEQAIJ,MAT_INITIAL_MATRIX,&A_IB);CHKERRQ(ierr);
501 ierr = MatConvert(tA_BI,MATSEQAIJ,MAT_INITIAL_MATRIX,&A_BI);CHKERRQ(ierr);
502 } else {
503 ierr = PetscObjectReference((PetscObject)tA_BB);CHKERRQ(ierr);
504 A_BB = tA_BB;
505 ierr = PetscObjectReference((PetscObject)tA_IB);CHKERRQ(ierr);
506 A_IB = tA_IB;
507 ierr = PetscObjectReference((PetscObject)tA_BI);CHKERRQ(ierr);
508 A_BI = tA_BI;
509 }
510 } else {
511 A_II = NULL;
512 A_IB = NULL;
513 A_BI = NULL;
514 A_BB = NULL;
515 }
516 S_all = NULL;
517
518 /* determine interior problems */
519 ierr = ISGetLocalSize(sub_schurs->is_I,&i);CHKERRQ(ierr);
520 if (nlayers >= 0 && i) { /* Interior problems can be different from the original one */
521 PetscBT touched;
522 const PetscInt* idx_B;
523 PetscInt n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering;
524
525 if (!xadj) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot request layering without adjacency");
526 /* get sizes */
527 ierr = ISGetLocalSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
528 ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr);
529
530 ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr);
531 ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr);
532 ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr);
533
534 /* all boundary dofs must be skipped when adding layers */
535 ierr = ISGetIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
536 for (j=0;j<n_B;j++) {
537 ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr);
538 }
539 ierr = PetscArraycpy(local_numbering,idx_B,n_B);CHKERRQ(ierr);
540 ierr = ISRestoreIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr);
541
542 /* add prescribed number of layers of dofs */
543 n_local_dofs = n_B;
544 n_prev_added = n_B;
545 for (layer=0;layer<nlayers;layer++) {
546 PetscInt n_added;
547 if (n_local_dofs == n_I+n_B) break;
548 if (n_local_dofs > n_I+n_B) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error querying layer %D. Out of bound access (%D > %D)",layer,n_local_dofs,n_I+n_B);
549 ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr);
550 n_prev_added = n_added;
551 n_local_dofs += n_added;
552 if (!n_added) break;
553 }
554 ierr = PetscBTDestroy(&touched);CHKERRQ(ierr);
555
556 /* IS for I layer dofs in original numbering */
557 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->is_I),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&is_I_layer);CHKERRQ(ierr);
558 ierr = PetscFree(local_numbering);CHKERRQ(ierr);
559 ierr = ISSort(is_I_layer);CHKERRQ(ierr);
560 /* IS for I layer dofs in I numbering */
561 if (!sub_schurs->schur_explicit) {
562 ISLocalToGlobalMapping ItoNmap;
563 ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_I,&ItoNmap);CHKERRQ(ierr);
564 ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,is_I_layer,&is_I);CHKERRQ(ierr);
565 ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr);
566
567 /* II block */
568 ierr = MatCreateSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr);
569 }
570 } else {
571 PetscInt n_I;
572
573 /* IS for I dofs in original numbering */
574 ierr = PetscObjectReference((PetscObject)sub_schurs->is_I);CHKERRQ(ierr);
575 is_I_layer = sub_schurs->is_I;
576
577 /* IS for I dofs in I numbering (strided 1) */
578 if (!sub_schurs->schur_explicit) {
579 ierr = ISGetSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr);
580 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),n_I,0,1,&is_I);CHKERRQ(ierr);
581
582 /* II block is the same */
583 ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr);
584 AE_II = A_II;
585 }
586 }
587
588 /* Get info on subset sizes and sum of all subsets sizes */
589 max_subset_size = 0;
590 local_size = 0;
591 for (i=0;i<sub_schurs->n_subs;i++) {
592 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
593 max_subset_size = PetscMax(subset_size,max_subset_size);
594 local_size += subset_size;
595 }
596
597 /* Work arrays for local indices */
598 extra = 0;
599 ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr);
600 if (sub_schurs->schur_explicit && is_I_layer) {
601 ierr = ISGetLocalSize(is_I_layer,&extra);CHKERRQ(ierr);
602 }
603 ierr = PetscMalloc1(n_B+extra,&all_local_idx_N);CHKERRQ(ierr);
604 if (extra) {
605 const PetscInt *idxs;
606 ierr = ISGetIndices(is_I_layer,&idxs);CHKERRQ(ierr);
607 ierr = PetscArraycpy(all_local_idx_N,idxs,extra);CHKERRQ(ierr);
608 ierr = ISRestoreIndices(is_I_layer,&idxs);CHKERRQ(ierr);
609 }
610 ierr = PetscMalloc1(sub_schurs->n_subs,&auxnum1);CHKERRQ(ierr);
611 ierr = PetscMalloc1(sub_schurs->n_subs,&auxnum2);CHKERRQ(ierr);
612
613 /* Get local indices in local numbering */
614 local_size = 0;
615 local_stash_size = 0;
616 for (i=0;i<sub_schurs->n_subs;i++) {
617 const PetscInt *idxs;
618
619 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
620 ierr = ISGetIndices(sub_schurs->is_subs[i],&idxs);CHKERRQ(ierr);
621 /* start (smallest in global ordering) and multiplicity */
622 auxnum1[i] = idxs[0];
623 auxnum2[i] = subset_size*subset_size;
624 /* subset indices in local numbering */
625 ierr = PetscArraycpy(all_local_idx_N+local_size+extra,idxs,subset_size);CHKERRQ(ierr);
626 ierr = ISRestoreIndices(sub_schurs->is_subs[i],&idxs);CHKERRQ(ierr);
627 local_size += subset_size;
628 local_stash_size += subset_size*subset_size;
629 }
630
631 /* allocate extra workspace needed only for GETRI or SYTRF */
632 use_potr = use_sytr = PETSC_FALSE;
633 if (benign_trick || (sub_schurs->is_hermitian && sub_schurs->is_posdef)) {
634 use_potr = PETSC_TRUE;
635 } else if (sub_schurs->is_symmetric) {
636 use_sytr = PETSC_TRUE;
637 }
638 if (local_size && !use_potr) {
639 PetscScalar lwork,dummyscalar = 0.;
640 PetscBLASInt dummyint = 0;
641
642 B_lwork = -1;
643 ierr = PetscBLASIntCast(local_size,&B_N);CHKERRQ(ierr);
644 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
645 if (use_sytr) {
646 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&B_N,&dummyscalar,&B_N,&dummyint,&lwork,&B_lwork,&B_ierr));
647 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SYTRF Lapack routine %d",(int)B_ierr);
648 } else {
649 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,&dummyscalar,&B_N,&dummyint,&lwork,&B_lwork,&B_ierr));
650 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GETRI Lapack routine %d",(int)B_ierr);
651 }
652 ierr = PetscFPTrapPop();CHKERRQ(ierr);
653 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(lwork),&B_lwork);CHKERRQ(ierr);
654 ierr = PetscMalloc2(B_lwork,&Bwork,B_N,&pivots);CHKERRQ(ierr);
655 } else {
656 Bwork = NULL;
657 pivots = NULL;
658 }
659
660 /* prepare data for summing up properly schurs on subsets */
661 ierr = ISCreateGeneral(comm_n,sub_schurs->n_subs,auxnum1,PETSC_OWN_POINTER,&all_subsets_n);CHKERRQ(ierr);
662 ierr = ISLocalToGlobalMappingApplyIS(sub_schurs->l2gmap,all_subsets_n,&all_subsets);CHKERRQ(ierr);
663 ierr = ISDestroy(&all_subsets_n);CHKERRQ(ierr);
664 ierr = ISCreateGeneral(comm_n,sub_schurs->n_subs,auxnum2,PETSC_OWN_POINTER,&all_subsets_mult);CHKERRQ(ierr);
665 ierr = ISRenumber(all_subsets,all_subsets_mult,&global_size,&all_subsets_n);CHKERRQ(ierr);
666 ierr = ISDestroy(&all_subsets);CHKERRQ(ierr);
667 ierr = ISDestroy(&all_subsets_mult);CHKERRQ(ierr);
668 ierr = ISGetLocalSize(all_subsets_n,&i);CHKERRQ(ierr);
669 if (i != local_stash_size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid size of new subset! %D != %D",i,local_stash_size);
670 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,local_stash_size,NULL,&lstash);CHKERRQ(ierr);
671 ierr = VecCreateMPI(comm_n,PETSC_DECIDE,global_size,&gstash);CHKERRQ(ierr);
672 ierr = VecScatterCreate(lstash,NULL,gstash,all_subsets_n,&sstash);CHKERRQ(ierr);
673 ierr = ISDestroy(&all_subsets_n);CHKERRQ(ierr);
674
675 /* subset indices in local boundary numbering */
676 if (!sub_schurs->is_Ej_all) {
677 PetscInt *all_local_idx_B;
678
679 ierr = PetscMalloc1(local_size,&all_local_idx_B);CHKERRQ(ierr);
680 ierr = ISGlobalToLocalMappingApply(sub_schurs->BtoNmap,IS_GTOLM_DROP,local_size,all_local_idx_N+extra,&subset_size,all_local_idx_B);CHKERRQ(ierr);
681 if (subset_size != local_size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in sub_schurs serial (BtoNmap)! %D != %D",subset_size,local_size);
682 ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_B,PETSC_OWN_POINTER,&sub_schurs->is_Ej_all);CHKERRQ(ierr);
683 }
684
685 if (change) {
686 ISLocalToGlobalMapping BtoS;
687 IS change_primal_B;
688 IS change_primal_all;
689
690 if (sub_schurs->change_primal_sub) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
691 if (sub_schurs->change) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
692 ierr = PetscMalloc1(sub_schurs->n_subs,&sub_schurs->change_primal_sub);CHKERRQ(ierr);
693 for (i=0;i<sub_schurs->n_subs;i++) {
694 ISLocalToGlobalMapping NtoS;
695 ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_subs[i],&NtoS);CHKERRQ(ierr);
696 ierr = ISGlobalToLocalMappingApplyIS(NtoS,IS_GTOLM_DROP,change_primal,&sub_schurs->change_primal_sub[i]);CHKERRQ(ierr);
697 ierr = ISLocalToGlobalMappingDestroy(&NtoS);CHKERRQ(ierr);
698 }
699 ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,change_primal,&change_primal_B);CHKERRQ(ierr);
700 ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_Ej_all,&BtoS);CHKERRQ(ierr);
701 ierr = ISGlobalToLocalMappingApplyIS(BtoS,IS_GTOLM_DROP,change_primal_B,&change_primal_all);CHKERRQ(ierr);
702 ierr = ISLocalToGlobalMappingDestroy(&BtoS);CHKERRQ(ierr);
703 ierr = ISDestroy(&change_primal_B);CHKERRQ(ierr);
704 ierr = PetscMalloc1(sub_schurs->n_subs,&sub_schurs->change);CHKERRQ(ierr);
705 for (i=0;i<sub_schurs->n_subs;i++) {
706 Mat change_sub;
707
708 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
709 ierr = KSPCreate(PETSC_COMM_SELF,&sub_schurs->change[i]);CHKERRQ(ierr);
710 ierr = KSPSetType(sub_schurs->change[i],KSPPREONLY);CHKERRQ(ierr);
711 if (!sub_schurs->change_with_qr) {
712 ierr = MatCreateSubMatrix(change,sub_schurs->is_subs[i],sub_schurs->is_subs[i],MAT_INITIAL_MATRIX,&change_sub);CHKERRQ(ierr);
713 } else {
714 Mat change_subt;
715 ierr = MatCreateSubMatrix(change,sub_schurs->is_subs[i],sub_schurs->is_subs[i],MAT_INITIAL_MATRIX,&change_subt);CHKERRQ(ierr);
716 ierr = MatConvert(change_subt,MATSEQDENSE,MAT_INITIAL_MATRIX,&change_sub);CHKERRQ(ierr);
717 ierr = MatDestroy(&change_subt);CHKERRQ(ierr);
718 }
719 ierr = KSPSetOperators(sub_schurs->change[i],change_sub,change_sub);CHKERRQ(ierr);
720 ierr = MatDestroy(&change_sub);CHKERRQ(ierr);
721 ierr = KSPSetOptionsPrefix(sub_schurs->change[i],sub_schurs->prefix);CHKERRQ(ierr);
722 ierr = KSPAppendOptionsPrefix(sub_schurs->change[i],"sub_schurs_change_");CHKERRQ(ierr);
723 }
724 ierr = ISDestroy(&change_primal_all);CHKERRQ(ierr);
725 }
726
727 /* Local matrix of all local Schur on subsets (transposed) */
728 if (!sub_schurs->S_Ej_all) {
729 Mat T;
730 PetscScalar *v;
731 PetscInt *ii,*jj;
732 PetscInt cum,i,j,k;
733
734 /* MatSeqAIJSetPreallocation + MatSetValues is slow for these kind of matrices (may have large blocks)
735 Allocate properly a representative matrix and duplicate */
736 ierr = PetscMalloc3(local_size+1,&ii,local_stash_size,&jj,local_stash_size,&v);CHKERRQ(ierr);
737 ii[0] = 0;
738 cum = 0;
739 for (i=0;i<sub_schurs->n_subs;i++) {
740 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
741 for (j=0;j<subset_size;j++) {
742 const PetscInt row = cum+j;
743 PetscInt col = cum;
744
745 ii[row+1] = ii[row] + subset_size;
746 for (k=ii[row];k<ii[row+1];k++) {
747 jj[k] = col;
748 col++;
749 }
750 }
751 cum += subset_size;
752 }
753 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,local_size,local_size,ii,jj,v,&T);CHKERRQ(ierr);
754 ierr = MatDuplicate(T,MAT_DO_NOT_COPY_VALUES,&sub_schurs->S_Ej_all);CHKERRQ(ierr);
755 ierr = MatDestroy(&T);CHKERRQ(ierr);
756 ierr = PetscFree3(ii,jj,v);CHKERRQ(ierr);
757 }
758 /* matrices for deluxe scaling and adaptive selection */
759 if (compute_Stilda) {
760 if (!sub_schurs->sum_S_Ej_tilda_all) {
761 ierr = MatDuplicate(sub_schurs->S_Ej_all,MAT_DO_NOT_COPY_VALUES,&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
762 }
763 if (!sub_schurs->sum_S_Ej_inv_all && deluxe) {
764 ierr = MatDuplicate(sub_schurs->S_Ej_all,MAT_DO_NOT_COPY_VALUES,&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
765 }
766 }
767
768 /* Compute Schur complements explicitly */
769 F = NULL;
770 if (!sub_schurs->schur_explicit) {
771 /* this code branch is used when MatFactor with Schur complement support is not present or when explicitly requested;
772 it is not efficient, unless the economic version of the scaling is used */
773 Mat S_Ej_expl;
774 PetscScalar *work;
775 PetscInt j,*dummy_idx;
776 PetscBool Sdense;
777
778 ierr = PetscMalloc2(max_subset_size,&dummy_idx,max_subset_size*max_subset_size,&work);CHKERRQ(ierr);
779 local_size = 0;
780 for (i=0;i<sub_schurs->n_subs;i++) {
781 IS is_subset_B;
782 Mat AE_EE,AE_IE,AE_EI,S_Ej;
783
784 /* subsets in original and boundary numbering */
785 ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,sub_schurs->is_subs[i],&is_subset_B);CHKERRQ(ierr);
786 /* EE block */
787 ierr = MatCreateSubMatrix(A_BB,is_subset_B,is_subset_B,MAT_INITIAL_MATRIX,&AE_EE);CHKERRQ(ierr);
788 /* IE block */
789 ierr = MatCreateSubMatrix(A_IB,is_I,is_subset_B,MAT_INITIAL_MATRIX,&AE_IE);CHKERRQ(ierr);
790 /* EI block */
791 if (sub_schurs->is_symmetric) {
792 ierr = MatCreateTranspose(AE_IE,&AE_EI);CHKERRQ(ierr);
793 } else if (sub_schurs->is_hermitian) {
794 ierr = MatCreateHermitianTranspose(AE_IE,&AE_EI);CHKERRQ(ierr);
795 } else {
796 ierr = MatCreateSubMatrix(A_BI,is_subset_B,is_I,MAT_INITIAL_MATRIX,&AE_EI);CHKERRQ(ierr);
797 }
798 ierr = ISDestroy(&is_subset_B);CHKERRQ(ierr);
799 ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE,AE_EI,AE_EE,&S_Ej);CHKERRQ(ierr);
800 ierr = MatDestroy(&AE_EE);CHKERRQ(ierr);
801 ierr = MatDestroy(&AE_IE);CHKERRQ(ierr);
802 ierr = MatDestroy(&AE_EI);CHKERRQ(ierr);
803 if (AE_II == A_II) { /* we can reuse the same ksp */
804 KSP ksp;
805 ierr = MatSchurComplementGetKSP(sub_schurs->S,&ksp);CHKERRQ(ierr);
806 ierr = MatSchurComplementSetKSP(S_Ej,ksp);CHKERRQ(ierr);
807 } else { /* build new ksp object which inherits ksp and pc types from the original one */
808 KSP origksp,schurksp;
809 PC origpc,schurpc;
810 KSPType ksp_type;
811 PetscInt n_internal;
812 PetscBool ispcnone;
813
814 ierr = MatSchurComplementGetKSP(sub_schurs->S,&origksp);CHKERRQ(ierr);
815 ierr = MatSchurComplementGetKSP(S_Ej,&schurksp);CHKERRQ(ierr);
816 ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr);
817 ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr);
818 ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr);
819 ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr);
820 ierr = PetscObjectTypeCompare((PetscObject)origpc,PCNONE,&ispcnone);CHKERRQ(ierr);
821 if (!ispcnone) {
822 PCType pc_type;
823 ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr);
824 ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr);
825 } else {
826 ierr = PCSetType(schurpc,PCLU);CHKERRQ(ierr);
827 }
828 ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr);
829 if (!n_internal) { /* UMFPACK gives error with 0 sized problems */
830 MatSolverType solver = NULL;
831 ierr = PCFactorGetMatSolverType(origpc,(MatSolverType*)&solver);CHKERRQ(ierr);
832 if (solver) {
833 ierr = PCFactorSetMatSolverType(schurpc,solver);CHKERRQ(ierr);
834 }
835 }
836 ierr = KSPSetUp(schurksp);CHKERRQ(ierr);
837 }
838 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
839 ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,work,&S_Ej_expl);CHKERRQ(ierr);
840 ierr = PCBDDCComputeExplicitSchur(S_Ej,sub_schurs->is_symmetric,MAT_REUSE_MATRIX,&S_Ej_expl);CHKERRQ(ierr);
841 ierr = PetscObjectTypeCompare((PetscObject)S_Ej_expl,MATSEQDENSE,&Sdense);CHKERRQ(ierr);
842 if (Sdense) {
843 for (j=0;j<subset_size;j++) {
844 dummy_idx[j]=local_size+j;
845 }
846 ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,work,INSERT_VALUES);CHKERRQ(ierr);
847 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet implemented for sparse matrices");
848 ierr = MatDestroy(&S_Ej);CHKERRQ(ierr);
849 ierr = MatDestroy(&S_Ej_expl);CHKERRQ(ierr);
850 local_size += subset_size;
851 }
852 ierr = PetscFree2(dummy_idx,work);CHKERRQ(ierr);
853 /* free */
854 ierr = ISDestroy(&is_I);CHKERRQ(ierr);
855 ierr = MatDestroy(&AE_II);CHKERRQ(ierr);
856 ierr = PetscFree(all_local_idx_N);CHKERRQ(ierr);
857 } else {
858 Mat A,cs_AIB_mat = NULL,benign_AIIm1_ones_mat = NULL;
859 Vec Dall = NULL;
860 IS is_A_all,*is_p_r = NULL;
861 MatType Stype;
862 PetscScalar *work,*S_data,*schur_factor,infty = PETSC_MAX_REAL;
863 PetscScalar *SEj_arr = NULL,*SEjinv_arr = NULL;
864 const PetscScalar *rS_data;
865 PetscInt n,n_I,size_schur,size_active_schur,cum,cum2;
866 PetscBool economic,solver_S,S_lower_triangular = PETSC_FALSE;
867 PetscBool schur_has_vertices,factor_workaround;
868 PetscBool use_cholesky;
869 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
870 PetscBool oldpin;
871 #endif
872
873 /* get sizes */
874 n_I = 0;
875 if (is_I_layer) {
876 ierr = ISGetLocalSize(is_I_layer,&n_I);CHKERRQ(ierr);
877 }
878 economic = PETSC_FALSE;
879 ierr = ISGetLocalSize(sub_schurs->is_I,&cum);CHKERRQ(ierr);
880 if (cum != n_I) economic = PETSC_TRUE;
881 ierr = MatGetLocalSize(sub_schurs->A,&n,NULL);CHKERRQ(ierr);
882 size_active_schur = local_size;
883
884 /* import scaling vector (wrong formulation if we have 3D edges) */
885 if (scaling && compute_Stilda) {
886 const PetscScalar *array;
887 PetscScalar *array2;
888 const PetscInt *idxs;
889 PetscInt i;
890
891 ierr = ISGetIndices(sub_schurs->is_Ej_all,&idxs);CHKERRQ(ierr);
892 ierr = VecCreateSeq(PETSC_COMM_SELF,size_active_schur,&Dall);CHKERRQ(ierr);
893 ierr = VecGetArrayRead(scaling,&array);CHKERRQ(ierr);
894 ierr = VecGetArray(Dall,&array2);CHKERRQ(ierr);
895 for (i=0;i<size_active_schur;i++) array2[i] = array[idxs[i]];
896 ierr = VecRestoreArray(Dall,&array2);CHKERRQ(ierr);
897 ierr = VecRestoreArrayRead(scaling,&array);CHKERRQ(ierr);
898 ierr = ISRestoreIndices(sub_schurs->is_Ej_all,&idxs);CHKERRQ(ierr);
899 deluxe = PETSC_FALSE;
900 }
901
902 /* size active schurs does not count any dirichlet or vertex dof on the interface */
903 factor_workaround = PETSC_FALSE;
904 schur_has_vertices = PETSC_FALSE;
905 cum = n_I+size_active_schur;
906 if (sub_schurs->is_dir) {
907 const PetscInt* idxs;
908 PetscInt n_dir;
909
910 ierr = ISGetLocalSize(sub_schurs->is_dir,&n_dir);CHKERRQ(ierr);
911 ierr = ISGetIndices(sub_schurs->is_dir,&idxs);CHKERRQ(ierr);
912 ierr = PetscArraycpy(all_local_idx_N+cum,idxs,n_dir);CHKERRQ(ierr);
913 ierr = ISRestoreIndices(sub_schurs->is_dir,&idxs);CHKERRQ(ierr);
914 cum += n_dir;
915 factor_workaround = PETSC_TRUE;
916 }
917 /* include the primal vertices in the Schur complement */
918 if (exact_schur && sub_schurs->is_vertices && (compute_Stilda || benign_n)) {
919 PetscInt n_v;
920
921 ierr = ISGetLocalSize(sub_schurs->is_vertices,&n_v);CHKERRQ(ierr);
922 if (n_v) {
923 const PetscInt* idxs;
924
925 ierr = ISGetIndices(sub_schurs->is_vertices,&idxs);CHKERRQ(ierr);
926 ierr = PetscArraycpy(all_local_idx_N+cum,idxs,n_v);CHKERRQ(ierr);
927 ierr = ISRestoreIndices(sub_schurs->is_vertices,&idxs);CHKERRQ(ierr);
928 cum += n_v;
929 factor_workaround = PETSC_TRUE;
930 schur_has_vertices = PETSC_TRUE;
931 }
932 }
933 size_schur = cum - n_I;
934 ierr = ISCreateGeneral(PETSC_COMM_SELF,cum,all_local_idx_N,PETSC_OWN_POINTER,&is_A_all);CHKERRQ(ierr);
935 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
936 oldpin = sub_schurs->A->boundtocpu;
937 ierr = MatBindToCPU(sub_schurs->A,PETSC_TRUE);CHKERRQ(ierr);
938 #endif
939 if (cum == n) {
940 ierr = ISSetPermutation(is_A_all);CHKERRQ(ierr);
941 ierr = MatPermute(sub_schurs->A,is_A_all,is_A_all,&A);CHKERRQ(ierr);
942 } else {
943 ierr = MatCreateSubMatrix(sub_schurs->A,is_A_all,is_A_all,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
944 }
945 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
946 ierr = MatBindToCPU(sub_schurs->A,oldpin);CHKERRQ(ierr);
947 #endif
948 ierr = MatSetOptionsPrefix(A,sub_schurs->prefix);CHKERRQ(ierr);
949 ierr = MatAppendOptionsPrefix(A,"sub_schurs_");CHKERRQ(ierr);
950
951 /* if we actually change the basis for the pressures, LDL^T factors will use a lot of memory
952 this is a workaround */
953 if (benign_n) {
954 Vec v,benign_AIIm1_ones;
955 ISLocalToGlobalMapping N_to_reor;
956 IS is_p0,is_p0_p;
957 PetscScalar *cs_AIB,*AIIm1_data;
958 PetscInt sizeA;
959
960 ierr = ISLocalToGlobalMappingCreateIS(is_A_all,&N_to_reor);CHKERRQ(ierr);
961 ierr = ISCreateGeneral(PETSC_COMM_SELF,benign_n,benign_p0_lidx,PETSC_COPY_VALUES,&is_p0);CHKERRQ(ierr);
962 ierr = ISGlobalToLocalMappingApplyIS(N_to_reor,IS_GTOLM_DROP,is_p0,&is_p0_p);CHKERRQ(ierr);
963 ierr = ISDestroy(&is_p0);CHKERRQ(ierr);
964 ierr = MatCreateVecs(A,&v,&benign_AIIm1_ones);CHKERRQ(ierr);
965 ierr = VecGetSize(v,&sizeA);CHKERRQ(ierr);
966 ierr = MatCreateSeqDense(PETSC_COMM_SELF,sizeA,benign_n,NULL,&benign_AIIm1_ones_mat);CHKERRQ(ierr);
967 ierr = MatCreateSeqDense(PETSC_COMM_SELF,size_schur,benign_n,NULL,&cs_AIB_mat);CHKERRQ(ierr);
968 ierr = MatDenseGetArray(cs_AIB_mat,&cs_AIB);CHKERRQ(ierr);
969 ierr = MatDenseGetArray(benign_AIIm1_ones_mat,&AIIm1_data);CHKERRQ(ierr);
970 ierr = PetscMalloc1(benign_n,&is_p_r);CHKERRQ(ierr);
971 /* compute colsum of A_IB restricted to pressures */
972 for (i=0;i<benign_n;i++) {
973 const PetscScalar *array;
974 const PetscInt *idxs;
975 PetscInt j,nz;
976
977 ierr = ISGlobalToLocalMappingApplyIS(N_to_reor,IS_GTOLM_DROP,benign_zerodiag_subs[i],&is_p_r[i]);CHKERRQ(ierr);
978 ierr = ISGetLocalSize(is_p_r[i],&nz);CHKERRQ(ierr);
979 ierr = ISGetIndices(is_p_r[i],&idxs);CHKERRQ(ierr);
980 for (j=0;j<nz;j++) AIIm1_data[idxs[j]+sizeA*i] = 1.;
981 ierr = ISRestoreIndices(is_p_r[i],&idxs);CHKERRQ(ierr);
982 ierr = VecPlaceArray(benign_AIIm1_ones,AIIm1_data+sizeA*i);CHKERRQ(ierr);
983 ierr = MatMult(A,benign_AIIm1_ones,v);CHKERRQ(ierr);
984 ierr = VecResetArray(benign_AIIm1_ones);CHKERRQ(ierr);
985 ierr = VecGetArrayRead(v,&array);CHKERRQ(ierr);
986 for (j=0;j<size_schur;j++) {
987 #if defined(PETSC_USE_COMPLEX)
988 cs_AIB[i*size_schur + j] = (PetscRealPart(array[j+n_I])/nz + PETSC_i*(PetscImaginaryPart(array[j+n_I])/nz));
989 #else
990 cs_AIB[i*size_schur + j] = array[j+n_I]/nz;
991 #endif
992 }
993 ierr = VecRestoreArrayRead(v,&array);CHKERRQ(ierr);
994 }
995 ierr = MatDenseRestoreArray(cs_AIB_mat,&cs_AIB);CHKERRQ(ierr);
996 ierr = MatDenseRestoreArray(benign_AIIm1_ones_mat,&AIIm1_data);CHKERRQ(ierr);
997 ierr = VecDestroy(&v);CHKERRQ(ierr);
998 ierr = VecDestroy(&benign_AIIm1_ones);CHKERRQ(ierr);
999 ierr = MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE);CHKERRQ(ierr);
1000 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
1001 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
1002 ierr = MatZeroRowsColumnsIS(A,is_p0_p,1.0,NULL,NULL);CHKERRQ(ierr);
1003 ierr = ISDestroy(&is_p0_p);CHKERRQ(ierr);
1004 ierr = ISLocalToGlobalMappingDestroy(&N_to_reor);CHKERRQ(ierr);
1005 }
1006 ierr = MatSetOption(A,MAT_SYMMETRIC,sub_schurs->is_symmetric);CHKERRQ(ierr);
1007 ierr = MatSetOption(A,MAT_HERMITIAN,sub_schurs->is_hermitian);CHKERRQ(ierr);
1008 ierr = MatSetOption(A,MAT_SPD,sub_schurs->is_posdef);CHKERRQ(ierr);
1009
1010 /* for complexes, symmetric and hermitian at the same time implies null imaginary part */
1011 use_cholesky = (PetscBool)((use_potr || use_sytr) && sub_schurs->is_hermitian && sub_schurs->is_symmetric);
1012
1013 /* when using the benign subspace trick, the local Schur complements are SPD */
1014 if (benign_trick) sub_schurs->is_posdef = PETSC_TRUE;
1015
1016 if (n_I) {
1017 IS is_schur;
1018 char stype[64];
1019 PetscBool gpu = PETSC_FALSE;
1020
1021 if (use_cholesky) {
1022 ierr = MatGetFactor(A,sub_schurs->mat_solver_type,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
1023 } else {
1024 ierr = MatGetFactor(A,sub_schurs->mat_solver_type,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
1025 }
1026 ierr = MatSetErrorIfFailure(A,PETSC_TRUE);CHKERRQ(ierr);
1027
1028 /* subsets ordered last */
1029 ierr = ISCreateStride(PETSC_COMM_SELF,size_schur,n_I,1,&is_schur);CHKERRQ(ierr);
1030 ierr = MatFactorSetSchurIS(F,is_schur);CHKERRQ(ierr);
1031 ierr = ISDestroy(&is_schur);CHKERRQ(ierr);
1032
1033 /* factorization step */
1034 if (use_cholesky) {
1035 ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr);
1036 #if defined(PETSC_HAVE_MUMPS) /* be sure that icntl 19 is not set by command line */
1037 ierr = MatMumpsSetIcntl(F,19,2);CHKERRQ(ierr);
1038 #endif
1039 ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr);
1040 S_lower_triangular = PETSC_TRUE;
1041 } else {
1042 ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr);
1043 #if defined(PETSC_HAVE_MUMPS) /* be sure that icntl 19 is not set by command line */
1044 ierr = MatMumpsSetIcntl(F,19,3);CHKERRQ(ierr);
1045 #endif
1046 ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr);
1047 }
1048 ierr = MatViewFromOptions(F,(PetscObject)A,"-mat_factor_view");CHKERRQ(ierr);
1049
1050 if (matl_dbg_viewer) {
1051 Mat S;
1052 IS is;
1053
1054 ierr = PetscObjectSetName((PetscObject)A,"A");CHKERRQ(ierr);
1055 ierr = MatView(A,matl_dbg_viewer);CHKERRQ(ierr);
1056 ierr = MatFactorCreateSchurComplement(F,&S,NULL);CHKERRQ(ierr);
1057 ierr = PetscObjectSetName((PetscObject)S,"S");CHKERRQ(ierr);
1058 ierr = MatView(S,matl_dbg_viewer);CHKERRQ(ierr);
1059 ierr = MatDestroy(&S);CHKERRQ(ierr);
1060 ierr = ISCreateStride(PETSC_COMM_SELF,n_I,0,1,&is);CHKERRQ(ierr);
1061 ierr = PetscObjectSetName((PetscObject)is,"I");CHKERRQ(ierr);
1062 ierr = ISView(is,matl_dbg_viewer);CHKERRQ(ierr);
1063 ierr = ISDestroy(&is);CHKERRQ(ierr);
1064 ierr = ISCreateStride(PETSC_COMM_SELF,size_schur,n_I,1,&is);CHKERRQ(ierr);
1065 ierr = PetscObjectSetName((PetscObject)is,"B");CHKERRQ(ierr);
1066 ierr = ISView(is,matl_dbg_viewer);CHKERRQ(ierr);
1067 ierr = ISDestroy(&is);CHKERRQ(ierr);
1068 ierr = PetscObjectSetName((PetscObject)is_A_all,"IA");CHKERRQ(ierr);
1069 ierr = ISView(is_A_all,matl_dbg_viewer);CHKERRQ(ierr);
1070 }
1071
1072 /* get explicit Schur Complement computed during numeric factorization */
1073 ierr = MatFactorGetSchurComplement(F,&S_all,NULL);CHKERRQ(ierr);
1074 ierr = PetscStrncpy(stype,MATSEQDENSE,sizeof(stype));CHKERRQ(ierr);
1075 #if defined(PETSC_HAVE_CUDA)
1076 ierr = PetscObjectTypeCompareAny((PetscObject)A,&gpu,MATSEQAIJVIENNACL,MATSEQAIJCUSPARSE,"");CHKERRQ(ierr);
1077 #endif
1078 if (gpu) {
1079 ierr = PetscStrncpy(stype,MATSEQDENSECUDA,sizeof(stype));CHKERRQ(ierr);
1080 }
1081 ierr = PetscOptionsGetString(NULL,sub_schurs->prefix,"-sub_schurs_schur_mat_type",stype,sizeof(stype),NULL);CHKERRQ(ierr);
1082 ierr = MatConvert(S_all,stype,MAT_INPLACE_MATRIX,&S_all);CHKERRQ(ierr);
1083 ierr = MatSetOption(S_all,MAT_SPD,sub_schurs->is_posdef);CHKERRQ(ierr);
1084 ierr = MatSetOption(S_all,MAT_HERMITIAN,sub_schurs->is_hermitian);CHKERRQ(ierr);
1085 ierr = MatGetType(S_all,&Stype);CHKERRQ(ierr);
1086
1087 /* we can reuse the solvers if we are not using the economic version */
1088 reuse_solvers = (PetscBool)(reuse_solvers && !economic);
1089 factor_workaround = (PetscBool)(reuse_solvers && factor_workaround);
1090 if (!sub_schurs->is_posdef && factor_workaround && compute_Stilda && size_active_schur)
1091 reuse_solvers = factor_workaround = PETSC_FALSE;
1092
1093 solver_S = PETSC_TRUE;
1094
1095 /* update the Schur complement with the change of basis on the pressures */
1096 if (benign_n) {
1097 const PetscScalar *cs_AIB;
1098 PetscScalar *S_data,*AIIm1_data;
1099 Mat S2 = NULL,S3 = NULL; /* dbg */
1100 PetscScalar *S2_data,*S3_data; /* dbg */
1101 Vec v,benign_AIIm1_ones;
1102 PetscInt sizeA;
1103
1104 ierr = MatDenseGetArray(S_all,&S_data);CHKERRQ(ierr);
1105 ierr = MatCreateVecs(A,&v,&benign_AIIm1_ones);CHKERRQ(ierr);
1106 ierr = VecGetSize(v,&sizeA);CHKERRQ(ierr);
1107 #if defined(PETSC_HAVE_MUMPS)
1108 ierr = MatMumpsSetIcntl(F,26,0);CHKERRQ(ierr);
1109 #endif
1110 #if defined(PETSC_HAVE_MKL_PARDISO)
1111 ierr = MatMkl_PardisoSetCntl(F,70,1);CHKERRQ(ierr);
1112 #endif
1113 ierr = MatDenseGetArrayRead(cs_AIB_mat,&cs_AIB);CHKERRQ(ierr);
1114 ierr = MatDenseGetArray(benign_AIIm1_ones_mat,&AIIm1_data);CHKERRQ(ierr);
1115 if (matl_dbg_viewer) {
1116 ierr = MatDuplicate(S_all,MAT_DO_NOT_COPY_VALUES,&S2);CHKERRQ(ierr);
1117 ierr = MatDuplicate(S_all,MAT_DO_NOT_COPY_VALUES,&S3);CHKERRQ(ierr);
1118 ierr = MatDenseGetArray(S2,&S2_data);CHKERRQ(ierr);
1119 ierr = MatDenseGetArray(S3,&S3_data);CHKERRQ(ierr);
1120 }
1121 for (i=0;i<benign_n;i++) {
1122 PetscScalar *array,sum = 0.,one = 1.,*sums;
1123 const PetscInt *idxs;
1124 PetscInt k,j,nz;
1125 PetscBLASInt B_k,B_n;
1126
1127 ierr = PetscCalloc1(benign_n,&sums);CHKERRQ(ierr);
1128 ierr = VecPlaceArray(benign_AIIm1_ones,AIIm1_data+sizeA*i);CHKERRQ(ierr);
1129 ierr = VecCopy(benign_AIIm1_ones,v);CHKERRQ(ierr);
1130 ierr = MatSolve(F,v,benign_AIIm1_ones);CHKERRQ(ierr);
1131 ierr = MatMult(A,benign_AIIm1_ones,v);CHKERRQ(ierr);
1132 ierr = VecResetArray(benign_AIIm1_ones);CHKERRQ(ierr);
1133 /* p0 dofs (eliminated) are excluded from the sums */
1134 for (k=0;k<benign_n;k++) {
1135 ierr = ISGetLocalSize(is_p_r[k],&nz);CHKERRQ(ierr);
1136 ierr = ISGetIndices(is_p_r[k],&idxs);CHKERRQ(ierr);
1137 for (j=0;j<nz-1;j++) sums[k] -= AIIm1_data[idxs[j]+sizeA*i];
1138 ierr = ISRestoreIndices(is_p_r[k],&idxs);CHKERRQ(ierr);
1139 }
1140 ierr = VecGetArrayRead(v,(const PetscScalar**)&array);CHKERRQ(ierr);
1141 if (matl_dbg_viewer) {
1142 Vec vv;
1143 char name[16];
1144
1145 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,size_schur,array+n_I,&vv);CHKERRQ(ierr);
1146 ierr = PetscSNPrintf(name,sizeof(name),"Pvs%D",i);CHKERRQ(ierr);
1147 ierr = PetscObjectSetName((PetscObject)vv,name);CHKERRQ(ierr);
1148 ierr = VecView(vv,matl_dbg_viewer);CHKERRQ(ierr);
1149 }
1150 /* perform sparse rank updates on symmetric Schur (TODO: move outside of the loop?) */
1151 /* cs_AIB already scaled by 1./nz */
1152 B_k = 1;
1153 for (k=0;k<benign_n;k++) {
1154 sum = sums[k];
1155 ierr = PetscBLASIntCast(size_schur,&B_n);CHKERRQ(ierr);
1156
1157 if (PetscAbsScalar(sum) == 0.0) continue;
1158 if (k == i) {
1159 PetscStackCallBLAS("BLASsyrk",BLASsyrk_("L","N",&B_n,&B_k,&sum,cs_AIB+i*size_schur,&B_n,&one,S_data,&B_n));
1160 if (matl_dbg_viewer) {
1161 PetscStackCallBLAS("BLASsyrk",BLASsyrk_("L","N",&B_n,&B_k,&sum,cs_AIB+i*size_schur,&B_n,&one,S3_data,&B_n));
1162 }
1163 } else { /* XXX Is it correct to use symmetric rank-2 update with half of the sum? */
1164 sum /= 2.0;
1165 PetscStackCallBLAS("BLASsyr2k",BLASsyr2k_("L","N",&B_n,&B_k,&sum,cs_AIB+k*size_schur,&B_n,cs_AIB+i*size_schur,&B_n,&one,S_data,&B_n));
1166 if (matl_dbg_viewer) {
1167 PetscStackCallBLAS("BLASsyr2k",BLASsyr2k_("L","N",&B_n,&B_k,&sum,cs_AIB+k*size_schur,&B_n,cs_AIB+i*size_schur,&B_n,&one,S3_data,&B_n));
1168 }
1169 }
1170 }
1171 sum = 1.;
1172 PetscStackCallBLAS("BLASsyr2k",BLASsyr2k_("L","N",&B_n,&B_k,&sum,array+n_I,&B_n,cs_AIB+i*size_schur,&B_n,&one,S_data,&B_n));
1173 if (matl_dbg_viewer) {
1174 PetscStackCallBLAS("BLASsyr2k",BLASsyr2k_("L","N",&B_n,&B_k,&sum,array+n_I,&B_n,cs_AIB+i*size_schur,&B_n,&one,S2_data,&B_n));
1175 }
1176 ierr = VecRestoreArrayRead(v,(const PetscScalar**)&array);CHKERRQ(ierr);
1177 /* set p0 entry of AIIm1_ones to zero */
1178 ierr = ISGetLocalSize(is_p_r[i],&nz);CHKERRQ(ierr);
1179 ierr = ISGetIndices(is_p_r[i],&idxs);CHKERRQ(ierr);
1180 for (j=0;j<benign_n;j++) AIIm1_data[idxs[nz-1]+sizeA*j] = 0.;
1181 ierr = ISRestoreIndices(is_p_r[i],&idxs);CHKERRQ(ierr);
1182 ierr = PetscFree(sums);CHKERRQ(ierr);
1183 }
1184 ierr = VecDestroy(&benign_AIIm1_ones);CHKERRQ(ierr);
1185 if (matl_dbg_viewer) {
1186 ierr = MatDenseRestoreArray(S2,&S2_data);CHKERRQ(ierr);
1187 ierr = MatDenseRestoreArray(S3,&S3_data);CHKERRQ(ierr);
1188 }
1189 if (!S_lower_triangular) { /* I need to expand the upper triangular data (column oriented) */
1190 PetscInt k,j;
1191 for (k=0;k<size_schur;k++) {
1192 for (j=k;j<size_schur;j++) {
1193 S_data[j*size_schur+k] = PetscConj(S_data[k*size_schur+j]);
1194 }
1195 }
1196 }
1197
1198 /* restore defaults */
1199 #if defined(PETSC_HAVE_MUMPS)
1200 ierr = MatMumpsSetIcntl(F,26,-1);CHKERRQ(ierr);
1201 #endif
1202 #if defined(PETSC_HAVE_MKL_PARDISO)
1203 ierr = MatMkl_PardisoSetCntl(F,70,0);CHKERRQ(ierr);
1204 #endif
1205 ierr = MatDenseRestoreArrayRead(cs_AIB_mat,&cs_AIB);CHKERRQ(ierr);
1206 ierr = MatDenseRestoreArray(benign_AIIm1_ones_mat,&AIIm1_data);CHKERRQ(ierr);
1207 ierr = VecDestroy(&v);CHKERRQ(ierr);
1208 ierr = MatDenseRestoreArray(S_all,&S_data);CHKERRQ(ierr);
1209 if (matl_dbg_viewer) {
1210 Mat S;
1211
1212 ierr = MatFactorRestoreSchurComplement(F,&S_all,MAT_FACTOR_SCHUR_UNFACTORED);CHKERRQ(ierr);
1213 ierr = MatFactorCreateSchurComplement(F,&S,NULL);CHKERRQ(ierr);
1214 ierr = PetscObjectSetName((PetscObject)S,"Sb");CHKERRQ(ierr);
1215 ierr = MatView(S,matl_dbg_viewer);CHKERRQ(ierr);
1216 ierr = MatDestroy(&S);CHKERRQ(ierr);
1217 ierr = PetscObjectSetName((PetscObject)S2,"S2P");CHKERRQ(ierr);
1218 ierr = MatView(S2,matl_dbg_viewer);CHKERRQ(ierr);
1219 ierr = PetscObjectSetName((PetscObject)S3,"S3P");CHKERRQ(ierr);
1220 ierr = MatView(S3,matl_dbg_viewer);CHKERRQ(ierr);
1221 ierr = PetscObjectSetName((PetscObject)cs_AIB_mat,"cs");CHKERRQ(ierr);
1222 ierr = MatView(cs_AIB_mat,matl_dbg_viewer);CHKERRQ(ierr);
1223 ierr = MatFactorGetSchurComplement(F,&S_all,NULL);CHKERRQ(ierr);
1224 }
1225 ierr = MatDestroy(&S2);CHKERRQ(ierr);
1226 ierr = MatDestroy(&S3);CHKERRQ(ierr);
1227 }
1228 if (!reuse_solvers) {
1229 for (i=0;i<benign_n;i++) {
1230 ierr = ISDestroy(&is_p_r[i]);CHKERRQ(ierr);
1231 }
1232 ierr = PetscFree(is_p_r);CHKERRQ(ierr);
1233 ierr = MatDestroy(&cs_AIB_mat);CHKERRQ(ierr);
1234 ierr = MatDestroy(&benign_AIIm1_ones_mat);CHKERRQ(ierr);
1235 }
1236 } else { /* we can't use MatFactor when size_schur == size_of_the_problem */
1237 ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&S_all);CHKERRQ(ierr);
1238 ierr = MatGetType(S_all,&Stype);CHKERRQ(ierr);
1239 reuse_solvers = PETSC_FALSE; /* TODO: why we can't reuse the solvers here? */
1240 factor_workaround = PETSC_FALSE;
1241 solver_S = PETSC_FALSE;
1242 }
1243
1244 if (reuse_solvers) {
1245 Mat A_II,Afake;
1246 Vec vec1_B;
1247 PCBDDCReuseSolvers msolv_ctx;
1248 PetscInt n_R;
1249
1250 if (sub_schurs->reuse_solver) {
1251 ierr = PCBDDCReuseSolversReset(sub_schurs->reuse_solver);CHKERRQ(ierr);
1252 } else {
1253 ierr = PetscNew(&sub_schurs->reuse_solver);CHKERRQ(ierr);
1254 }
1255 msolv_ctx = sub_schurs->reuse_solver;
1256 ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1257 ierr = PetscObjectReference((PetscObject)F);CHKERRQ(ierr);
1258 msolv_ctx->F = F;
1259 ierr = MatCreateVecs(F,&msolv_ctx->sol,NULL);CHKERRQ(ierr);
1260 /* currently PETSc has no support for MatSolve(F,x,x), so cheat and let rhs and sol share the same memory */
1261 {
1262 PetscScalar *array;
1263 PetscInt n;
1264
1265 ierr = VecGetLocalSize(msolv_ctx->sol,&n);CHKERRQ(ierr);
1266 ierr = VecGetArray(msolv_ctx->sol,&array);CHKERRQ(ierr);
1267 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)msolv_ctx->sol),1,n,array,&msolv_ctx->rhs);CHKERRQ(ierr);
1268 ierr = VecRestoreArray(msolv_ctx->sol,&array);CHKERRQ(ierr);
1269 }
1270 msolv_ctx->has_vertices = schur_has_vertices;
1271
1272 /* interior solver */
1273 ierr = PCCreate(PetscObjectComm((PetscObject)A_II),&msolv_ctx->interior_solver);CHKERRQ(ierr);
1274 ierr = PCSetOperators(msolv_ctx->interior_solver,A_II,A_II);CHKERRQ(ierr);
1275 ierr = PCSetType(msolv_ctx->interior_solver,PCSHELL);CHKERRQ(ierr);
1276 ierr = PCShellSetName(msolv_ctx->interior_solver,"Interior solver (w/o Schur factorization)");CHKERRQ(ierr);
1277 ierr = PCShellSetContext(msolv_ctx->interior_solver,msolv_ctx);CHKERRQ(ierr);
1278 ierr = PCShellSetView(msolv_ctx->interior_solver,PCBDDCReuseSolvers_View);CHKERRQ(ierr);
1279 ierr = PCShellSetApply(msolv_ctx->interior_solver,PCBDDCReuseSolvers_Interior);CHKERRQ(ierr);
1280 ierr = PCShellSetApplyTranspose(msolv_ctx->interior_solver,PCBDDCReuseSolvers_InteriorTranspose);CHKERRQ(ierr);
1281
1282 /* correction solver */
1283 ierr = PCCreate(PetscObjectComm((PetscObject)A_II),&msolv_ctx->correction_solver);CHKERRQ(ierr);
1284 ierr = PCSetType(msolv_ctx->correction_solver,PCSHELL);CHKERRQ(ierr);
1285 ierr = PCShellSetName(msolv_ctx->correction_solver,"Correction solver (with Schur factorization)");CHKERRQ(ierr);
1286 ierr = PCShellSetContext(msolv_ctx->correction_solver,msolv_ctx);CHKERRQ(ierr);
1287 ierr = PCShellSetView(msolv_ctx->interior_solver,PCBDDCReuseSolvers_View);CHKERRQ(ierr);
1288 ierr = PCShellSetApply(msolv_ctx->correction_solver,PCBDDCReuseSolvers_Correction);CHKERRQ(ierr);
1289 ierr = PCShellSetApplyTranspose(msolv_ctx->correction_solver,PCBDDCReuseSolvers_CorrectionTranspose);CHKERRQ(ierr);
1290
1291 /* scatter and vecs for Schur complement solver */
1292 ierr = MatCreateVecs(S_all,&msolv_ctx->sol_B,&msolv_ctx->rhs_B);CHKERRQ(ierr);
1293 ierr = MatCreateVecs(sub_schurs->S,&vec1_B,NULL);CHKERRQ(ierr);
1294 if (!schur_has_vertices) {
1295 ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,is_A_all,&msolv_ctx->is_B);CHKERRQ(ierr);
1296 ierr = VecScatterCreate(vec1_B,msolv_ctx->is_B,msolv_ctx->sol_B,NULL,&msolv_ctx->correction_scatter_B);CHKERRQ(ierr);
1297 ierr = PetscObjectReference((PetscObject)is_A_all);CHKERRQ(ierr);
1298 msolv_ctx->is_R = is_A_all;
1299 } else {
1300 IS is_B_all;
1301 const PetscInt* idxs;
1302 PetscInt dual,n_v,n;
1303
1304 ierr = ISGetLocalSize(sub_schurs->is_vertices,&n_v);CHKERRQ(ierr);
1305 dual = size_schur - n_v;
1306 ierr = ISGetLocalSize(is_A_all,&n);CHKERRQ(ierr);
1307 ierr = ISGetIndices(is_A_all,&idxs);CHKERRQ(ierr);
1308 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is_A_all),dual,idxs+n_I,PETSC_COPY_VALUES,&is_B_all);CHKERRQ(ierr);
1309 ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,is_B_all,&msolv_ctx->is_B);CHKERRQ(ierr);
1310 ierr = ISDestroy(&is_B_all);CHKERRQ(ierr);
1311 ierr = ISCreateStride(PetscObjectComm((PetscObject)is_A_all),dual,0,1,&is_B_all);CHKERRQ(ierr);
1312 ierr = VecScatterCreate(vec1_B,msolv_ctx->is_B,msolv_ctx->sol_B,is_B_all,&msolv_ctx->correction_scatter_B);CHKERRQ(ierr);
1313 ierr = ISDestroy(&is_B_all);CHKERRQ(ierr);
1314 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is_A_all),n-n_v,idxs,PETSC_COPY_VALUES,&msolv_ctx->is_R);CHKERRQ(ierr);
1315 ierr = ISRestoreIndices(is_A_all,&idxs);CHKERRQ(ierr);
1316 }
1317 ierr = ISGetLocalSize(msolv_ctx->is_R,&n_R);CHKERRQ(ierr);
1318 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,n_R,n_R,0,NULL,&Afake);CHKERRQ(ierr);
1319 ierr = MatAssemblyBegin(Afake,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1320 ierr = MatAssemblyEnd(Afake,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1321 ierr = PCSetOperators(msolv_ctx->correction_solver,Afake,Afake);CHKERRQ(ierr);
1322 ierr = MatDestroy(&Afake);CHKERRQ(ierr);
1323 ierr = VecDestroy(&vec1_B);CHKERRQ(ierr);
1324
1325 /* communicate benign info to solver context */
1326 if (benign_n) {
1327 PetscScalar *array;
1328
1329 msolv_ctx->benign_n = benign_n;
1330 msolv_ctx->benign_zerodiag_subs = is_p_r;
1331 ierr = PetscMalloc1(benign_n,&msolv_ctx->benign_save_vals);CHKERRQ(ierr);
1332 msolv_ctx->benign_csAIB = cs_AIB_mat;
1333 ierr = MatCreateVecs(cs_AIB_mat,&msolv_ctx->benign_corr_work,NULL);CHKERRQ(ierr);
1334 ierr = VecGetArray(msolv_ctx->benign_corr_work,&array);CHKERRQ(ierr);
1335 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,size_schur,array,&msolv_ctx->benign_dummy_schur_vec);CHKERRQ(ierr);
1336 ierr = VecRestoreArray(msolv_ctx->benign_corr_work,&array);CHKERRQ(ierr);
1337 msolv_ctx->benign_AIIm1ones = benign_AIIm1_ones_mat;
1338 }
1339 } else {
1340 if (sub_schurs->reuse_solver) {
1341 ierr = PCBDDCReuseSolversReset(sub_schurs->reuse_solver);CHKERRQ(ierr);
1342 }
1343 ierr = PetscFree(sub_schurs->reuse_solver);CHKERRQ(ierr);
1344 }
1345 ierr = MatDestroy(&A);CHKERRQ(ierr);
1346 ierr = ISDestroy(&is_A_all);CHKERRQ(ierr);
1347
1348 /* Work arrays */
1349 ierr = PetscMalloc1(max_subset_size*max_subset_size,&work);CHKERRQ(ierr);
1350
1351 /* S_Ej_all */
1352 cum = cum2 = 0;
1353 ierr = MatDenseGetArrayRead(S_all,&rS_data);CHKERRQ(ierr);
1354 ierr = MatSeqAIJGetArray(sub_schurs->S_Ej_all,&SEj_arr);CHKERRQ(ierr);
1355 if (compute_Stilda) {
1356 ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_inv_all,&SEjinv_arr);CHKERRQ(ierr);
1357 }
1358 for (i=0;i<sub_schurs->n_subs;i++) {
1359 PetscInt j;
1360
1361 /* get S_E */
1362 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
1363 if (S_lower_triangular) { /* I need to expand the upper triangular data (column oriented) */
1364 PetscInt k;
1365 for (k=0;k<subset_size;k++) {
1366 for (j=k;j<subset_size;j++) {
1367 work[k*subset_size+j] = rS_data[cum2+k*size_schur+j];
1368 work[j*subset_size+k] = PetscConj(rS_data[cum2+k*size_schur+j]);
1369 }
1370 }
1371 } else { /* just copy to workspace */
1372 PetscInt k;
1373 for (k=0;k<subset_size;k++) {
1374 for (j=0;j<subset_size;j++) {
1375 work[k*subset_size+j] = rS_data[cum2+k*size_schur+j];
1376 }
1377 }
1378 }
1379 /* insert S_E values */
1380 if (sub_schurs->change) {
1381 Mat change_sub,SEj,T;
1382
1383 /* change basis */
1384 ierr = KSPGetOperators(sub_schurs->change[i],&change_sub,NULL);CHKERRQ(ierr);
1385 ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,work,&SEj);CHKERRQ(ierr);
1386 if (!sub_schurs->change_with_qr) { /* currently there's no support for PtAP with P SeqAIJ */
1387 Mat T2;
1388 ierr = MatTransposeMatMult(change_sub,SEj,MAT_INITIAL_MATRIX,1.0,&T2);CHKERRQ(ierr);
1389 ierr = MatMatMult(T2,change_sub,MAT_INITIAL_MATRIX,1.0,&T);CHKERRQ(ierr);
1390 ierr = MatConvert(T,MATSEQDENSE,MAT_INPLACE_MATRIX,&T);CHKERRQ(ierr);
1391 ierr = MatDestroy(&T2);CHKERRQ(ierr);
1392 } else {
1393 ierr = MatPtAP(SEj,change_sub,MAT_INITIAL_MATRIX,1.0,&T);CHKERRQ(ierr);
1394 }
1395 ierr = MatCopy(T,SEj,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
1396 ierr = MatDestroy(&T);CHKERRQ(ierr);
1397 ierr = MatZeroRowsColumnsIS(SEj,sub_schurs->change_primal_sub[i],1.0,NULL,NULL);CHKERRQ(ierr);
1398 ierr = MatDestroy(&SEj);CHKERRQ(ierr);
1399 }
1400 if (deluxe) {
1401 ierr = PetscArraycpy(SEj_arr,work,subset_size*subset_size);CHKERRQ(ierr);
1402 /* if adaptivity is requested, invert S_E blocks */
1403 if (compute_Stilda) {
1404 Mat M;
1405 const PetscScalar *vals;
1406 PetscBool isdense,isdensecuda;
1407
1408 ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,work,&M);CHKERRQ(ierr);
1409 ierr = MatSetOption(M,MAT_SPD,sub_schurs->is_posdef);CHKERRQ(ierr);
1410 ierr = MatSetOption(M,MAT_HERMITIAN,sub_schurs->is_hermitian);CHKERRQ(ierr);
1411 if (!PetscBTLookup(sub_schurs->is_edge,i)) {
1412 ierr = MatSetType(M,Stype);CHKERRQ(ierr);
1413 }
1414 ierr = PetscObjectTypeCompare((PetscObject)M,MATSEQDENSE,&isdense);CHKERRQ(ierr);
1415 ierr = PetscObjectTypeCompare((PetscObject)M,MATSEQDENSECUDA,&isdensecuda);CHKERRQ(ierr);
1416 if (use_cholesky) {
1417 ierr = MatCholeskyFactor(M,NULL,NULL);CHKERRQ(ierr);
1418 } else {
1419 ierr = MatLUFactor(M,NULL,NULL,NULL);CHKERRQ(ierr);
1420 }
1421 if (isdense) {
1422 ierr = MatSeqDenseInvertFactors_Private(M);CHKERRQ(ierr);
1423 #if defined(PETSC_HAVE_CUDA)
1424 } else if (isdensecuda) {
1425 ierr = MatSeqDenseCUDAInvertFactors_Private(M);CHKERRQ(ierr);
1426 #endif
1427 } else SETERRQ1(PetscObjectComm((PetscObject)M),PETSC_ERR_SUP,"Not implemented for type %s",Stype);
1428 ierr = MatDenseGetArrayRead(M,&vals);CHKERRQ(ierr);
1429 ierr = PetscArraycpy(SEjinv_arr,vals,subset_size*subset_size);CHKERRQ(ierr);
1430 ierr = MatDenseRestoreArrayRead(M,&vals);CHKERRQ(ierr);
1431 ierr = MatDestroy(&M);CHKERRQ(ierr);
1432 }
1433 } else if (compute_Stilda) { /* not using deluxe */
1434 Mat SEj;
1435 Vec D;
1436 PetscScalar *array;
1437
1438 ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,work,&SEj);CHKERRQ(ierr);
1439 ierr = VecGetArray(Dall,&array);CHKERRQ(ierr);
1440 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,subset_size,array+cum,&D);CHKERRQ(ierr);
1441 ierr = VecRestoreArray(Dall,&array);CHKERRQ(ierr);
1442 ierr = VecShift(D,-1.);CHKERRQ(ierr);
1443 ierr = MatDiagonalScale(SEj,D,D);CHKERRQ(ierr);
1444 ierr = MatDestroy(&SEj);CHKERRQ(ierr);
1445 ierr = VecDestroy(&D);CHKERRQ(ierr);
1446 ierr = PetscArraycpy(SEj_arr,work,subset_size*subset_size);CHKERRQ(ierr);
1447 }
1448 cum += subset_size;
1449 cum2 += subset_size*(size_schur + 1);
1450 SEj_arr += subset_size*subset_size;
1451 if (SEjinv_arr) SEjinv_arr += subset_size*subset_size;
1452 }
1453 ierr = MatDenseRestoreArrayRead(S_all,&rS_data);CHKERRQ(ierr);
1454 ierr = MatSeqAIJRestoreArray(sub_schurs->S_Ej_all,&SEj_arr);CHKERRQ(ierr);
1455 if (compute_Stilda) {
1456 ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_inv_all,&SEjinv_arr);CHKERRQ(ierr);
1457 }
1458 if (solver_S) {
1459 ierr = MatFactorRestoreSchurComplement(F,&S_all,MAT_FACTOR_SCHUR_UNFACTORED);CHKERRQ(ierr);
1460 }
1461
1462 /* may prevent from unneeded copies, since MUMPS or MKL_Pardiso always use CPU memory
1463 however, preliminary tests indicate using GPUs is still faster in the solve phase */
1464 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1465 if (reuse_solvers) {
1466 Mat St;
1467 MatFactorSchurStatus st;
1468 PetscBool flg = PETSC_FALSE;
1469
1470 ierr = PetscOptionsGetBool(NULL,sub_schurs->prefix,"-sub_schurs_schur_pin_to_cpu",&flg,NULL);CHKERRQ(ierr);
1471 ierr = MatFactorGetSchurComplement(F,&St,&st);CHKERRQ(ierr);
1472 ierr = MatBindToCPU(St,flg);CHKERRQ(ierr);
1473 ierr = MatFactorRestoreSchurComplement(F,&St,st);CHKERRQ(ierr);
1474 }
1475 #endif
1476
1477 schur_factor = NULL;
1478 if (compute_Stilda && size_active_schur) {
1479
1480 ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all,&SEjinv_arr);CHKERRQ(ierr);
1481 if (sub_schurs->n_subs == 1 && size_schur == size_active_schur && deluxe) { /* we already computed the inverse */
1482 ierr = PetscArraycpy(SEjinv_arr,work,size_schur*size_schur);CHKERRQ(ierr);
1483 } else {
1484 Mat S_all_inv=NULL;
1485
1486 if (solver_S) {
1487 /* for adaptive selection we need S^-1; for solver reusage we need S_\Delta\Delta^-1.
1488 The latter is not the principal subminor for S^-1. However, the factors can be reused since S_\Delta\Delta is the leading principal submatrix of S */
1489 if (factor_workaround) {/* invert without calling MatFactorInvertSchurComplement, since we are hacking */
1490 PetscScalar *data;
1491 PetscInt nd = 0;
1492
1493 if (!use_potr) {
1494 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor update not yet implemented for non SPD matrices");
1495 }
1496 ierr = MatFactorGetSchurComplement(F,&S_all_inv,NULL);CHKERRQ(ierr);
1497 ierr = MatDenseGetArray(S_all_inv,&data);CHKERRQ(ierr);
1498 if (sub_schurs->is_dir) { /* dirichlet dofs could have different scalings */
1499 ierr = ISGetLocalSize(sub_schurs->is_dir,&nd);CHKERRQ(ierr);
1500 }
1501
1502 /* factor and invert activedofs and vertices (dirichlet dofs does not contribute) */
1503 if (schur_has_vertices) {
1504 Mat M;
1505 PetscScalar *tdata;
1506 PetscInt nv = 0, news;
1507
1508 ierr = ISGetLocalSize(sub_schurs->is_vertices,&nv);CHKERRQ(ierr);
1509 news = size_active_schur + nv;
1510 ierr = PetscCalloc1(news*news,&tdata);CHKERRQ(ierr);
1511 for (i=0;i<size_active_schur;i++) {
1512 ierr = PetscArraycpy(tdata+i*(news+1),data+i*(size_schur+1),size_active_schur-i);CHKERRQ(ierr);
1513 ierr = PetscArraycpy(tdata+i*(news+1)+size_active_schur-i,data+i*size_schur+size_active_schur+nd,nv);CHKERRQ(ierr);
1514 }
1515 for (i=0;i<nv;i++) {
1516 PetscInt k = i+size_active_schur;
1517 ierr = PetscArraycpy(tdata+k*(news+1),data+(k+nd)*(size_schur+1),nv-i);CHKERRQ(ierr);
1518 }
1519
1520 ierr = MatCreateSeqDense(PETSC_COMM_SELF,news,news,tdata,&M);CHKERRQ(ierr);
1521 ierr = MatSetOption(M,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
1522 ierr = MatCholeskyFactor(M,NULL,NULL);CHKERRQ(ierr);
1523 /* save the factors */
1524 cum = 0;
1525 ierr = PetscMalloc1((size_active_schur*(size_active_schur +1))/2+nd,&schur_factor);CHKERRQ(ierr);
1526 for (i=0;i<size_active_schur;i++) {
1527 ierr = PetscArraycpy(schur_factor+cum,tdata+i*(news+1),size_active_schur-i);CHKERRQ(ierr);
1528 cum += size_active_schur - i;
1529 }
1530 for (i=0;i<nd;i++) schur_factor[cum+i] = PetscSqrtReal(PetscRealPart(data[(i+size_active_schur)*(size_schur+1)]));
1531 ierr = MatSeqDenseInvertFactors_Private(M);CHKERRQ(ierr);
1532 /* move back just the active dofs to the Schur complement */
1533 for (i=0;i<size_active_schur;i++) {
1534 ierr = PetscArraycpy(data+i*size_schur,tdata+i*news,size_active_schur);CHKERRQ(ierr);
1535 }
1536 ierr = PetscFree(tdata);CHKERRQ(ierr);
1537 ierr = MatDestroy(&M);CHKERRQ(ierr);
1538 } else { /* we can factorize and invert just the activedofs */
1539 Mat M;
1540 PetscScalar *aux;
1541
1542 ierr = PetscMalloc1(nd,&aux);CHKERRQ(ierr);
1543 for (i=0;i<nd;i++) aux[i] = 1.0/data[(i+size_active_schur)*(size_schur+1)];
1544 ierr = MatCreateSeqDense(PETSC_COMM_SELF,size_active_schur,size_active_schur,data,&M);CHKERRQ(ierr);
1545 ierr = MatDenseSetLDA(M,size_schur);CHKERRQ(ierr);
1546 ierr = MatSetOption(M,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
1547 ierr = MatCholeskyFactor(M,NULL,NULL);CHKERRQ(ierr);
1548 ierr = MatSeqDenseInvertFactors_Private(M);CHKERRQ(ierr);
1549 ierr = MatDestroy(&M);CHKERRQ(ierr);
1550 ierr = MatCreateSeqDense(PETSC_COMM_SELF,size_schur,nd,data+size_active_schur*size_schur,&M);CHKERRQ(ierr);
1551 ierr = MatZeroEntries(M);CHKERRQ(ierr);
1552 ierr = MatDestroy(&M);CHKERRQ(ierr);
1553 ierr = MatCreateSeqDense(PETSC_COMM_SELF,nd,size_schur,data+size_active_schur,&M);CHKERRQ(ierr);
1554 ierr = MatDenseSetLDA(M,size_schur);CHKERRQ(ierr);
1555 ierr = MatZeroEntries(M);CHKERRQ(ierr);
1556 ierr = MatDestroy(&M);CHKERRQ(ierr);
1557 for (i=0;i<nd;i++) data[(i+size_active_schur)*(size_schur+1)] = aux[i];
1558 ierr = PetscFree(aux);CHKERRQ(ierr);
1559 }
1560 ierr = MatDenseRestoreArray(S_all_inv,&data);CHKERRQ(ierr);
1561 } else { /* use MatFactor calls to invert S */
1562 ierr = MatFactorInvertSchurComplement(F);CHKERRQ(ierr);
1563 ierr = MatFactorGetSchurComplement(F,&S_all_inv,NULL);CHKERRQ(ierr);
1564 }
1565 } else { /* we need to invert explicitly since we are not using MatFactor for S */
1566 ierr = PetscObjectReference((PetscObject)S_all);CHKERRQ(ierr);
1567 S_all_inv = S_all;
1568 ierr = MatDenseGetArray(S_all_inv,&S_data);CHKERRQ(ierr);
1569 ierr = PetscBLASIntCast(size_schur,&B_N);CHKERRQ(ierr);
1570 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
1571 if (use_potr) {
1572 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,S_data,&B_N,&B_ierr));
1573 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr);
1574 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,S_data,&B_N,&B_ierr));
1575 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr);
1576 } else if (use_sytr) {
1577 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&B_N,S_data,&B_N,pivots,Bwork,&B_lwork,&B_ierr));
1578 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYTRF Lapack routine %d",(int)B_ierr);
1579 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&B_N,S_data,&B_N,pivots,Bwork,&B_ierr));
1580 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYTRI Lapack routine %d",(int)B_ierr);
1581 } else {
1582 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,S_data,&B_N,pivots,&B_ierr));
1583 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
1584 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,S_data,&B_N,pivots,Bwork,&B_lwork,&B_ierr));
1585 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
1586 }
1587 ierr = PetscLogFlops(1.0*size_schur*size_schur*size_schur);CHKERRQ(ierr);
1588 ierr = PetscFPTrapPop();CHKERRQ(ierr);
1589 ierr = MatDenseRestoreArray(S_all_inv,&S_data);CHKERRQ(ierr);
1590 }
1591 /* S_Ej_tilda_all */
1592 cum = cum2 = 0;
1593 ierr = MatDenseGetArrayRead(S_all_inv,&rS_data);CHKERRQ(ierr);
1594 for (i=0;i<sub_schurs->n_subs;i++) {
1595 PetscInt j;
1596
1597 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
1598 /* get (St^-1)_E */
1599 /* Unless we are changing the variables, I don't need to expand to upper triangular since St^-1
1600 will be properly accessed later during adaptive selection */
1601 if (S_lower_triangular) {
1602 PetscInt k;
1603 if (sub_schurs->change) {
1604 for (k=0;k<subset_size;k++) {
1605 for (j=k;j<subset_size;j++) {
1606 work[k*subset_size+j] = rS_data[cum2+k*size_schur+j];
1607 work[j*subset_size+k] = work[k*subset_size+j];
1608 }
1609 }
1610 } else {
1611 for (k=0;k<subset_size;k++) {
1612 for (j=k;j<subset_size;j++) {
1613 work[k*subset_size+j] = rS_data[cum2+k*size_schur+j];
1614 }
1615 }
1616 }
1617 } else {
1618 PetscInt k;
1619 for (k=0;k<subset_size;k++) {
1620 for (j=0;j<subset_size;j++) {
1621 work[k*subset_size+j] = rS_data[cum2+k*size_schur+j];
1622 }
1623 }
1624 }
1625 if (sub_schurs->change) {
1626 Mat change_sub,SEj,T;
1627
1628 /* change basis */
1629 ierr = KSPGetOperators(sub_schurs->change[i],&change_sub,NULL);CHKERRQ(ierr);
1630 ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,work,&SEj);CHKERRQ(ierr);
1631 if (!sub_schurs->change_with_qr) { /* currently there's no support for PtAP with P SeqAIJ */
1632 Mat T2;
1633 ierr = MatTransposeMatMult(change_sub,SEj,MAT_INITIAL_MATRIX,1.0,&T2);CHKERRQ(ierr);
1634 ierr = MatMatMult(T2,change_sub,MAT_INITIAL_MATRIX,1.0,&T);CHKERRQ(ierr);
1635 ierr = MatDestroy(&T2);CHKERRQ(ierr);
1636 ierr = MatConvert(T,MATSEQDENSE,MAT_INPLACE_MATRIX,&T);CHKERRQ(ierr);
1637 } else {
1638 ierr = MatPtAP(SEj,change_sub,MAT_INITIAL_MATRIX,1.0,&T);CHKERRQ(ierr);
1639 }
1640 ierr = MatCopy(T,SEj,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
1641 ierr = MatDestroy(&T);CHKERRQ(ierr);
1642 /* set diagonal entry to a very large value to pick the basis we are eliminating as the first eigenvectors with adaptive selection */
1643 ierr = MatZeroRowsColumnsIS(SEj,sub_schurs->change_primal_sub[i],1./PETSC_SMALL,NULL,NULL);CHKERRQ(ierr);
1644 ierr = MatDestroy(&SEj);CHKERRQ(ierr);
1645 }
1646 ierr = PetscArraycpy(SEjinv_arr,work,subset_size*subset_size);CHKERRQ(ierr);
1647 cum += subset_size;
1648 cum2 += subset_size*(size_schur + 1);
1649 SEjinv_arr += subset_size*subset_size;
1650 }
1651 ierr = MatDenseRestoreArrayRead(S_all_inv,&rS_data);CHKERRQ(ierr);
1652 if (solver_S) {
1653 if (schur_has_vertices) {
1654 ierr = MatFactorRestoreSchurComplement(F,&S_all_inv,MAT_FACTOR_SCHUR_FACTORED);CHKERRQ(ierr);
1655 } else {
1656 ierr = MatFactorRestoreSchurComplement(F,&S_all_inv,MAT_FACTOR_SCHUR_INVERTED);CHKERRQ(ierr);
1657 }
1658 }
1659 ierr = MatDestroy(&S_all_inv);CHKERRQ(ierr);
1660 }
1661 ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all,&SEjinv_arr);CHKERRQ(ierr);
1662
1663 /* move back factors if needed */
1664 if (schur_has_vertices) {
1665 Mat S_tmp;
1666 PetscInt nd = 0;
1667
1668 if (!solver_S) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1669 ierr = MatFactorGetSchurComplement(F,&S_tmp,NULL);CHKERRQ(ierr);
1670 if (use_potr) {
1671 PetscScalar *data;
1672
1673 ierr = MatDenseGetArray(S_tmp,&data);CHKERRQ(ierr);
1674 ierr = PetscArrayzero(data,size_schur*size_schur);CHKERRQ(ierr);
1675
1676 if (S_lower_triangular) {
1677 cum = 0;
1678 for (i=0;i<size_active_schur;i++) {
1679 ierr = PetscArraycpy(data+i*(size_schur+1),schur_factor+cum,size_active_schur-i);CHKERRQ(ierr);
1680 cum += size_active_schur-i;
1681 }
1682 } else {
1683 ierr = PetscArraycpy(data,schur_factor,size_schur*size_schur);CHKERRQ(ierr);
1684 }
1685 if (sub_schurs->is_dir) {
1686 ierr = ISGetLocalSize(sub_schurs->is_dir,&nd);CHKERRQ(ierr);
1687 for (i=0;i<nd;i++) {
1688 data[(i+size_active_schur)*(size_schur+1)] = schur_factor[cum+i];
1689 }
1690 }
1691 /* workaround: since I cannot modify the matrices used inside the solvers for the forward and backward substitutions,
1692 set the diagonal entry of the Schur factor to a very large value */
1693 for (i=size_active_schur+nd;i<size_schur;i++) {
1694 data[i*(size_schur+1)] = infty;
1695 }
1696 ierr = MatDenseRestoreArray(S_tmp,&data);CHKERRQ(ierr);
1697 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor update not yet implemented for non SPD matrices");
1698 ierr = MatFactorRestoreSchurComplement(F,&S_tmp,MAT_FACTOR_SCHUR_FACTORED);CHKERRQ(ierr);
1699 }
1700 } else if (factor_workaround) { /* we need to eliminate any unneeded coupling */
1701 PetscScalar *data;
1702 PetscInt nd = 0;
1703
1704 if (sub_schurs->is_dir) { /* dirichlet dofs could have different scalings */
1705 ierr = ISGetLocalSize(sub_schurs->is_dir,&nd);CHKERRQ(ierr);
1706 }
1707 ierr = MatFactorGetSchurComplement(F,&S_all,NULL);CHKERRQ(ierr);
1708 ierr = MatDenseGetArray(S_all,&data);CHKERRQ(ierr);
1709 for (i=0;i<size_active_schur;i++) {
1710 ierr = PetscArrayzero(data+i*size_schur+size_active_schur,size_schur-size_active_schur);CHKERRQ(ierr);
1711 }
1712 for (i=size_active_schur+nd;i<size_schur;i++) {
1713 ierr = PetscArrayzero(data+i*size_schur+size_active_schur,size_schur-size_active_schur);CHKERRQ(ierr);
1714 data[i*(size_schur+1)] = infty;
1715 }
1716 ierr = MatDenseRestoreArray(S_all,&data);CHKERRQ(ierr);
1717 ierr = MatFactorRestoreSchurComplement(F,&S_all,MAT_FACTOR_SCHUR_UNFACTORED);CHKERRQ(ierr);
1718 }
1719 ierr = PetscFree(work);CHKERRQ(ierr);
1720 ierr = PetscFree(schur_factor);CHKERRQ(ierr);
1721 ierr = VecDestroy(&Dall);CHKERRQ(ierr);
1722 }
1723 ierr = ISDestroy(&is_I_layer);CHKERRQ(ierr);
1724 ierr = MatDestroy(&S_all);CHKERRQ(ierr);
1725 ierr = MatDestroy(&A_BB);CHKERRQ(ierr);
1726 ierr = MatDestroy(&A_IB);CHKERRQ(ierr);
1727 ierr = MatDestroy(&A_BI);CHKERRQ(ierr);
1728 ierr = MatDestroy(&F);CHKERRQ(ierr);
1729
1730 ierr = PetscMalloc1(sub_schurs->n_subs,&nnz);CHKERRQ(ierr);
1731 for (i=0;i<sub_schurs->n_subs;i++) {
1732 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&nnz[i]);CHKERRQ(ierr);
1733 }
1734 ierr = ISCreateGeneral(PETSC_COMM_SELF,sub_schurs->n_subs,nnz,PETSC_OWN_POINTER,&is_I_layer);CHKERRQ(ierr);
1735 ierr = MatSetVariableBlockSizes(sub_schurs->S_Ej_all,sub_schurs->n_subs,nnz);CHKERRQ(ierr);
1736 ierr = MatAssemblyBegin(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1737 ierr = MatAssemblyEnd(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1738 if (compute_Stilda) {
1739 ierr = MatSetVariableBlockSizes(sub_schurs->sum_S_Ej_tilda_all,sub_schurs->n_subs,nnz);CHKERRQ(ierr);
1740 ierr = MatAssemblyBegin(sub_schurs->sum_S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1741 ierr = MatAssemblyEnd(sub_schurs->sum_S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1742 if (deluxe) {
1743 ierr = MatSetVariableBlockSizes(sub_schurs->sum_S_Ej_inv_all,sub_schurs->n_subs,nnz);CHKERRQ(ierr);
1744 ierr = MatAssemblyBegin(sub_schurs->sum_S_Ej_inv_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1745 ierr = MatAssemblyEnd(sub_schurs->sum_S_Ej_inv_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1746 }
1747 }
1748 ierr = ISDestroy(&is_I_layer);CHKERRQ(ierr);
1749
1750 /* Get local part of (\sum_j S_Ej) */
1751 if (!sub_schurs->sum_S_Ej_all) {
1752 ierr = MatDuplicate(sub_schurs->S_Ej_all,MAT_DO_NOT_COPY_VALUES,&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
1753 }
1754 ierr = VecSet(gstash,0.0);CHKERRQ(ierr);
1755 ierr = MatSeqAIJGetArray(sub_schurs->S_Ej_all,&stasharray);CHKERRQ(ierr);
1756 ierr = VecPlaceArray(lstash,stasharray);CHKERRQ(ierr);
1757 ierr = VecScatterBegin(sstash,lstash,gstash,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1758 ierr = VecScatterEnd(sstash,lstash,gstash,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1759 ierr = MatSeqAIJRestoreArray(sub_schurs->S_Ej_all,&stasharray);CHKERRQ(ierr);
1760 ierr = VecResetArray(lstash);CHKERRQ(ierr);
1761 ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_all,&stasharray);CHKERRQ(ierr);
1762 ierr = VecPlaceArray(lstash,stasharray);CHKERRQ(ierr);
1763 ierr = VecScatterBegin(sstash,gstash,lstash,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1764 ierr = VecScatterEnd(sstash,gstash,lstash,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1765 ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_all,&stasharray);CHKERRQ(ierr);
1766 ierr = VecResetArray(lstash);CHKERRQ(ierr);
1767
1768 /* Get local part of (\sum_j S^-1_Ej) (\sum_j St^-1_Ej) */
1769 if (compute_Stilda) {
1770 ierr = VecSet(gstash,0.0);CHKERRQ(ierr);
1771 ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all,&stasharray);CHKERRQ(ierr);
1772 ierr = VecPlaceArray(lstash,stasharray);CHKERRQ(ierr);
1773 ierr = VecScatterBegin(sstash,lstash,gstash,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1774 ierr = VecScatterEnd(sstash,lstash,gstash,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1775 ierr = VecScatterBegin(sstash,gstash,lstash,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1776 ierr = VecScatterEnd(sstash,gstash,lstash,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1777 ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all,&stasharray);CHKERRQ(ierr);
1778 ierr = VecResetArray(lstash);CHKERRQ(ierr);
1779 if (deluxe) {
1780 ierr = VecSet(gstash,0.0);CHKERRQ(ierr);
1781 ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_inv_all,&stasharray);CHKERRQ(ierr);
1782 ierr = VecPlaceArray(lstash,stasharray);CHKERRQ(ierr);
1783 ierr = VecScatterBegin(sstash,lstash,gstash,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1784 ierr = VecScatterEnd(sstash,lstash,gstash,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1785 ierr = VecScatterBegin(sstash,gstash,lstash,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1786 ierr = VecScatterEnd(sstash,gstash,lstash,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1787 ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_inv_all,&stasharray);CHKERRQ(ierr);
1788 ierr = VecResetArray(lstash);CHKERRQ(ierr);
1789 } else {
1790 PetscScalar *array;
1791 PetscInt cum;
1792
1793 ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all,&array);CHKERRQ(ierr);
1794 cum = 0;
1795 for (i=0;i<sub_schurs->n_subs;i++) {
1796 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
1797 ierr = PetscBLASIntCast(subset_size,&B_N);CHKERRQ(ierr);
1798 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
1799 if (use_potr) {
1800 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,array+cum,&B_N,&B_ierr));
1801 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr);
1802 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,array+cum,&B_N,&B_ierr));
1803 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr);
1804 } else if (use_sytr) {
1805 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&B_N,array+cum,&B_N,pivots,Bwork,&B_lwork,&B_ierr));
1806 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYTRF Lapack routine %d",(int)B_ierr);
1807 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&B_N,array+cum,&B_N,pivots,Bwork,&B_ierr));
1808 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYTRI Lapack routine %d",(int)B_ierr);
1809 } else {
1810 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,array+cum,&B_N,pivots,&B_ierr));
1811 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr);
1812 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,array+cum,&B_N,pivots,Bwork,&B_lwork,&B_ierr));
1813 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr);
1814 }
1815 ierr = PetscLogFlops(1.0*subset_size*subset_size*subset_size);CHKERRQ(ierr);
1816 ierr = PetscFPTrapPop();CHKERRQ(ierr);
1817 cum += subset_size*subset_size;
1818 }
1819 ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all,&array);CHKERRQ(ierr);
1820 ierr = PetscObjectReference((PetscObject)sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
1821 ierr = MatDestroy(&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
1822 sub_schurs->sum_S_Ej_inv_all = sub_schurs->sum_S_Ej_all;
1823 }
1824 }
1825 ierr = VecDestroy(&lstash);CHKERRQ(ierr);
1826 ierr = VecDestroy(&gstash);CHKERRQ(ierr);
1827 ierr = VecScatterDestroy(&sstash);CHKERRQ(ierr);
1828
1829 if (matl_dbg_viewer) {
1830 PetscInt cum;
1831
1832 if (sub_schurs->S_Ej_all) {
1833 ierr = PetscObjectSetName((PetscObject)sub_schurs->S_Ej_all,"SE");CHKERRQ(ierr);
1834 ierr = MatView(sub_schurs->S_Ej_all,matl_dbg_viewer);CHKERRQ(ierr);
1835 }
1836 if (sub_schurs->sum_S_Ej_all) {
1837 ierr = PetscObjectSetName((PetscObject)sub_schurs->sum_S_Ej_all,"SSE");CHKERRQ(ierr);
1838 ierr = MatView(sub_schurs->sum_S_Ej_all,matl_dbg_viewer);CHKERRQ(ierr);
1839 }
1840 if (sub_schurs->sum_S_Ej_inv_all) {
1841 ierr = PetscObjectSetName((PetscObject)sub_schurs->sum_S_Ej_inv_all,"SSEm");CHKERRQ(ierr);
1842 ierr = MatView(sub_schurs->sum_S_Ej_inv_all,matl_dbg_viewer);CHKERRQ(ierr);
1843 }
1844 if (sub_schurs->sum_S_Ej_tilda_all) {
1845 ierr = PetscObjectSetName((PetscObject)sub_schurs->sum_S_Ej_tilda_all,"SSEt");CHKERRQ(ierr);
1846 ierr = MatView(sub_schurs->sum_S_Ej_tilda_all,matl_dbg_viewer);CHKERRQ(ierr);
1847 }
1848 for (i=0,cum=0;i<sub_schurs->n_subs;i++) {
1849 IS is;
1850 char name[16];
1851
1852 ierr = PetscSNPrintf(name,sizeof(name),"IE%D",i);CHKERRQ(ierr);
1853 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr);
1854 ierr = ISCreateStride(PETSC_COMM_SELF,subset_size,cum,1,&is);CHKERRQ(ierr);
1855 ierr = PetscObjectSetName((PetscObject)is,name);CHKERRQ(ierr);
1856 ierr = ISView(is,matl_dbg_viewer);CHKERRQ(ierr);
1857 ierr = ISDestroy(&is);CHKERRQ(ierr);
1858 cum += subset_size;
1859 }
1860 }
1861
1862 /* free workspace */
1863 ierr = PetscViewerDestroy(&matl_dbg_viewer);CHKERRQ(ierr);
1864 ierr = PetscFree2(Bwork,pivots);CHKERRQ(ierr);
1865 ierr = PetscCommDestroy(&comm_n);CHKERRQ(ierr);
1866 PetscFunctionReturn(0);
1867 }
1868
PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs,const char * prefix,IS is_I,IS is_B,PCBDDCGraph graph,ISLocalToGlobalMapping BtoNmap,PetscBool copycc)1869 PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, const char* prefix, IS is_I, IS is_B, PCBDDCGraph graph, ISLocalToGlobalMapping BtoNmap, PetscBool copycc)
1870 {
1871 IS *faces,*edges,*all_cc,vertices;
1872 PetscInt i,n_faces,n_edges,n_all_cc;
1873 PetscBool is_sorted,ispardiso,ismumps;
1874 PetscErrorCode ierr;
1875
1876 PetscFunctionBegin;
1877 ierr = ISSorted(is_I,&is_sorted);CHKERRQ(ierr);
1878 if (!is_sorted) SETERRQ(PetscObjectComm((PetscObject)is_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted");
1879 ierr = ISSorted(is_B,&is_sorted);CHKERRQ(ierr);
1880 if (!is_sorted) SETERRQ(PetscObjectComm((PetscObject)is_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted");
1881
1882 /* reset any previous data */
1883 ierr = PCBDDCSubSchursReset(sub_schurs);CHKERRQ(ierr);
1884
1885 /* get index sets for faces and edges (already sorted by global ordering) */
1886 ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,&vertices);CHKERRQ(ierr);
1887 n_all_cc = n_faces+n_edges;
1888 ierr = PetscBTCreate(n_all_cc,&sub_schurs->is_edge);CHKERRQ(ierr);
1889 ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr);
1890 for (i=0;i<n_faces;i++) {
1891 if (copycc) {
1892 ierr = ISDuplicate(faces[i],&all_cc[i]);CHKERRQ(ierr);
1893 } else {
1894 ierr = PetscObjectReference((PetscObject)faces[i]);CHKERRQ(ierr);
1895 all_cc[i] = faces[i];
1896 }
1897 }
1898 for (i=0;i<n_edges;i++) {
1899 if (copycc) {
1900 ierr = ISDuplicate(edges[i],&all_cc[n_faces+i]);CHKERRQ(ierr);
1901 } else {
1902 ierr = PetscObjectReference((PetscObject)edges[i]);CHKERRQ(ierr);
1903 all_cc[n_faces+i] = edges[i];
1904 }
1905 ierr = PetscBTSet(sub_schurs->is_edge,n_faces+i);CHKERRQ(ierr);
1906 }
1907 ierr = PetscObjectReference((PetscObject)vertices);CHKERRQ(ierr);
1908 sub_schurs->is_vertices = vertices;
1909 ierr = PCBDDCGraphRestoreCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,&vertices);CHKERRQ(ierr);
1910 sub_schurs->is_dir = NULL;
1911 ierr = PCBDDCGraphGetDirichletDofsB(graph,&sub_schurs->is_dir);CHKERRQ(ierr);
1912
1913 /* Determine if MatFactor can be used */
1914 ierr = PetscStrallocpy(prefix,&sub_schurs->prefix);CHKERRQ(ierr);
1915 #if defined(PETSC_HAVE_MUMPS)
1916 ierr = PetscStrncpy(sub_schurs->mat_solver_type,MATSOLVERMUMPS,sizeof(sub_schurs->mat_solver_type));CHKERRQ(ierr);
1917 #elif defined(PETSC_HAVE_MKL_PARDISO)
1918 ierr = PetscStrncpy(sub_schurs->mat_solver_type,MATSOLVERMKL_PARDISO,sizeof(sub_schurs->mat_solver_type));CHKERRQ(ierr);
1919 #else
1920 ierr = PetscStrncpy(sub_schurs->mat_solver_type,MATSOLVERPETSC,sizeof(sub_schurs->mat_solver_type));CHKERRQ(ierr);
1921 #endif
1922 #if defined(PETSC_USE_COMPLEX)
1923 sub_schurs->is_hermitian = PETSC_FALSE; /* Hermitian Cholesky is not supported by PETSc and external packages */
1924 #else
1925 sub_schurs->is_hermitian = PETSC_TRUE;
1926 #endif
1927 sub_schurs->is_posdef = PETSC_TRUE;
1928 sub_schurs->is_symmetric = PETSC_TRUE;
1929 sub_schurs->debug = PETSC_FALSE;
1930 sub_schurs->restrict_comm = PETSC_FALSE;
1931 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)graph->l2gmap),sub_schurs->prefix,"BDDC sub_schurs options","PC");CHKERRQ(ierr);
1932 ierr = PetscOptionsString("-sub_schurs_mat_solver_type","Specific direct solver to use",NULL,sub_schurs->mat_solver_type,sub_schurs->mat_solver_type,sizeof(sub_schurs->mat_solver_type),NULL);CHKERRQ(ierr);
1933 ierr = PetscOptionsBool("-sub_schurs_symmetric","Symmetric problem",NULL,sub_schurs->is_symmetric,&sub_schurs->is_symmetric,NULL);CHKERRQ(ierr);
1934 ierr = PetscOptionsBool("-sub_schurs_hermitian","Hermitian problem",NULL,sub_schurs->is_hermitian,&sub_schurs->is_hermitian,NULL);CHKERRQ(ierr);
1935 ierr = PetscOptionsBool("-sub_schurs_posdef","Positive definite problem",NULL,sub_schurs->is_posdef,&sub_schurs->is_posdef,NULL);CHKERRQ(ierr);
1936 ierr = PetscOptionsBool("-sub_schurs_restrictcomm","Restrict communicator on active processes only",NULL,sub_schurs->restrict_comm,&sub_schurs->restrict_comm,NULL);CHKERRQ(ierr);
1937 ierr = PetscOptionsBool("-sub_schurs_debug","Debug output",NULL,sub_schurs->debug,&sub_schurs->debug,NULL);CHKERRQ(ierr);
1938 ierr = PetscOptionsEnd();CHKERRQ(ierr);
1939 ierr = PetscStrcmp(sub_schurs->mat_solver_type,MATSOLVERMUMPS,&ismumps);CHKERRQ(ierr);
1940 ierr = PetscStrcmp(sub_schurs->mat_solver_type,MATSOLVERMKL_PARDISO,&ispardiso);CHKERRQ(ierr);
1941 sub_schurs->schur_explicit = (PetscBool)(ispardiso || ismumps);
1942
1943 /* for reals, symmetric and hermitian are synonims */
1944 #if !defined(PETSC_USE_COMPLEX)
1945 sub_schurs->is_symmetric = (PetscBool)(sub_schurs->is_symmetric && sub_schurs->is_hermitian);
1946 sub_schurs->is_hermitian = sub_schurs->is_symmetric;
1947 #endif
1948
1949 ierr = PetscObjectReference((PetscObject)is_I);CHKERRQ(ierr);
1950 sub_schurs->is_I = is_I;
1951 ierr = PetscObjectReference((PetscObject)is_B);CHKERRQ(ierr);
1952 sub_schurs->is_B = is_B;
1953 ierr = PetscObjectReference((PetscObject)graph->l2gmap);CHKERRQ(ierr);
1954 sub_schurs->l2gmap = graph->l2gmap;
1955 ierr = PetscObjectReference((PetscObject)BtoNmap);CHKERRQ(ierr);
1956 sub_schurs->BtoNmap = BtoNmap;
1957 sub_schurs->n_subs = n_all_cc;
1958 sub_schurs->is_subs = all_cc;
1959 sub_schurs->S_Ej_all = NULL;
1960 sub_schurs->sum_S_Ej_all = NULL;
1961 sub_schurs->sum_S_Ej_inv_all = NULL;
1962 sub_schurs->sum_S_Ej_tilda_all = NULL;
1963 sub_schurs->is_Ej_all = NULL;
1964 PetscFunctionReturn(0);
1965 }
1966
PCBDDCSubSchursCreate(PCBDDCSubSchurs * sub_schurs)1967 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs)
1968 {
1969 PCBDDCSubSchurs schurs_ctx;
1970 PetscErrorCode ierr;
1971
1972 PetscFunctionBegin;
1973 ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr);
1974 schurs_ctx->n_subs = 0;
1975 *sub_schurs = schurs_ctx;
1976 PetscFunctionReturn(0);
1977 }
1978
PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs)1979 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs)
1980 {
1981 PetscInt i;
1982 PetscErrorCode ierr;
1983
1984 PetscFunctionBegin;
1985 if (!sub_schurs) PetscFunctionReturn(0);
1986 ierr = PetscFree(sub_schurs->prefix);CHKERRQ(ierr);
1987 ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr);
1988 ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr);
1989 ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr);
1990 ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr);
1991 ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->l2gmap);CHKERRQ(ierr);
1992 ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->BtoNmap);CHKERRQ(ierr);
1993 ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr);
1994 ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
1995 ierr = MatDestroy(&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr);
1996 ierr = MatDestroy(&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr);
1997 ierr = ISDestroy(&sub_schurs->is_Ej_all);CHKERRQ(ierr);
1998 ierr = ISDestroy(&sub_schurs->is_vertices);CHKERRQ(ierr);
1999 ierr = ISDestroy(&sub_schurs->is_dir);CHKERRQ(ierr);
2000 ierr = PetscBTDestroy(&sub_schurs->is_edge);CHKERRQ(ierr);
2001 for (i=0;i<sub_schurs->n_subs;i++) {
2002 ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr);
2003 }
2004 if (sub_schurs->n_subs) {
2005 ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr);
2006 }
2007 if (sub_schurs->reuse_solver) {
2008 ierr = PCBDDCReuseSolversReset(sub_schurs->reuse_solver);CHKERRQ(ierr);
2009 }
2010 ierr = PetscFree(sub_schurs->reuse_solver);CHKERRQ(ierr);
2011 if (sub_schurs->change) {
2012 for (i=0;i<sub_schurs->n_subs;i++) {
2013 ierr = KSPDestroy(&sub_schurs->change[i]);CHKERRQ(ierr);
2014 ierr = ISDestroy(&sub_schurs->change_primal_sub[i]);CHKERRQ(ierr);
2015 }
2016 }
2017 ierr = PetscFree(sub_schurs->change);CHKERRQ(ierr);
2018 ierr = PetscFree(sub_schurs->change_primal_sub);CHKERRQ(ierr);
2019 sub_schurs->n_subs = 0;
2020 PetscFunctionReturn(0);
2021 }
2022
PCBDDCSubSchursDestroy(PCBDDCSubSchurs * sub_schurs)2023 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs* sub_schurs)
2024 {
2025 PetscErrorCode ierr;
2026
2027 PetscFunctionBegin;
2028 ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr);
2029 ierr = PetscFree(*sub_schurs);CHKERRQ(ierr);
2030 PetscFunctionReturn(0);
2031 }
2032
PCBDDCAdjGetNextLayer_Private(PetscInt * queue_tip,PetscInt n_prev,PetscBT touched,PetscInt * xadj,PetscInt * adjncy,PetscInt * n_added)2033 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added)
2034 {
2035 PetscInt i,j,n;
2036 PetscErrorCode ierr;
2037
2038 PetscFunctionBegin;
2039 n = 0;
2040 for (i=-n_prev;i<0;i++) {
2041 PetscInt start_dof = queue_tip[i];
2042 for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
2043 PetscInt dof = adjncy[j];
2044 if (!PetscBTLookup(touched,dof)) {
2045 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
2046 queue_tip[n] = dof;
2047 n++;
2048 }
2049 }
2050 }
2051 *n_added = n;
2052 PetscFunctionReturn(0);
2053 }
2054