PetscSolverFeti.cc 54.9 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/PetscSolverGlobalMatrix.h"
19
#include "io/VtkWriter.h"
20
21
22
23
24

namespace AMDiS {

  using namespace std;

25
26
27
28
29
30
  double FetiTimings::fetiSolve = 0.0;
  double FetiTimings::fetiSolve01 = 0.0;
  double FetiTimings::fetiSolve02 = 0.0;
  double FetiTimings::fetiPreconditioner = 0.0;


31
32
33
34
35
36
37
  // 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);
38
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
39

40
41
42
43
    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);
44
45
46
47
48
49
50
51
52
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
53
54
    FUNCNAME("petscMultMatFeti()");

55
56
    //    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)
57

58
59
    double wtime = MPI::Wtime();

60
61
    void *ctx;
    MatShellGetContext(mat, &ctx);
62
    FetiData* data = static_cast<FetiData*>(ctx);
63

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
65
66

    double wtime01 = MPI::Wtime();
67
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
68
69
70

    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
71
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
72

73
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
74
75

    wtime01 = MPI::Wtime();
76
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
77
78
    FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);

79
    MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
80
81

    wtime01 = MPI::Wtime();
82
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
83
84
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
85
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
86

Thomas Witkowski's avatar
Thomas Witkowski committed
87
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
88

89
90
    FetiTimings::fetiSolve += (MPI::Wtime() - wtime);

91
92
93
94
    return 0;
  }


95
  // y = PC * x
96
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
97
  {
98
99
    double wtime = MPI::Wtime();

100
    // Get data for the preconditioner
101
102
    void *ctx;
    PCShellGetContext(pc, &ctx);
103
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
104

105
    // Multiply with scaled Lagrange constraint matrix.
106
107
108
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


109
    // === Restriction of the B nodes to the boundary nodes. ===
110

111
112
113
114
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
115

116
117
118
119
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

120
121
122
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
123
124

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


128
    // === K_DD - K_DI inv(K_II) K_ID ===
129

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

132
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
133
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
134
135
136
137
138
139
140
141
142
143
    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);

144
145
146
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
147
148
149
150
151
152
153
154

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

155
156
    FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);

157
158
159
160
161
162
163
164
165
166
167
168
169
170
    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);
171
172


173
    // === Restriction of the B nodes to the boundary nodes. ===
174

175
176
177
178
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
179

180
    PetscScalar *local_b, *local_duals;
181
    VecGetArray(data->tmp_vec_b, &local_b);
182
    VecGetArray(data->tmp_vec_duals0, &local_duals);
183

Thomas Witkowski's avatar
Thomas Witkowski committed
184
185
186
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
187
188

    VecRestoreArray(data->tmp_vec_b, &local_b);
189
190
191
192
193
194
195
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

196

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

199
200
201
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

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

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

209
210

    // Multiply with scaled Lagrange constraint matrix.
211
212
213
214
215
216
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


217
218
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
219
      schurPrimalSolver(0),
220
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
221
      subdomain(NULL),
222
223
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
      nGlobalOverallInterior(0),
      printTimings(false)
226
227
228
229
230
231
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
232
      MSG("Create FETI-DP solver with no preconditioner!\n");
233
234
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
235
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
236
237
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
238
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
239
240
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
241
242
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
243
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
244
245
246
247
248

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

    Parameters::get("parallel->multi level test", multiLevelTest);
251
252
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
253
254

    Parameters::get("parallel->print timings", printTimings);
255
256
257
  }


258
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
259
  {
260
261
    FUNCNAME("PetscSolverFeti::initialize()");

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

265
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
266
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
267

Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
270

271
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
272
	subdomain->setMeshDistributor(meshDistributor, 
273
				      mpiCommGlobal, mpiCommLocal);
274
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
275
	subdomain->setMeshDistributor(meshDistributor, 
276
277
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
278
	subdomain->setLevel(meshLevel);
279
280
      }
    }
281

282
283
284
285
    primalDofMap.init(levelData, feSpaces, uniqueFe);   
    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
    lagrangeMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
286

287
288
    if (fullInterface != NULL)
      interfaceDofMap.init(levelData, fullInterface);
Thomas Witkowski's avatar
Thomas Witkowski committed
289

290
    if (fetiPreconditioner != FETI_NONE) {
291
292
293
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

294
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
295
    }
296
297
298
299
300
301
302
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
    double timeCounter = MPI::Wtime();

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

309
310
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

311
312
313
314
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
315
    if (fetiPreconditioner != FETI_NONE)
316
317
      interiorDofMap.clear();

318
319
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
320

321
322
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
323
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
324
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
325
326
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
327

328
329
330
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
331
332
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

333
334
335
336
337
338
    if (fullInterface != NULL) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

339
340
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
341

342
343
      createPrimals(feSpace);  

344
345
346
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
347

348
349
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
352
353

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
354
    if (fetiPreconditioner != FETI_NONE)
355
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
356

357
358
359
    if (fullInterface != NULL)
      interfaceDofMap.update();

360
361
362
363
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
364

365
366
367
    lagrangeMap.update();


368
369
370
371
372
373
374
375
376
377
378
379
    // === ===

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

380
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
381
382
383
384
385
386
387
388
389
			   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);
    }

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

395
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
396
397
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
398
      
399
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
400
401
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
402

403
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
404
405
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
406

407
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
408
409
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
410
    }
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
414
    subdomain->setDofMapping(&localDofMap, &primalDofMap);

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
415
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
416
      timeCounter = MPI::Wtime() - timeCounter;
417
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
418
419
	  timeCounter);
    }
420
421
422
  }


423
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
424
  {
425
    FUNCNAME("PetscSolverFeti::createPrimals()");  
426

427
428
429
    if (feSpace == fullInterface)
      return;

430
431
432
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    // Set of DOF indices that are considered to be primal variables.
434
    DofContainerSet& vertices = 
435
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
436
437

    DofIndexSet primals;
438
    for (DofContainerSet::iterator it = vertices.begin(); 
439
440
	 it != vertices.end(); ++it) {
      double e = 1e-8;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
441
      if (meshLevel == 0) {
442
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
443
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
446
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
447
448
449
450
451
452
453
	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);
	}
      }
454
    }
455

456
457
458
459

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

460
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
461
462
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
463
      else
464
  	primalDofMap[feSpace].insertNonRankDof(*it);
465
466
467
  }


468
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
469
470
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
471

472
473
474
    if (feSpace == fullInterface)
      return;

475
476
477
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
478
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
479
480
481
482

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
483
484
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
485
486
487
488
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
489
490
491
  }

  
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

    if (feSpace != fullInterface)
      return;

    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	interfaceDofMap[feSpace].insertRankDof(*it);
      else
	interfaceDofMap[feSpace].insertNonRankDof(*it);
  }


511
512
513
514
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

515
516
517
    if (feSpace == fullInterface)
      return;

518
519
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
522
    // 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
523
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
524

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
525
    if (meshLevel > 0) {
526
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
527
528
529
530
531
532
533
534
535

      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()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
536
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
537
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
538
	  else
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
	    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
554
555
556
557
558
559
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex()))
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
560

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
561
562
563
564
565
566
567
568
569
570
571
572
573
574
    if (dualDofMap[feSpace].nLocalDofs == 0)
      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)).         ===

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

575
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
576
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
577
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
578
579
580
581
	}
      }
    }

582
583
584
585

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

586
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
587

588
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
589
590
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
591
	if (!isPrimal(feSpace, it.getDofIndex()))
592
593
594
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
595
596
597

    stdMpi.updateSendDataSize();

598
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
599
	 !it.end(); it.nextRank()) {
600
      bool recvFromRank = false;
601
      for (; !it.endDofIter(); it.nextDof()) {
602
	if (!isPrimal(feSpace, it.getDofIndex())) {
603
604
605
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
606
607
608
	    recvFromRank = true;
	    break;
	  }
609
	}
610
      }
611
612

      if (recvFromRank)
613
	stdMpi.recv(it.getRank());
614
    }
615

616
617
    stdMpi.startCommunication();

618
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
619
	 !it.end(); it.nextRank()) {
620
      int i = 0;
621
      for (; !it.endDofIter(); it.nextDof())
622
	if (!isPrimal(feSpace, it.getDofIndex()))
623
624
625
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
626
627
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
628
629
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
630
631
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
632

633
634
635
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

636
    int nRankLagrange = 0;
637
    DofMap& dualMap = dualDofMap[feSpace].getMap();
638
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
639
640
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
641
	int degree = boundaryDofRanks[feSpace][it->first].size();
642
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
643
      } else {
644
	lagrangeMap[feSpace].insertNonRankDof(it->first);
645
646
      }
    }
647
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
648
649
650
  }


651
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
652
  {
653
    FUNCNAME("PetscSolverFeti::createIndexB()");
654

655
    DOFAdmin* admin = feSpace->getAdmin();
656
657
658
659

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

661
    int nLocalInterior = 0;    
662
    for (int i = 0; i < admin->getUsedSize(); i++) {
663
      if (admin->isDofFree(i) == false && 
664
	  isPrimal(feSpace, i) == false &&
665
666
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
667
668
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
669

670
671
672
673
674
675
676
677
678
	  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);
679
680
681

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
682
	}
683
      }
684
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
685
    
686
687
    // === And finally, add the global indicies of all dual nodes. ===

688
689
690
691
692
693
694
695
696
697
    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);
      }
698
699
700
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
705
706
    double wtime = MPI::Wtime();

707
708
    // === Create distributed matrix for Lagrange constraints. ===

709
710
711
712
713
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
714
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
715

716
717
718
719
720
721
722
    // === 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
723
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
724
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
725

726
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
727
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
728
729
730
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
731
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
732

Thomas Witkowski's avatar
Thomas Witkowski committed
733
	// Copy set of all ranks that contain this dual node.
734
735
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
736
737
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
738
739

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
740
741
742
743
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
744
745
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
746

747
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
748
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
749
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
750
	    index++;	      
751
752
753
754
755
756
757
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
759
760
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
761
762
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
763
    }
764
765
766
  }


767
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
768
  {
769
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
770

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

Thomas Witkowski's avatar
Thomas Witkowski committed
774
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
775
      
776
      VecCreateMPI(mpiCommGlobal, 
777
		   localDofMap.getRankDofs(), 
778
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
779
		   &(schurPrimalData.tmp_vec_b));
780
      VecCreateMPI(mpiCommGlobal, 
781
782
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
783
784
		   &(schurPrimalData.tmp_vec_primal));

785
      MatCreateShell(mpiCommGlobal,
786
787
788
789
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
790
791
792
793
794
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
795
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
796
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
797
798
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
799
800
      KSPSetFromOptions(ksp_schur_primal);
    } else {
801
802
      MSG("Create direct schur primal solver!\n");

803
804
      double wtime = MPI::Wtime();

805
806
807
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

808
809
810
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
811

Thomas Witkowski's avatar
Thomas Witkowski committed
812
      Mat matBPi;
813
814
815
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
816
817
818
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

819

820
821
822
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
823

824
825
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

826
      for (int i = 0; i < nRowsRankB; i++) {
827
	PetscInt row = localDofMap.getStartDofs() + i;
Thomas Witkowski's avatar
Thomas Witkowski committed
828
	MatGetRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
829
830
831
832
833

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

Thomas Witkowski's avatar
Thomas Witkowski committed
834
	MatRestoreRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
835
836
      }

837
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
838
		primalDofMap.getLocalDofs())
839
	("Should not happen!\n");
840

841
842
843
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
844
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
845
846
847
848
849
850
851
852

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

Thomas Witkowski's avatar
Thomas Witkowski committed
853
	subdomain->solve(tmpVec, tmpVec);
854

Thomas Witkowski's avatar
Thomas Witkowski committed
855
	PetscScalar *tmpValues;
856
	VecGetArray(tmpVec, &tmpValues);
857
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
858
	  MatSetValue(matBPi, 
859
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
860
861
862
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
863
864
865
866
867
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
868
869
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi<