PetscSolverFeti.cc 41 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"
14
#include "parallel/PetscSolverFeti.h"
15
#include "parallel/PetscSolverFetiStructs.h"
16
17
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
18
#include "parallel/SubDomainSolver.h"
19
#include "io/VtkWriter.h"
20
21
22
23
24

namespace AMDiS {

  using namespace std;

25
26
27
28
29
30
31
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
32
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
33

34
35
36
37
    MatMult(data->subSolver->getMatIntCoarse(), x, data->tmp_vec_b);
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarseCoarse(), x, y);
38
39
40
41
42
43
44
45
46
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
47
48
    //    F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
49
50
51

    void *ctx;
    MatShellGetContext(mat, &ctx);
52
    FetiData* data = static_cast<FetiData*>(ctx);
53

Thomas Witkowski's avatar
Thomas Witkowski committed
54
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
55
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
57

58
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
59
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
60
   MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
61
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
63

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
65
66
67
68
69

    return 0;
  }


70
  // y = PC * x
71
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
72
  {
73
    // Get data for the preconditioner
74
75
    void *ctx;
    PCShellGetContext(pc, &ctx);
76
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
77

78
    // Multiply with scaled Lagrange constraint matrix.
79
80
81
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


82
    // === Restriction of the B nodes to the boundary nodes. ===
83

84
85
86
87
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
88

89
90
91
92
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

93
94
95
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
96
97

    VecRestoreArray(data->tmp_vec_b, &local_b);
98
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
99
100


101
    // === K_DD - K_DI inv(K_II) K_ID ===
102

103
    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
104

105
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
106
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
107
108
109
110
111
112
113
114
115
116
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

    VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);


    // === Prolongation from local dual nodes to B nodes.

    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

117
118
119
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
142
143


144
    // === Restriction of the B nodes to the boundary nodes. ===
145

146
147
148
149
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
150

151
    PetscScalar *local_b, *local_duals;
152
    VecGetArray(data->tmp_vec_b, &local_b);
153
    VecGetArray(data->tmp_vec_duals0, &local_duals);
154

Thomas Witkowski's avatar
Thomas Witkowski committed
155
156
157
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
158
159

    VecRestoreArray(data->tmp_vec_b, &local_b);
160
161
162
163
164
165
166
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);

167

168
    // === Prolongation from local dual nodes to B nodes.
169

170
171
172
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
175
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
176

177
178
    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
179

180
181

    // Multiply with scaled Lagrange constraint matrix.
182
183
184
185
186
187
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


188
189
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
190
      schurPrimalSolver(0),
191
      multiLevelTest(false),
192
193
      subDomainSolver(NULL),
      meshLevel(0)
194
195
196
197
198
199
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
200
      MSG("Create FETI-DP solver with no preconditioner!\n");
201
202
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
203
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
204
205
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
206
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
207
208
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
209
210
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
211
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
212
213
214
215
216

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

    Parameters::get("parallel->multi level test", multiLevelTest);
219
220
    if (multiLevelTest)
      meshLevel = 1;
221
222
223
  }


224
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
225
  {
226
227
    FUNCNAME("PetscSolverFeti::initialize()");

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

231
232
233
234
235
236
237
238
239
240
241
242
243
    if (subDomainSolver == NULL) {
      if (meshLevel == 0) {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);
      } else {
	MeshLevelData& levelData = meshDistributor->getMeshLevelData();

	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, 
			      levelData.getMpiComm(meshLevel - 1),
			      levelData.getMpiComm(meshLevel));
      }
    }
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

    primalDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(), 
		      true, true);
    dualDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		    false, false);
    lagrangeMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		     true, true);
    localDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		     false, false);
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.init(mpiComm, feSpaces,  meshDistributor->getFeSpaces(),
			  false, false);
  }


  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
266

267
268
269
270
271
272
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
    interiorDofMap.clear();

273
274
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
275

276
277
      createPrimals(feSpace);      
      createDuals(feSpace);
278
      createLagrange(feSpace);      
279
280
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
281

282
283
284
285
286
    primalDofMap.setDofComm(meshDistributor->getSendDofs(),
			    meshDistributor->getRecvDofs());
    lagrangeMap.setDofComm(meshDistributor->getSendDofs(), 
			   meshDistributor->getRecvDofs());

Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
289
290
    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();
291
    if (fetiPreconditioner != FETI_NONE)
292
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
295

    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
296
      MSG("FETI-DP data for %d-ith FE space:\n", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
297

298
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
      
301
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
302
	  dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
303

304
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
305
	  lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
306
307
308
309

      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
310
    }
311
312
313


    // If multi level test, inform sub domain solver about coarse space.
314
    subDomainSolver->setDofMapping(&primalDofMap, &localDofMap);
315
316
317
  }


318
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
319
  {
320
    FUNCNAME("PetscSolverFeti::createPrimals()");  
321

322
323
324
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

325
326
327
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
328
    DofContainerSet& vertices = 
329
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
330
331
332
    for (DofContainerSet::iterator it = vertices.begin(); 
	 it != vertices.end(); ++it)
      primals.insert(**it);
333

334
335
336
337

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

338
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
339
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
340
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
341
342
      else
  	primalDofMap[feSpace].insert(*it);
343
344
345
  }


346
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
347
348
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
349

350
351
352
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
353
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
354
355
356
357
358
359
360
361
362
363
364
365

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
	dualDofMap[feSpace].insertRankDof(**it);
  }

  
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

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

369
    boundaryDofRanks[feSpace].clear();
370

371
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), meshLevel, feSpace); 
372
373
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
374
	if (!isPrimal(feSpace, it.getDofIndex())) {
375
376
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
377
378
	}
      }
379
	
380
381
382
383

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

384
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
385

386
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), meshLevel, feSpace);
387
388
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
389
	if (!isPrimal(feSpace, it.getDofIndex()))
390
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
391
392
393

    stdMpi.updateSendDataSize();

394
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), meshLevel, feSpace); 
395
	 !it.end(); it.nextRank()) {
396
      bool recvFromRank = false;
397
      for (; !it.endDofIter(); it.nextDof()) {
398
	if (!isPrimal(feSpace, it.getDofIndex())) {
399
400
401
	  recvFromRank = true;
	  break;
	}
402
      }
403
404

      if (recvFromRank)
405
	stdMpi.recv(it.getRank());
406
    }
407

408
409
    stdMpi.startCommunication();

410
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), meshLevel, feSpace); 
411
	 !it.end(); it.nextRank()) {
412
      int i = 0;
413
      for (; !it.endDofIter(); it.nextDof())
414
	if (!isPrimal(feSpace, it.getDofIndex()))
415
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
416
	    stdMpi.getRecvData(it.getRank())[i++];
417
418
    }

419
420
421
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

422
    int nRankLagrange = 0;
423
424
    DofMap& dualMap = dualDofMap[feSpace].getMap();
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
425
426

      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
427
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
428
	int degree = boundaryDofRanks[feSpace][it->first].size();
429
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
      } else {
	lagrangeMap[feSpace].insert(it->first);
432
433
      }
    }
434
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
435
436
437
  }


438
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
439
  {
440
    FUNCNAME("PetscSolverFeti::createIndexB()");
441

442
    DOFAdmin* admin = feSpace->getAdmin();
443
444
445
446

    // === 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.                                 ===
447

448
    int nLocalInterior = 0;    
449
    for (int i = 0; i < admin->getUsedSize(); i++) {
450
      if (admin->isDofFree(i) == false && 
451
	  isPrimal(feSpace, i) == false &&
452
	  dualDofMap[feSpace].isSet(i) == false) {
453
454
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);

455
	if (fetiPreconditioner != FETI_NONE)
456
457
458
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);

	nLocalInterior++;
459
      }
460
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
461
    
462
463
    // === And finally, add the global indicies of all dual nodes. ===

464
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
465
	 it != dualDofMap[feSpace].getMap().end(); ++it)
466
      localDofMap[feSpace].insertRankDof(it->first);
467
468
469
  }


470
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
471
472
473
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

474
475
    // === Create distributed matrix for Lagrange constraints. ===

476
    MatCreateMPIAIJ(mpiComm,
477
478
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
479
480
481
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

482
483
484
485
486
487
488
    // === 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.                                                     ===

Thomas Witkowski's avatar
Thomas Witkowski committed
489
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
490
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
491

492
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
493
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
494
495
496
497
498
499
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
	int index = lagrangeMap.getMatIndex(k, it->first);
	
	// Copy set of all ranks that contain this dual node.
500
501
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
504
505
506
507
	// Number of ranks that contain this dual node.
	int degree = W.size();
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
508
	      int colIndex = localDofMap.getMatIndex(k, it->first);
509
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
510
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
511
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
512
	    index++;	      
513
514
515
516
517
518
519
520
521
522
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
  }


523
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
524
525
526
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

530
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
531
      
532
      VecCreateMPI(mpiComm, 
533
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
534
		   &(schurPrimalData.tmp_vec_b));
535
      VecCreateMPI(mpiComm, 
536
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
537
538
		   &(schurPrimalData.tmp_vec_primal));

539
      MatCreateShell(mpiComm,
540
541
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
542
543
544
545
546
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
547
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
548
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
549
550
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
      KSPSetFromOptions(ksp_schur_primal);
    } else {
553
554
      MSG("Create direct schur primal solver!\n");

555
556
      double wtime = MPI::Wtime();

557
558
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
559
560
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
561

Thomas Witkowski's avatar
Thomas Witkowski committed
562
      Mat matBPi;
563
      MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
564
565
566
567
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
568

569
570
571
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
572

573
574
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

575
576
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
577
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
578
579
580
581
582

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

583
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
584
585
      }

586
587
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
		primalDofMap.getLocalDofs())
588
	("Should not happen!\n");
589

590
591
592
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
593
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
594
595
596
597
598
599
600
601

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);

602
	subDomainSolver->solve(tmpVec, tmpVec);
603

Thomas Witkowski's avatar
Thomas Witkowski committed
604
	PetscScalar *tmpValues;
605
	VecGetArray(tmpVec, &tmpValues);
606
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
607
	  MatSetValue(matBPi, 
608
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
609
610
611
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
612
613
614
615
616
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
617
618
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
619
620

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
621
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
622
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
623
624

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
625
      MatDestroy(&matBPi);
626
627

      MatInfo minfo;
628
      MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
629
630
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

631
      KSPCreate(mpiComm, &ksp_schur_primal);
632
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
633
		      SAME_NONZERO_PATTERN);
634
635
636
637
638
639
      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);
640
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
641

642
643
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
644
    }
645
646
647
648
649
650
651
  }


  void PetscSolverFeti::destroySchurPrimalKsp()
  {
    FUNCNAME("PetscSolverFeti::destroySchurPrimal()");

652
    if (schurPrimalSolver == 0) {
653
      schurPrimalData.subSolver = NULL;
654

655
656
657
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
658
659
660
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
661
662
663
  }


664
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
665
666
667
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

668
669
    // === Create FETI-DP solver object. ===

670
    fetiData.mat_lagrange = &mat_lagrange;
671
    fetiData.subSolver = subDomainSolver;
672
    fetiData.ksp_schur_primal = &ksp_schur_primal;
673

674
    VecCreateMPI(mpiComm, 
675
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
676
		 &(fetiData.tmp_vec_b));
677
    VecCreateMPI(mpiComm,
678
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
679
		 &(fetiData.tmp_vec_lagrange));
680
    VecCreateMPI(mpiComm, 
681
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
682
		 &(fetiData.tmp_vec_primal));
683

684
    MatCreateShell(mpiComm,
685
686
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
687
		   &fetiData, &mat_feti);
688
689
690
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


691
    KSPCreate(mpiComm, &ksp_feti);
692
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
693
694
695
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
696
    KSPSetFromOptions(ksp_feti);
697
698


699
    // === Create FETI-DP preconditioner object. ===
700

701
702
703
704
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
705

706
707
708
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
709
710
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
711
712
713
714
715
716
      KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
      KSPSetType(ksp_interior, KSPPREONLY);
      PC pc_interior;
      KSPGetPC(ksp_interior, &pc_interior);
      PCSetType(pc_interior, PCLU);
      PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK);
717
718
719
720
721
722
723
724
725
      KSPSetFromOptions(ksp_interior);
            
      fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
      fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
      fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
      fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
      fetiDirichletPreconData.ksp_interior = &ksp_interior;
      
726
      VecCreateMPI(mpiComm, 
727
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
728
		   &(fetiDirichletPreconData.tmp_vec_b));      
729
730
731
732
733
734
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_interior));
735
736

      for (unsigned int i = 0; i < feSpaces.size(); i++) {
737
738
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
739
740
741
742
743
744
745
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

746
747
748
749
750
751
752
753
754
755
756
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
      
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;

Thomas Witkowski's avatar
Thomas Witkowski committed
757
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
758
759
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
760
761
762
763
764
765
766
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

767
      VecCreateMPI(mpiComm, 
768
769
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
770
		   &(fetiLumpedPreconData.tmp_vec_b));
771
772
773
774
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
775
776
777
778
779
780

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
781
782
      break;
    default:
783
784
      break;
    }
785
786
787
788
789
790
791
  }
  

  void PetscSolverFeti::destroyFetiKsp()
  {
    FUNCNAME("PetscSolverFeti::destroyFetiKsp()");

792
793
    // === Destroy FETI-DP solver object. ===

794
    fetiData.mat_lagrange = PETSC_NULL;
795
    fetiData.subSolver = NULL;
796
    fetiData.ksp_schur_primal = PETSC_NULL;
797

Thomas Witkowski's avatar
Thomas Witkowski committed
798
799
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
800
    VecDestroy(&fetiData.tmp_vec_primal);
801
802
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
803
804


805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
    // === Destroy FETI-DP preconditioner object. ===

    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPDestroy(&ksp_interior);

      fetiDirichletPreconData.mat_lagrange_scaled = NULL;
      fetiDirichletPreconData.mat_interior_interior = NULL;
      fetiDirichletPreconData.mat_duals_duals = NULL;
      fetiDirichletPreconData.mat_interior_duals = NULL;
      fetiDirichletPreconData.mat_duals_interior = NULL;
      fetiDirichletPreconData.ksp_interior = NULL;
      
      VecDestroy(&fetiDirichletPreconData.tmp_vec_b);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals0);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals1);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_interior);
      MatDestroy(&mat_lagrange_scaled);
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = NULL;
      fetiLumpedPreconData.mat_duals_duals = NULL;

      VecDestroy(&fetiLumpedPreconData.tmp_vec_b);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1);
      break;
833
834
    default:
      break;
835
    }
836
837
838
839
840
841
842
843
844
  }


  void PetscSolverFeti::recoverSolution(Vec &vec_sol_b,
					Vec &vec_sol_primal,
					SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::recoverSolution()");

845
    // === Get local part of the solution for B variables. ===
846
847
848
849
850

    PetscScalar *localSolB;
    VecGetArray(vec_sol_b, &localSolB);


851
852
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
853
    
Thomas Witkowski's avatar
Thomas Witkowski committed
854
855
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);

856
    vector<PetscInt> globalIsIndex, localIsIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
857
858
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());
859
860

    {
Thomas Witkowski's avatar
Thomas Witkowski committed
861
862
      int cnt = 0;
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
863
864
	DofMap& dofMap = primalDofMap[feSpaces[i]].getMap();
	for (DofMap::iterator it = dofMap.begin();it != dofMap.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
	  globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
	  localIsIndex.push_back(cnt++);
867
868
	}
      }
869
      
Thomas Witkowski's avatar
Thomas Witkowski committed
870
      TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs())("Should not happen!\n");
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
    }
    
    IS globalIs, localIs;
    ISCreateGeneral(PETSC_COMM_SELF, 
		    globalIsIndex.size(), 
		    &(globalIsIndex[0]),
		    PETSC_USE_POINTER,
		    &globalIs);

    ISCreateGeneral(PETSC_COMM_SELF, 
		    localIsIndex.size(), 
		    &(localIsIndex[0]),
		    PETSC_USE_POINTER,
		    &localIs);

    Vec local_sol_primal;
    VecCreateSeq(PETSC_COMM_SELF, localIsIndex.size(), &local_sol_primal);

    VecScatter primalScatter;
    VecScatterCreate(vec_sol_primal, globalIs, local_sol_primal, localIs, &primalScatter);
    VecScatterBegin(primalScatter, vec_sol_primal, local_sol_primal, 
		    INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(primalScatter, vec_sol_primal, local_sol_primal, 
		  INSERT_VALUES, SCATTER_FORWARD);

896
897
898
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
899
900
901
902

    PetscScalar *localSolPrimal;
    VecGetArray(local_sol_primal, &localSolPrimal);