PetscSolverFeti.cc 77.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscHelper.h"
16
#include "parallel/PetscSolverFeti.h"
17
18
#include "parallel/PetscSolverFetiDebug.h"
#include "parallel/PetscSolverFetiMonitor.h"
19
#include "parallel/PetscSolverFetiStructs.h"
20
21
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
22
23
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
24
#include "parallel/PetscSolverGlobalMatrix.h"
25
#include "io/VtkWriter.h"
26
27
28
29
30

namespace AMDiS {

  using namespace std;

31
32
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
33
34
35
36
37
38
      primalDofMap(COMPONENT_WISE),
      dualDofMap(COMPONENT_WISE),
      interfaceDofMap(COMPONENT_WISE),
      localDofMap(COMPONENT_WISE),
      lagrangeMap(COMPONENT_WISE),
      interiorDofMap(COMPONENT_WISE),
39
      schurPrimalSolver(0),
40
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
41
      subdomain(NULL),
42
      massMatrixSolver(NULL),
43
44
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
45
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
46
      printTimings(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
47
      stokesMode(false),
48
      augmentedLagrange(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
49
      pressureComponent(-1)
50
51
52
53
54
55
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
56
      MSG("Create FETI-DP solver with no preconditioner!\n");
57
58
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
59
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
60
61
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
62
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
63
64
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
65
66
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
67
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
68
69
70
71
72

    Parameters::get("parallel->feti->schur primal solver", schurPrimalSolver);
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
      ("Wrong solver \"%d\"for the Schur primal complement!\n", 
       schurPrimalSolver);
73
74

    Parameters::get("parallel->multi level test", multiLevelTest);
75
76
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
77
78

    Parameters::get("parallel->print timings", printTimings);
79
80

    Parameters::get("parallel->feti->augmented lagrange", augmentedLagrange);
81
82

    Parameters::get("parallel->feti->symmetric", isSymmetric);
83
84
85
  }


86
  void PetscSolverFeti::initialize()
87
  {
88
89
    FUNCNAME("PetscSolverFeti::initialize()");

90
91
92
    TEST_EXIT_DBG(meshLevel + 1 == meshDistributor->getMeshLevelData().getLevelNumber())
      ("Mesh hierarchy does not contain %d levels!\n", meshLevel + 1);

93
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
Thomas Witkowski's avatar
Thomas Witkowski committed
94

Thomas Witkowski's avatar
Thomas Witkowski committed
95
96
97
98
99
100
    stokesMode = false;
    Parameters::get("parallel->feti->stokes mode", stokesMode);
    if (stokesMode) {
      Parameters::get("parallel->feti->pressure component", pressureComponent);
      TEST_EXIT(pressureComponent >= 0)
	("FETI-DP in Stokes mode, no pressure component defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
103
104
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
105
      subdomain->setSymmetric(isSymmetric);
106
      subdomain->setHandleDirichletRows(false);
107

108
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
109
	subdomain->setMeshDistributor(meshDistributor, 
110
				      mpiCommGlobal, mpiCommLocal);
111
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
112
	subdomain->setMeshDistributor(meshDistributor, 
113
114
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
115
	subdomain->setLevel(meshLevel);
116
117
      }
    }
118

119
120
121
122
    primalDofMap.init(levelData, componentSpaces, feSpaces);   
    dualDofMap.init(levelData, componentSpaces, feSpaces, false);
    localDofMap.init(levelData, componentSpaces, feSpaces, meshLevel != 0);
    lagrangeMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
123

Thomas Witkowski's avatar
Thomas Witkowski committed
124
    if (stokesMode)
125
      interfaceDofMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
126

127
    if (fetiPreconditioner == FETI_DIRICHLET) {
128
129
130
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

131
      interiorDofMap.init(levelData, componentSpaces, feSpaces, false);
132
    }
133
134
135
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
136
137
138
139
140
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

    int nComponents = mat.getSize();
141
142
    for (int component = 0; component < nComponents; component++) {
      DOFMatrix* dofMat = mat[component][component];
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
143
144
      if (!dofMat)
	continue;
145
146

      dirichletRows[component] = dofMat->getDirichletRows();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
147
148
149
150
    }
  }


151
152
153
154
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
155
156
    double timeCounter = MPI::Wtime();

157
158
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

159
160
161
162
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
163
    if (fetiPreconditioner == FETI_DIRICHLET)
164
165
      interiorDofMap.clear();

166
167
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
168

169
170
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
171
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
172
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
173
    if (fetiPreconditioner == FETI_DIRICHLET)
174
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
175

176
177
178
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
179
180
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
181
    if (stokesMode) {
182
183
184
185
186
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

187
188
189
190
191
192
    int nComponents = componentSpaces.size();
    for (int component = 0; component < nComponents; component++) {
      createPrimals(component);  
      createDuals(component);
      createInterfaceNodes(component);
      createIndexB(component);
193
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
196
197

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
198
    if (fetiPreconditioner == FETI_DIRICHLET)
199
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
200

Thomas Witkowski's avatar
Thomas Witkowski committed
201
    if (stokesMode)
202
203
      interfaceDofMap.update();

204
205
206
    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
207
    }
208

209
210
211
    lagrangeMap.update();


212
213
214
215
216
217
218
219
220
221
222
223
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

224
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
225
226
227
228
229
230
231
232
233
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	tmp = rStartInterior;

      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }

234
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
235
236
    for (unsigned int i = 0; i < componentSpaces.size(); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];
237
238

      MSG("FETI-DP data for %d-ith component (FE space %p):\n", i, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
239

240
      if (i == pressureComponent) {
Thomas Witkowski's avatar
Thomas Witkowski committed
241
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
242
243
	    interfaceDofMap[i].nRankDofs, 
	    interfaceDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
244
      } else {
245
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
246
247
248
	    primalDofMap[i].nRankDofs, 
	    primalDofMap[i].nLocalDofs,
	    primalDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
251
252
	    dualDofMap[i].nRankDofs, 
	    dualDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
253
254
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
255
256
	    lagrangeMap[i].nRankDofs, 
	    lagrangeMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
257
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
258
    }
259

Thomas Witkowski's avatar
Thomas Witkowski committed
260
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
263
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
266
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
267
      timeCounter = MPI::Wtime() - timeCounter;
268
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
	  timeCounter);
    }
271
272
273
  }


274
  void PetscSolverFeti::createPrimals(int component)
Thomas Witkowski's avatar
Thomas Witkowski committed
275
  {
276
    FUNCNAME("PetscSolverFeti::createPrimals()");  
277

278
    if (component == pressureComponent)
279
280
      return;

281
282
    const FiniteElemSpace *feSpace = componentSpaces[component];

283
284
285
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
286
    // Set of DOF indices that are considered to be primal variables.
287
    DofContainerSet& vertices = 
288
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
289
290

    DofIndexSet primals;
291
    for (DofContainerSet::iterator it = vertices.begin(); 
292
	 it != vertices.end(); ++it) {
293
      
294
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
295
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
296

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
297
298
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
299
      } else {
300
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
301
302
303
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
304
305
306
307
308
309
310
	if (fabs(c[0]) < e || fabs(c[1]) < e ||
	    fabs(c[0] - 25.0) < e || fabs(c[1] - 25.0) < e ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 12.5) < e)) {
	  MSG("PRIMAL COORD %f %f\n", c[0], c[1]);
	  primals.insert(**it);
	}
      }
311
    }
312

313
314
315
316

    // === Calculate the number of primals that are owned by the rank and ===
    // === create local indices of the primals starting at zero.          ===

317
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
318
      if (dofMap[feSpace].isRankDof(*it)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
319
	primalDofMap[component].insertRankDof(*it);
320
      } else
Thomas Witkowski's avatar
Thomas Witkowski committed
321
  	primalDofMap[component].insertNonRankDof(*it);
322
323
324
  }


325
  void PetscSolverFeti::createDuals(int component)
326
327
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
328

329
    if (component == pressureComponent)
330
331
      return;

332
333
    const FiniteElemSpace *feSpace = componentSpaces[component];

334
335
336
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
337
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
338
    
339
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
340
	 it != allBoundaryDofs.end(); ++it) {
341
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
342
343
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
344
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
345
346
347
	continue;

      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
348
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
349
      } else {
350
	if (dofMapSd[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	  dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
352
353
      }	  
    }
354
355
356
  }

  
357
  void PetscSolverFeti::createInterfaceNodes(int component)
358
359
360
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

361
    if (component != pressureComponent)
362
363
      return;

364
    const FiniteElemSpace *feSpace = componentSpaces[component];
365
366
367
368
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
369
	 it != allBoundaryDofs.end(); ++it) {
370
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
371
372
	continue;      
      
373
      if (dofMap[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
374
	interfaceDofMap[component].insertRankDof(**it);
375
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
376
	interfaceDofMap[component].insertNonRankDof(**it);      
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
377
    }
378
379
380
  }


381
  void PetscSolverFeti::createLagrange(int component)
382
383
384
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

385
    if (component == pressureComponent)
386
387
      return;

388
    const FiniteElemSpace *feSpace = componentSpaces[component];
389
390
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
393
    // Stores for all rank owned communication DOFs, if the counterpart is
    // a rank owned DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFs that fulfill this requirenment.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
394
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
395

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
396
    if (meshLevel > 0) {
397
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
398
399
400
401
402
403
404
405
406

      for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank()) {

	vector<int> subdomainRankDofs;
	subdomainRankDofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
407
	  if (dofMapSd[feSpace].isRankDof(it.getDofIndex()))
408
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
409
	  else
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
	    subdomainRankDofs.push_back(0);
	}

	stdMpi.send(it.getRank(), subdomainRankDofs);
      }	     

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace); 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
425
426
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
427
	  if (!isPrimal(component, it.getDofIndex()))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
428
429
430
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
431

Thomas Witkowski's avatar
Thomas Witkowski committed
432
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
433
434
435
436
437
438
      return;


    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===

439
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
440
441
442
443
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
444
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
445
446
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

447
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
448
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
449
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
450
451
452
453
	}
      }
    }

454
455
456
457

    // === Communicate these sets for all rank owned dual nodes to other ===
    // === ranks that also have this node.                               ===

458
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
459

460
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
461
462
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
463
	if (!isPrimal(component, it.getDofIndex()))
464
465
466
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
467
468
469

    stdMpi.updateSendDataSize();

470
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
471
	 !it.end(); it.nextRank()) {
472
      bool recvFromRank = false;
473
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
474
	if (!isPrimal(component, it.getDofIndex())) {
475
476
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
477
	       dofMapSd[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
478
479
480
	    recvFromRank = true;
	    break;
	  }
481
	}
482
      }
483
484

      if (recvFromRank)
485
	stdMpi.recv(it.getRank());
486
    }
487

488
489
    stdMpi.startCommunication();

490
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
491
	 !it.end(); it.nextRank()) {
492
      int i = 0;
493
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
494
	if (!isPrimal(component, it.getDofIndex()))
495
496
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
497
	       dofMapSd[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
498
499
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
500
	  else
Thomas Witkowski's avatar
Thomas Witkowski committed
501
	    lagrangeMap[component].insertNonRankDof(it.getDofIndex());
502
503
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
504

505
506
507
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

508
    int nRankLagrange = 0;
509
    DofMap& dualMap = dualDofMap[component].getMap();
510
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
511
      if (dofMap[feSpace].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
512
	lagrangeMap[component].insertRankDof(it->first, nRankLagrange);
513
	int degree = boundaryDofRanks[feSpace][it->first].size();
514
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
515
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
516
	lagrangeMap[component].insertNonRankDof(it->first);
517
518
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
519
    lagrangeMap[component].nRankDofs = nRankLagrange;
520
521
522
  }


523
  void PetscSolverFeti::createAugmentedLagrange(int component)
524
525
526
527
528
529
530
531
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


532
  void PetscSolverFeti::createIndexB(int component)
533
  {
534
    FUNCNAME("PetscSolverFeti::createIndexB()");
535

536
537

    const FiniteElemSpace *feSpace = componentSpaces[component];
538
    DOFAdmin* admin = feSpace->getAdmin();
539
540
541
542

    // === To ensure that all interior node on each rank are listen first in ===
    // === the global index of all B nodes, insert all interior nodes first, ===
    // === without defining a correct index.                                 ===
543

544
    int nLocalInterior = 0;    
545
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
546
      if (admin->isDofFree(i) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
547
548
549
	  isPrimal(component, i) ||
	  isDual(component, i) ||
	  isInterface(component, i) ||
550
	  dirichletRows[component].count(i))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
551
	continue;      
552

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
553
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
554
	localDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
555
556
	
	if (fetiPreconditioner == FETI_DIRICHLET)
557
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
558
559
560
	
	nLocalInterior++;	
      } else {
561
	if (dofMapSd[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
562
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
563
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
564
	  localDofMap[component].insertNonRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
565
566
567
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
568
      }
569
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
570
    
571
572
    // === And finally, add the global indicies of all dual nodes. ===

573
574
    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
	 it != dualDofMap[component].getMap().end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
575
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
576
	localDofMap[component].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
577
      } else {
578
	if (dofMapSd[feSpace].isRankDof(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
579
	  localDofMap[component].insertRankDof(it->first);
580
	else 
Thomas Witkowski's avatar
Thomas Witkowski committed
581
	  localDofMap[component].insertNonRankDof(it->first);
582
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
583
    }
584
585
586
  }


587
  void PetscSolverFeti::createMatLagrange()
588
589
590
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
591
    double wtime = MPI::Wtime();
592
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
593

594
595
    // === Create distributed matrix for Lagrange constraints. ===

596
597
598
599
600
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
601
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
602

603
604
605
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

606
607
608
609
610
611
612
    // === Create for all duals the corresponding Lagrange constraints. On ===
    // === each rank we traverse all pairs (n, m) of ranks, with n < m,    ===
    // === that contain this node. If the current rank number is r, and    ===
    // === n == r, the rank sets 1.0 for the corresponding constraint, if  ===
    // === m == r, than the rank sets -1.0 for the corresponding           ===
    // === constraint.                                                     ===

613
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
614
      DofMap &dualMap = dualDofMap[component].getMap();
615

616
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
617
	TEST_EXIT_DBG(boundaryDofRanks[componentSpaces[component]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
620
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
621
	int index = lagrangeMap.getMatIndex(component, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
622

Thomas Witkowski's avatar
Thomas Witkowski committed
623
	// Copy set of all ranks that contain this dual node.
624
625
	vector<int> W(boundaryDofRanks[componentSpaces[component]][it->first].begin(), 
		      boundaryDofRanks[componentSpaces[component]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
626
627
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
628
629

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
630
	
631
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
634
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
635
636
	      MatSetValue(mat_lagrange, 
			  index + counter, 
637
			  localDofMap.getMatIndex(component, it->first) + rStartInterior,
638
639
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
640
	    }
641
642
643
644
645
646
647
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
648
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
649
650
651
652
653
654
	  int nConstraints = (degree * (degree - 1)) / 2;
	  for (int i = 0; i < nConstraints; i++) {
	    VecSetValue(vec_scale_lagrange,
			index + i,
			1.0 / static_cast<double>(degree),
			INSERT_VALUES);
655
656
657
658
659
660
661
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
662

663

Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
664
665
666
667
668
669
670
671
672
673
674
675
    int nZeroRows = PetscSolverFetiDebug::testZeroRows(mat_lagrange);
    int m,n;
    MatGetSize(mat_lagrange, &m ,&n);
    MSG("Lagrange matrix has %d zero rows and global size of %d %d!\n", nZeroRows, m, n);

    PetscViewer petscView;
    PetscViewerBinaryOpen(PETSC_COMM_WORLD, "lagrange.mat", 
			  FILE_MODE_WRITE, &petscView);
    MatView(mat_lagrange, petscView);
    PetscViewerDestroy(&petscView);


676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

    if (fetiPreconditioner != FETI_NONE || stokesMode) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatDiagonalScale(mat_lagrange_scaled, vec_scale_lagrange, PETSC_NULL);
    }

    VecDestroy(&vec_scale_lagrange);


    // === Print final timings. ===

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
691
692
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
693
694
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
695
    }
696
697
  }

698

Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
  {
    Element *el = edge.el;
    int i0 = el->getVertexOfEdge(edge.ithObj, 0);
    int i1 = el->getVertexOfEdge(edge.ithObj, 1);
    DegreeOfFreedom d0 = el->getDof(i0, 0);
    DegreeOfFreedom d1 = el->getDof(i1, 0);
    WorldVector<double> c0, c1;
    el->getMesh()->getDofIndexCoords(d0, feSpace, c0);
    el->getMesh()->getDofIndexCoords(d1, feSpace, c1);
    bool xe = fabs(c0[0] - c1[0]) < 1e-8;
    bool ye = fabs(c0[1] - c1[1]) < 1e-8;
    bool ze = fabs(c0[2] - c1[2]) < 1e-8;
    int counter = static_cast<int>(xe) + static_cast<int>(ye) + static_cast<int>(ze);
    return (counter == 2);
  }
715

716

717
  void PetscSolverFeti::createMatAugmentedLagrange()
718
719
720
721
722
723
724
725
726
727
  {
    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

    nOverallEdges = 0;
    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
Thomas Witkowski's avatar
Thomas Witkowski committed
728
    std::set<BoundaryObject> allEdges;
729
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
730
731
732
      if ((it->rankObj.subObj == FACE ||
	   (it->rankObj.subObj == EDGE && 
	    testWirebasketEdge(it->rankObj, feSpaces[0]))) && 
733
734
735
736
737
	  allEdges.count(it->rankObj) == 0) {
	bool dirichletOnlyEdge = true;

	DofContainer edgeDofs;
	it->rankObj.el->getAllDofs(feSpaces[0], it->rankObj, edgeDofs);
738
	
739
740
	for (DofContainer::iterator dit = edgeDofs.begin();
	     dit != edgeDofs.end(); ++dit) {
741
	  if (dirichletRows[0].count(**dit) == 0) {
742
743
744
745
746
747
748
749
750
	    dirichletOnlyEdge = false;
	    break;
	  }
	}

	if (!dirichletOnlyEdge)
	  allEdges.insert(it->rankObj);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
751
752

    nRankEdges = allEdges.size();    
753
754
755
756
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

    MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);
757
    
758
759
760
    nRankEdges *= componentSpaces.size();
    rStartEdges *= componentSpaces.size();
    nOverallEdges *= componentSpaces.size();
761
762

    MatCreateAIJ(mpiCommGlobal,
763
764
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
765
		 2, PETSC_NULL, 2, PETSC_NULL, 
766
767
768
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

769
    int rowCounter = rStartEdges;
Thomas Witkowski's avatar
Thomas Witkowski committed
770
771
    for (std::set<BoundaryObject>::iterator edgeIt = allEdges.begin(); 
	 edgeIt != allEdges.end(); ++edgeIt) {
772
      for (int component = 0; component < componentSpaces.size(); component++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
773
	DofContainer edgeDofs;
774
	edgeIt->el->getAllDofs(componentSpaces[component], *edgeIt, edgeDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
777
	
	for (DofContainer::iterator it = edgeDofs.begin();
	     it != edgeDofs.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
778
	  TEST_EXIT_DBG(isPrimal(component, **it) == false)
Thomas Witkowski's avatar
Thomas Witkowski committed
779
780
	    ("Should not be primal!\n");
	  
781
	  if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
782
783
	    continue;
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
784
	  int col = lagrangeMap.getMatIndex(component, **it);
Thomas Witkowski's avatar
Thomas Witkowski committed
785
786
	  double value = 1.0;
	  MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
787
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
790
	
	rowCounter++;
      }      
791
792
793
794
795
    }

    MatAssemblyBegin(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);

Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
796
797
798
799
800
    int nZeroRows = PetscSolverFetiDebug::testZeroRows(mat_augmented_lagrange);
    int m,n;
    MatGetSize(mat_augmented_lagrange, &m ,&n);
    MSG("Augmented lagrange matrix has %d zero rows and global size of %d %d!\n", nZeroRows, m, n);

801
802
803
804
805
806
807
808
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 05a: %.5f seconds (creation of augmented lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
    }
  }


809
  void PetscSolverFeti::createSchurPrimalKsp()
810
  {
811
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
812

Thomas Witkowski's avatar
Thomas Witkowski committed
813
    if (schurPrimalSolver == 0) {
814
815
      MSG("Create iterative schur primal solver!\n");

816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
      if (augmentedLagrange == false) {
	schurPrimalData.subSolver = subdomain;
	
	localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
	
	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getOverallDofs(), 
		       primalDofMap.getOverallDofs(),
		       &schurPrimalData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimal);	
      } else {
	schurPrimalAugmentedData.subSolver = subdomain;
833
	schurPrimalAugmentedData.nestedVec = true;
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852

	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b0, nGlobalOverallInterior);
	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b1, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalAugmentedData.tmp_vec_primal);
	lagrangeMap.createVec(schurPrimalAugmentedData.tmp_vec_lagrange);

	schurPrimalAugmentedData.mat_lagrange = &mat_lagrange;
	schurPrimalAugmentedData.mat_augmented_lagrange = &mat_augmented_lagrange;

	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges,
		       &schurPrimalAugmentedData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimalAugmented);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
853

854
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
855
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
856
857
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
858
859
      KSPSetFromOptions(ksp_schur_primal);
    } else {
860
861
      MSG("Create direct schur primal solver!\n");

862
863
      double wtime = MPI::Wtime();

864
865
866
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

867

868
      // === Create explicit matrix representation of the Schur primal system. ===
869

870
871
872
873
      if (!augmentedLagrange)
	createMatExplicitSchurPrimal();
      else
	createMatExplicitAugmentedSchurPrimal();
874
875
876


      // === Create KSP solver object and set appropriate solver options. ====
877

878
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
879
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
880
		      SAME_NONZERO_PATTERN);
881
882
883
884
885
886
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPPREONLY);
      PC pc_schur_primal;      
      KSPGetPC(ksp_schur_primal, &pc_schur_primal);
      PCSetType(pc_schur_primal, PCLU);
      PCFactorSetMatSolverPackage(pc_schur_primal, MATSOLVERMUMPS);
887
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
888

889
890
891

      // === And finally print timings, if required. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
892
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
893
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
894
895
896
897
898
899
900
901
902
903
	MatInfo minfo;
	MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
	MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);
	
	MSG("FETI-DP timing 06: %.5f seconds (creation of schur primal matrix)\n",
	    MPI::Wtime() - wtime);

	wtime = MPI::Wtime();
	KSPSetUp(ksp_schur_primal);
	KSPSetUpOnBlocks(ksp_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
904
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
905
906
907
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
908
    }
909
910
911
  }


912
913
914
915
916
917
918
919
920
921
922
923
  void PetscSolverFeti::createMatExplicitSchurPrimal()
  {
    FUNCNAME("PetscSolverFeti::createMatExplicitSchurPrimal()");

    int creationMode = 0;
    Parameters::get("parallel->feti->schur primal creation mode", creationMode);
    if (creationMode == 0) {
      // matK = inv(A_BB) A_BPi
      Mat matK;
      petsc_helper::blockMatMatSolve(subdomain->getSolver(), 
				     subdomain->getMatInteriorCoarse(),
				     matK);
924
     
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
      
      // mat_schur_primal = A_PiPi - A_PiB inv(A_BB) A_BPi