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


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/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
303
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
      //      interfaceDofMap.init(levelData, fullInterface);
Thomas Witkowski's avatar
Thomas Witkowski committed
304

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

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


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

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

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

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

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

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

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

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

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

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

359
360
361
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
362

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

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

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

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

380
381
382
    lagrangeMap.update();


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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
      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
431
    }
432

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

    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
443
444

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


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

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

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

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

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

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

486
487
488
489

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

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


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

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

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

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

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

  
522
523
524
525
526
527
528
529
530
531
532
533
  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
534
535
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
	interfaceDofMap[feSpace].insertRankDof(**it);
536
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
537
	interfaceDofMap[feSpace].insertNonRankDof(**it);
538
539
540
  }


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

545
546
547
    if (feSpace == fullInterface)
      return;

548
549
    boundaryDofRanks[feSpace].clear();

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

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

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

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

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

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

613
614
615
616

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

617
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
618

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
644
	stdMpi.recv(it.getRank());
645
    }
646

647
648
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
663

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

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


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

686
    DOFAdmin* admin = feSpace->getAdmin();
687
688
689
690

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

692
    int nLocalInterior = 0;    
693
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
696
697
      if (i == 0) {
	MSG("JETZT: %d %d %d %d\n", admin->isDofFree(i), isPrimal(feSpace, i), isDual(feSpace, i), isInterface(feSpace, i));
      }

698
      if (admin->isDofFree(i) == false && 
699
	  isPrimal(feSpace, i) == false &&
700
701
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
702
703
704
705
	if (i == 0) {
	  MSG("DRIN\n");
	}

706
707
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
708

709
710
711
712
713
714
715
716
717
	  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);
718
719
720

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
721
	}
722
      }
723
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
724
    
725
726
    // === And finally, add the global indicies of all dual nodes. ===

727
728
729
730
731
732
733
734
735
736
    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);
      }
737
738
739
  }


740
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
741
742
743
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
744
    double wtime = MPI::Wtime();
745
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
746

747
748
    // === Create distributed matrix for Lagrange constraints. ===

749
750
751
752
753
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
754
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
755

756
757
758
759
760
761
762
    // === 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
763
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
764
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
765

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

Thomas Witkowski's avatar
Thomas Witkowski committed
773
	// Copy set of all ranks that contain this dual node.
774
775
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
776
777
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
778
779

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
780
781
782
783
	
	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
784
785
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
786

787
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
788
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
789
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
790
	    index++;	      
791
792
793
794
795
796
797
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
799
800
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
801
802
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
803
    }
804
805
806
  }


807
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
808
  {
809
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
810

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

Thomas Witkowski's avatar
Thomas Witkowski committed
814
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
815
      
816
      VecCreateMPI(mpiCommGlobal, 
817
		   localDofMap.getRankDofs(), 
818
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
819
		   &(schurPrimalData.tmp_vec_b));
820
      VecCreateMPI(mpiCommGlobal, 
821
822
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
		   &(schurPrimalData.tmp_vec_primal));

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

843
844
      double wtime = MPI::Wtime();

845
846
847
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

848
849
850
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
851

Thomas Witkowski's avatar
Thomas Witkowski committed
852
      Mat matBPi;
853
854
855
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
856
857
858
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

859

860
861
862
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
863

864
865
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

866
      for (int i = 0; i < nRowsRankB; i++) {
867
	PetscInt row = localDofMap.getStartDofs() + i;
868
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
869
870
871
872
873

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

874
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
875
876
      }

877
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
878
		primalDofMap.getLocalDofs())
879
	("Should not happen!\n");
880

881
882
883
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;