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

    return 0;
  }


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

92
93
94
95
    return 0;
  }


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

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

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


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

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

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

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

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


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

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

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

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

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

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

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


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

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

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

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

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


    // === K_DD ===

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

197

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

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

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

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

210
211

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

    return 0;
  }


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

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

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

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

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


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

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

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

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

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

283
284
285
286
    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
287

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

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

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


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

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

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

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

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

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

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

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

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

340
341
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
342
    
343
344
      createPrimals(feSpace);  

345
346
347
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
348

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

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

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

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

366
367
368
    lagrangeMap.update();


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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
    subdomain->setDofMapping(&localDofMap);
    subdomain->setCoarseSpaceDofMapping(&primalDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
415
416

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


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

429
430
431
    if (feSpace == fullInterface)
      return;

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

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

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

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

458
459
460
461

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

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


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

474
475
476
    if (feSpace == fullInterface)
      return;

477
478
479
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
480
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
481
482
483
484

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

  
494
495
496
497
498
499
500
501
502
503
504
505
  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
506
507
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
	interfaceDofMap[feSpace].insertRankDof(**it);
508
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
509
	interfaceDofMap[feSpace].insertNonRankDof(**it);
510
511
512
  }


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

517
518
519
    if (feSpace == fullInterface)
      return;

520
521
    boundaryDofRanks[feSpace].clear();

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

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

      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
538
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
539
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
540
	  else
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
	    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
556
557
558
559
560
561
	   !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());
    }
562

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
563
564
565
566
567
568
569
570
571
572
573
574
575
576
    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);

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

584
585
586
587

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

588
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
589

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
615
	stdMpi.recv(it.getRank());
616
    }
617

618
619
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
634

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

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


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

657
    DOFAdmin* admin = feSpace->getAdmin();
658
659
660
661

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

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

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

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

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


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

Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
    double wtime = MPI::Wtime();

709
710
    // === Create distributed matrix for Lagrange constraints. ===

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

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

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
742
743
744
745
	
	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
746
747
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
748

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

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

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


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

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

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

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

805
806
      double wtime = MPI::Wtime();

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

810
811
812
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
813

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

821

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

826
827
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

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

	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
836
	MatRestoreRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
837
838
      }

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

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

 	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
855
	subdomain->solve(tmpVec, tmpVec);
856

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
870
871
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
872

Thomas Witkowski's avatar
Thomas Witkowski committed
873
      MatDuplicate(subdomain->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
874
875

      Mat matPrimal;
Thomas Witkowski's avatar
Thomas Witkowski committed
876
      MatMatMult(subdomain->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
877
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
878
879

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
880
      MatDestroy(&matBPi);
881

882
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
883
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
884
		      SAME_NONZERO_PATTERN);
885
886
887
888
889
890
      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);
891
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
892

Thomas Witkowski's avatar
Thomas Witkowski committed
893
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
894
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
895
896
897
898
899
900
901
902
903
904
	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
905
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
906
907
908
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
909
    }
910
911
912
913
914
915
916
  }


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

917
    if (schurPrimalSolver == 0) {
918
      schurPrimalData.subSolver = NULL;
919

920
921
922
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
923
924
925