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


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscSolverFeti.h"
16
#include "parallel/PetscSolverFetiStructs.h"
17
18
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
19
#include "parallel/PetscSolverGlobalMatrix.h"
20
#include "io/VtkWriter.h"
21
22
23
24
25

namespace AMDiS {

  using namespace std;

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


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

41
    MatMult(data->subSolver->getMatInteriorCoarse(), x, data->tmp_vec_b);
42
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
43
44
45
    MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b, 
	    data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarse(), x, y);
46
47
48
49
50
51
52
53
54
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


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

57
58
    //    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)
59

60
61
    double wtime = MPI::Wtime();

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

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

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

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

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

75
76
    MatMult(data->subSolver->getMatCoarseInterior(), 
	    data->tmp_vec_b, data->tmp_vec_primal);
77
78

    wtime01 = MPI::Wtime();
79
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
80
81
    FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);

82
83
    MatMult(data->subSolver->getMatInteriorCoarse(), 
	    data->tmp_vec_primal, data->tmp_vec_b);
84
85

    wtime01 = MPI::Wtime();
86
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
87
88
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

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

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

93
94
    FetiTimings::fetiSolve += (MPI::Wtime() - wtime);

95
96
97
98
    return 0;
  }


99
  // y = PC * x
100
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
101
  {
102
103
    double wtime = MPI::Wtime();

104
    // Get data for the preconditioner
105
106
    void *ctx;
    PCShellGetContext(pc, &ctx);
107
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
108

109
    // Multiply with scaled Lagrange constraint matrix.
110
111
112
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


113
    // === Restriction of the B nodes to the boundary nodes. ===
114

115
116
117
118
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
119

120
121
122
123
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

124
125
126
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
127
128

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


132
    // === K_DD - K_DI inv(K_II) K_ID ===
133

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

136
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
137
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
138
139
140
141
142
143
144
145
146
147
    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);

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

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

159
160
    FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);

161
162
163
164
165
166
167
168
169
170
171
172
173
174
    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);
175
176


177
    // === Restriction of the B nodes to the boundary nodes. ===
178

179
180
181
182
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
183

184
    PetscScalar *local_b, *local_duals;
185
    VecGetArray(data->tmp_vec_b, &local_b);
186
    VecGetArray(data->tmp_vec_duals0, &local_duals);
187

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

    VecRestoreArray(data->tmp_vec_b, &local_b);
193
194
195
196
197
198
199
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

200

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

203
204
205
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

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

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

213
214

    // Multiply with scaled Lagrange constraint matrix.
215
216
217
218
219
220
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


221
222
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
223
      schurPrimalSolver(0),
224
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
225
      subdomain(NULL),
226
227
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
228
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
229
230
      printTimings(false),
      enableStokesMode(false)
231
232
233
234
235
236
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

    Parameters::get("parallel->multi level test", multiLevelTest);
256
257
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
258
259

    Parameters::get("parallel->print timings", printTimings);
260
261
262
  }


263
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
264
  {
265
266
    FUNCNAME("PetscSolverFeti::initialize()");

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

270
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
271
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
272

Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
275
276
277
278
279
280
281

    enableStokesMode = false;
    Parameters::get("parallel->stokes mode", enableStokesMode);
    if (enableStokesMode) {
      MSG("========= !!!! SUPER WARNING !!! ========\n");
      MSG("Use make use of stokes mode which is work in progress! We guarantee for nothing and even less than this!\n");
      fullInterface = feSpaces[2];
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
282
283
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
284

285
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
286
	subdomain->setMeshDistributor(meshDistributor, 
287
				      mpiCommGlobal, mpiCommLocal);
288
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
289
	subdomain->setMeshDistributor(meshDistributor, 
290
291
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
292
	subdomain->setLevel(meshLevel);
293
294
      }
    }
295

296
297
298
299
    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
300

301
    if (fullInterface != NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
302
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
303

304
    if (fetiPreconditioner != FETI_NONE) {
305
306
307
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

308
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
309
    }
310
311
312
313
314
315
316
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
317
318
    double timeCounter = MPI::Wtime();

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

323
324
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

325
326
327
328
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
329
    if (fetiPreconditioner != FETI_NONE)
330
331
      interiorDofMap.clear();

332
333
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
334

335
336
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
337
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
338
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
339
340
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
341

342
343
344
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
345
346
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

347
348
349
350
351
352
    if (fullInterface != NULL) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

353
354
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
355
    
356
357
      createPrimals(feSpace);  

358
359
360
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
361

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
362
      createIndexB(feSpace);     
363
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
366
367

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
368
    if (fetiPreconditioner != FETI_NONE)
369
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
370

371
372
373
    if (fullInterface != NULL)
      interfaceDofMap.update();

374
375
376
377
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
378

379
380
381
    lagrangeMap.update();


382
383
384
385
386
387
388
389
390
391
392
393
    // === ===

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

394
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
395
396
397
398
399
400
401
402
403
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
      if (feSpace == fullInterface) {
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
	MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
	    primalDofMap[feSpace].nRankDofs, 
	    primalDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
	    dualDofMap[feSpace].nRankDofs, 
	    dualDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
	    lagrangeMap[feSpace].nRankDofs, 
	    lagrangeMap[feSpace].nOverallDofs);

	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
430
    }
431

Thomas Witkowski's avatar
Thomas Witkowski committed
432
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
433
434
435
436
437
438
439
440
441

    if (enableStokesMode) {
      MSG("WARNING: MAKE THIS MORE GENERAL!!!!!\n");
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 0);
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 1);
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, 2);
    } else {	  
      subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
442
443

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
444
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
445
      timeCounter = MPI::Wtime() - timeCounter;
446
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
	  timeCounter);
    }
449
450
451
  }


452
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
453
  {
454
    FUNCNAME("PetscSolverFeti::createPrimals()");  
455

456
457
458
    if (feSpace == fullInterface)
      return;

459
460
461
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
462
    // Set of DOF indices that are considered to be primal variables.
463
    DofContainerSet& vertices = 
464
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
465
466

    DofIndexSet primals;
467
    for (DofContainerSet::iterator it = vertices.begin(); 
468
469
	 it != vertices.end(); ++it) {
      double e = 1e-8;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
470
      if (meshLevel == 0) {
471
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
472
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
473
474
475
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
476
477
478
479
480
481
482
	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);
	}
      }
483
    }
484

485
486
487
488

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

489
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
490
491
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
492
      else
493
  	primalDofMap[feSpace].insertNonRankDof(*it);
494
495
496
  }


497
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
498
499
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
500

501
502
503
    if (feSpace == fullInterface)
      return;

504
505
506
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
507
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
508
509
510
511

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
512
513
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
514
515
516
517
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
518
519
520
  }

  
521
522
523
524
525
526
527
528
529
530
531
532
  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)
Thomas Witkowski's avatar
Thomas Witkowski committed
533
534
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it)) {
	MSG("INSERT INTERFACE RANK DOF %d\n", **it);
Thomas Witkowski's avatar
Thomas Witkowski committed
535
	interfaceDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Thomas Witkowski committed
536
537
      } else {
	MSG("INSERT INTERFACE NON-RANK DOF %d\n", **it);
Thomas Witkowski's avatar
Thomas Witkowski committed
538
	interfaceDofMap[feSpace].insertNonRankDof(**it);
Thomas Witkowski's avatar
Thomas Witkowski committed
539
      }
540
541
542
  }


543
544
545
546
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

547
548
549
    if (feSpace == fullInterface)
      return;

550
551
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
552
553
554
    // 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
555
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
556

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
557
    if (meshLevel > 0) {
558
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
559
560
561
562
563
564
565
566
567

      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
568
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
569
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
570
	  else
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
	    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
586
587
588
589
590
591
	   !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());
    }
592

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
593
594
595
596
597
598
599
    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)).         ===

600
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
601
602
603
604
605
606
607
    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);

608
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
609
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
610
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
611
612
613
614
	}
      }
    }

615
616
617
618

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

619
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
620

621
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
622
623
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
624
	if (!isPrimal(feSpace, it.getDofIndex()))
625
626
627
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
628
629
630

    stdMpi.updateSendDataSize();

631
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
632
	 !it.end(); it.nextRank()) {
633
      bool recvFromRank = false;
634
      for (; !it.endDofIter(); it.nextDof()) {
635
	if (!isPrimal(feSpace, it.getDofIndex())) {
636
637
638
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
639
640
641
	    recvFromRank = true;
	    break;
	  }
642
	}
643
      }
644
645

      if (recvFromRank)
646
	stdMpi.recv(it.getRank());
647
    }
648

649
650
    stdMpi.startCommunication();

651
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
652
	 !it.end(); it.nextRank()) {
653
      int i = 0;
654
      for (; !it.endDofIter(); it.nextDof())
655
	if (!isPrimal(feSpace, it.getDofIndex()))
656
657
658
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
659
660
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
661
662
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
663
664
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
665

666
667
668
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

669
    int nRankLagrange = 0;
670
    DofMap& dualMap = dualDofMap[feSpace].getMap();
671
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
672
673
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
674
	int degree = boundaryDofRanks[feSpace][it->first].size();
675
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
676
      } else {
677
	lagrangeMap[feSpace].insertNonRankDof(it->first);
678
679
      }
    }
680
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
681
682
683
  }


684
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
685
  {
686
    FUNCNAME("PetscSolverFeti::createIndexB()");
687

688
    DOFAdmin* admin = feSpace->getAdmin();
689
690
691
692

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

694
    int nLocalInterior = 0;    
695
    for (int i = 0; i < admin->getUsedSize(); i++) {
696
      if (admin->isDofFree(i) == false && 
697
	  isPrimal(feSpace, i) == false &&
698
699
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
700
701
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
702

703
704
705
706
707
708
709
710
711
	  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);
712
713
714

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
715
	}
716
      }
717
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
718
    
719
720
    // === And finally, add the global indicies of all dual nodes. ===

721
722
723
724
725
726
727
728
729
730
    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);
      }
731
732
733
  }


734
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
735
736
737
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
738
    double wtime = MPI::Wtime();
739
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
740

741
742
    // === Create distributed matrix for Lagrange constraints. ===

743
744
745
746
747
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
748
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
749

750
751
752
753
754
755
756
    // === 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
757
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
758
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
759

760
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
761
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
762
763
764
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
765
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
766

Thomas Witkowski's avatar
Thomas Witkowski committed
767
	// Copy set of all ranks that contain this dual node.
768
769
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
770
771
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
772
773

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
774
775
776
777
	
	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
778
779
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
780

781
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
782
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
783
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
784
	    index++;	      
785
786
787
788
789
790
791
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
793
794
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
795
796
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
797
    }
798
799
800
  }


801
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
802
  {
803
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
804

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

Thomas Witkowski's avatar
Thomas Witkowski committed
808
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
809
      
810
      VecCreateMPI(mpiCommGlobal, 
811
		   localDofMap.getRankDofs(), 
812
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
813
		   &(schurPrimalData.tmp_vec_b));
814
      VecCreateMPI(mpiCommGlobal, 
815
816
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
817
818
		   &(schurPrimalData.tmp_vec_primal));

819
      MatCreateShell(mpiCommGlobal,
820
821
822
823
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
824
825
826
827
828
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
829
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
830
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
831
832
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
833
834
      KSPSetFromOptions(ksp_schur_primal);
    } else {
835
836
      MSG("Create direct schur primal solver!\n");

837
838
      double wtime = MPI::Wtime();

839
840
841
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

842
843
844
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
845

Thomas Witkowski's avatar
Thomas Witkowski committed
846
      Mat matBPi;
847
848
849
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
850
851
852
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

853

854
855
856
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
857

858
859
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

860
      for (int i = 0; i < nRowsRankB; i++) {
861
	PetscInt row = localDofMap.getStartDofs() + i;
862
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
863
864
865
866
867

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

868
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
869
870
      }

871
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
872
		primalDofMap.getLocalDofs())
873
	("Should not happen!\n");
874

875
876
877
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
878
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
879
880
881
882
883
884
885
886

 	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
887
	subdomain->solve(tmpVec, tmpVec);
888

Thomas Witkowski's avatar
Thomas Witkowski committed
889
	PetscScalar *tmpValues;
890
	VecGetArray(tmpVec, &tmpValues);
891
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
892
	  MatSetValue(matBPi, 
893
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
894
895
896
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
897
898
899
900
901
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
902
903
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
904

905
906
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
907
908

      Mat matPrimal;
909
910
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
911
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
912
913

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
914
      MatDestroy(&matBPi);
915

916
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
917
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
918
		      SAME_NONZERO_PATTERN);
919
920
921
922
923
924
      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);
925
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
926

Thomas Witkowski's avatar
Thomas Witkowski committed
927
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
928
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
929
930
931
932
933
934
935
936
937
938
	MatInfo minfo;
	MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
	MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);
	
	MSG("FETI-DP timing 06: %.5f seconds (creation of schur primal matrix)\n",
	    MPI::Wtime() - wtime);

	wtime = MPI::Wtime();
	KSPSetUp(ksp_schur_primal);
	KSPSetUpOnBlocks(ksp_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
939
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
940
941
942
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }