PetscSolverFeti.cc 46.4 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
      subDomainSolver(NULL),
193
194
195
      meshLevel(0),
      rStartInterior(0),
      nGlobalOverallInterior(0)
196
197
198
199
200
201
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

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


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

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

233
234
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

235
236
237
238
239
240
241
242
243
    if (subDomainSolver == NULL) {
      if (meshLevel == 0) {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);
      } else {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, 
			      levelData.getMpiComm(meshLevel - 1),
			      levelData.getMpiComm(meshLevel));
244
	subDomainSolver->setLevel(meshLevel);
245
246
      }
    }
247

248
249
    primalDofMap.init(levelData,
		      feSpaces, meshDistributor->getFeSpaces(), 
250
		      true, true);
251
252
    primalDofMap.setComputeMatIndex(true);
    
253
254
    dualDofMap.init(levelData, 
		    feSpaces, meshDistributor->getFeSpaces(),
255
		    false, false);
256
257
    dualDofMap.setComputeMatIndex(true);

258
259
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
260
		     true, true);
261
    lagrangeMap.setComputeMatIndex(true);
262
263
264
265
266
267
268
269
270
271
272
273

    if (meshLevel == 0) {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       false, false);
      localDofMap.setComputeMatIndex(true);
    } else {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       true, true);
      localDofMap.setComputeMatIndex(true);
    }
274
275

    if (fetiPreconditioner == FETI_DIRICHLET) {
276
277
278
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

279
280
      interiorDofMap.init(levelData, 
			  feSpaces,  meshDistributor->getFeSpaces(),
281
			  false, false);
282
283
      interiorDofMap.setComputeMatIndex(true);
    }
284
285
286
287
288
289
290
291
292
293
  }


  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");
294

295
296
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

297
298
299
300
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
301
302
303
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.clear();

304
305
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
306

307
308
309
310
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);    
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
311
312
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
313
314
315
316

    if (meshLevel > 0)
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

317
318
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
319

320
321
      createPrimals(feSpace);  

322
      createDuals(feSpace);      
323

324
325
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
328
329

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
330
    if (fetiPreconditioner != FETI_NONE)
331
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
332

333
334
335
336
337
338
339
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
    lagrangeMap.update();


340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372

    // === ===

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

      MSG("RANK %d FROM %d\n",
	  levelData.getMpiComm(1).Get_rank(),
	  levelData.getMpiComm(1).Get_size());

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

      mpi::getDofNumbering(mpiComm, groupRowsInterior,
			   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);

      MSG("COMM TEST FETI-DP: %d %d %d %d %d\n",
	  levelData.getMpiComm(1).Get_size(),
	  localDofMap.getRankDofs(),
	  localDofMap.getOverallDofs(),
	  nGlobalOverallInterior, rStartInterior);
    }

373
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
Thomas Witkowski's avatar
Thomas Witkowski committed
374
375
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
376
      MSG("FETI-DP data for %d-ith FE space:\n", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
377

378
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
379
380
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
381
      
382
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
383
384
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
385

386
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
387
388
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
389

390
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
391
392
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
393
    }
394
395
396


    // If multi level test, inform sub domain solver about coarse space.
397
    subDomainSolver->setDofMapping(&primalDofMap, &localDofMap);
398
399
400
  }


401
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
402
  {
403
    FUNCNAME("PetscSolverFeti::createPrimals()");  
404

405
406
407
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

408
409
    /// Set of DOF indices that are considered to be primal variables.
    
410
    DofContainerSet& vertices = 
411
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
412
413

    DofIndexSet primals;
414
    for (DofContainerSet::iterator it = vertices.begin(); 
415
416
417
418
419
420
421
422
423
424
	 it != vertices.end(); ++it) {
      WorldVector<double> c;
      feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
      MSG("COORD %f %f\n", c[0], c[1]);
      double e = 1e-8;
      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))
	primals.insert(**it);
    }
425

426
427
428
429

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

430
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
431
432
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
433
      else
434
  	primalDofMap[feSpace].insertNonRankDof(*it);
435
436
437
  }


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

442
443
444
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
445
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
446
447
448
449

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
450
451
452
453
454
455
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
	    dualDofMap[feSpace].insertRankDof(**it);
	}	  
456
457
458
459
460
461
462
  }

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

463
464
    boundaryDofRanks[feSpace].clear();

465
    if (dualDofMap[feSpace].nLocalDofs == 0)
466
467
      return;

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

471
472
473
474
475
476
477
478
479
    if (meshLevel == 0) {
      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				feSpace); 
	   !it.end(); it.nextRank()) {
	for (; !it.endDofIter(); it.nextDof()) {
	  if (!isPrimal(feSpace, it.getDofIndex())) {
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
	  }
480
481
	}
      }
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
    } else {
      StdMpi<vector<int> > stdMpi(mpiComm);

      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()) {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())) {
	    MSG("IN SD DOF\n");
	    subdomainRankDofs.push_back(1);
	  } else {
	    MSG("NOT IN SD DOF\n");
	    subdomainRankDofs.push_back(0);
	  }
	}
	
	MSG("SEND SIZE %d TO RANK %d\n", subdomainRankDofs.size(), it.getRank());

	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); 
	   !it.end(); it.nextRank()) {
	for (; !it.endDofIter(); it.nextDof()) {
	  if (!isPrimal(feSpace, it.getDofIndex())) {
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

	    if (it.getDofCounter() >= stdMpi.getRecvData(it.getRank()).size())
	      MSG("VERY BIG SHIT: %d from %d\n", stdMpi.getRecvData(it.getRank()).size(), it.getRank());

	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()])
	      boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
	  }
	}
      }
    }

    for (DofIndexToPartitions::iterator it = boundaryDofRanks[feSpace].begin();
	 it != boundaryDofRanks[feSpace].end(); ++it)
      MSG("DOF %d IS DUAL IN %d RANKS\n", it->first, it->second.size());
534
	
535
536
537
538

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

539
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
540

541
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
542
543
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
544
	if (!isPrimal(feSpace, it.getDofIndex()))
545
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
546
547
548

    stdMpi.updateSendDataSize();

549
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
550
	 !it.end(); it.nextRank()) {
551
      bool recvFromRank = false;
552
      for (; !it.endDofIter(); it.nextDof()) {
553
	if (!isPrimal(feSpace, it.getDofIndex())) {
554
555
556
	  recvFromRank = true;
	  break;
	}
557
      }
558
559

      if (recvFromRank)
560
	stdMpi.recv(it.getRank());
561
    }
562

563
564
    stdMpi.startCommunication();

565
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
566
	 !it.end(); it.nextRank()) {
567
      int i = 0;
568
      for (; !it.endDofIter(); it.nextDof())
569
	if (!isPrimal(feSpace, it.getDofIndex()))
570
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
571
	    stdMpi.getRecvData(it.getRank())[i++];
572
573
    }

574
575
576
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

577
    int nRankLagrange = 0;
578
    DofMap& dualMap = dualDofMap[feSpace].getMap();
579
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
580

581
582
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
583
	int degree = boundaryDofRanks[feSpace][it->first].size();
584
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
585
      } else {
586
	lagrangeMap[feSpace].insertNonRankDof(it->first);
587
588
      }
    }
589
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
590
591
592
  }


593
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
594
  {
595
    FUNCNAME("PetscSolverFeti::createIndexB()");
596

597
    DOFAdmin* admin = feSpace->getAdmin();
598
599
600
601

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

603
    int nLocalInterior = 0;    
604
    for (int i = 0; i < admin->getUsedSize(); i++) {
605
      if (admin->isDofFree(i) == false && 
606
	  isPrimal(feSpace, i) == false &&
607
	  isDual(feSpace, i) == false) {
608
609
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
610

611
612
613
614
615
616
617
618
619
	  if (fetiPreconditioner != FETI_NONE)
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
620
621
622

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
623
	}
624
      }
625
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
626
    
627
628
    // === And finally, add the global indicies of all dual nodes. ===

629
630
631
632
633
634
635
636
637
638
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      if (meshLevel == 0)
	localDofMap[feSpace].insertRankDof(it->first);
      else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
639
640
641
  }


642
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
643
644
645
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

646
647
    // === Create distributed matrix for Lagrange constraints. ===

648
    MatCreateMPIAIJ(mpiComm,
649
650
651
		    lagrangeMap.getRankDofs(), 
		    localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), 
652
		    nGlobalOverallInterior,
653
654
655
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

656
657
658
659
660
661
662
    // === 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
663
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
664
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
665

666
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
667
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
668
669
670
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
671
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Thomas Witkowski committed
672
673
	
	// Copy set of all ranks that contain this dual node.
674
675
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
678
679
680
681
	// 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) {
682
	      int colIndex = localDofMap.getMatIndex(k, it->first);
683
684
685
686

	      if (meshLevel > 0)
		colIndex += rStartInterior;

687
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
688
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
689
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
690
	    index++;	      
691
692
693
694
695
696
697
698
699
700
	  }
	}
      }
    }

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


701
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
702
  {
703
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
704

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

708
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
709
      
710
      VecCreateMPI(mpiComm, 
711
		   localDofMap.getRankDofs(), 
712
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
713
		   &(schurPrimalData.tmp_vec_b));
714
      VecCreateMPI(mpiComm, 
715
716
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
717
718
		   &(schurPrimalData.tmp_vec_primal));

719
      MatCreateShell(mpiComm,
720
721
722
723
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
724
725
726
727
728
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
729
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
730
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
731
732
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
733
734
      KSPSetFromOptions(ksp_schur_primal);
    } else {
735
736
      MSG("Create direct schur primal solver!\n");

737
738
      double wtime = MPI::Wtime();

739
740
741
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

742
743
744
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
745

Thomas Witkowski's avatar
Thomas Witkowski committed
746
      Mat matBPi;
747
      MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
748
		      nRowsRankB, nRowsRankPrimal, 
749
		      nGlobalOverallInterior, nRowsOverallPrimal,
Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
752

753
754
755
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
756

757
758
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

759
      for (int i = 0; i < nRowsRankB; i++) {
760
	PetscInt row = localDofMap.getStartDofs() + i;
761
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
762
763
764
765
766

	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]));

767
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
768
769
      }

770
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
771
		primalDofMap.getLocalDofs())
772
	("Should not happen!\n");
773

774
775
776
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
777
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
778
779
780
781
782
783
784
785

 	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);

786
	subDomainSolver->solve(tmpVec, tmpVec);
787

Thomas Witkowski's avatar
Thomas Witkowski committed
788
	PetscScalar *tmpValues;
789
	VecGetArray(tmpVec, &tmpValues);
790
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
791
	  MatSetValue(matBPi, 
792
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
793
794
795
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
796
797
798
799
800
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
801
802
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
803
804

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
805
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
806
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
807
808

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
809
      MatDestroy(&matBPi);
810
811

      MatInfo minfo;
812
      MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
813
814
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

815
      KSPCreate(mpiComm, &ksp_schur_primal);
816
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
817
		      SAME_NONZERO_PATTERN);
818
819
820
821
822
823
      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);
824
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
825

826
827
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
828
    }
829
830
831
832
833
834
835
  }


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

836
    if (schurPrimalSolver == 0) {
837
      schurPrimalData.subSolver = NULL;
838

839
840
841
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
842
843
844
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
845
846
847
  }


848
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
849
850
851
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

852
853
    // === Create FETI-DP solver object. ===

854
    fetiData.mat_lagrange = &mat_lagrange;
855
    fetiData.subSolver = subDomainSolver;
856
    fetiData.ksp_schur_primal = &ksp_schur_primal;
857

858
    VecCreateMPI(mpiComm, 
859
		 localDofMap.getRankDofs(), 
860
		 nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
861
		 &(fetiData.tmp_vec_b));
862
    VecCreateMPI(mpiComm,
863
864
		 lagrangeMap.getRankDofs(), 
		 lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
865
		 &(fetiData.tmp_vec_lagrange));
866
    VecCreateMPI(mpiComm, 
867
868
		 primalDofMap.getRankDofs(), 
		 primalDofMap.getOverallDofs(),
869
		 &(fetiData.tmp_vec_primal));
870

871
    MatCreateShell(mpiComm,
872
873
874
875
		   lagrangeMap.getRankDofs(), 
		   lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), 
		   lagrangeMap.getOverallDofs(),
876
		   &fetiData, &mat_feti);
877
878
879
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


880
    KSPCreate(mpiComm, &ksp_feti);
881
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
882
883
884
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
885
    KSPSetFromOptions(ksp_feti);
886
887


888
    // === Create FETI-DP preconditioner object. ===
889

890
891
892
893
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
894

895
    switch (fetiPreconditioner) {
896
897
898
899
    case FETI_DIRICHLET:      
      TEST_EXIT_DBG(meshLevel == 0)
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

900
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
901
902
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
903
904
905
906
907
908
      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);
909
910
911
912
913
914
915
916
917
      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;
      
918
      VecCreateMPI(mpiComm, 
919
		   localDofMap.getRankDofs(),
920
		   nGlobalOverallInterior,
921
		   &(fetiDirichletPreconData.tmp_vec_b));      
922
923
924
925
926
927
      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));
928

929
930
931
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

932
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
933
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
934
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
935
	  DegreeOfFreedom d = it->first;	  
936
937
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
938
939
940
941
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

942
943
944
945
946
947
948
949
950
951
952
      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
953
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
954
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
955
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
956
	  DegreeOfFreedom d = it->first;
957
958
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
Thomas Witkowski's avatar
Thomas Witkowski committed
959
960
961
962
	  fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

963
      VecCreateMPI(mpiComm, 
964
965
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
966
		   &(fetiLumpedPreconData.tmp_vec_b));