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


13
#include "AMDiS.h"
14
#include "parallel/PetscSolverFeti.h"
15
#include "parallel/PetscSolverFetiStructs.h"
16
17
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
18
#include "parallel/PetscSolverGlobalMatrix.h"
19
#include "io/VtkWriter.h"
20
21
22
23
24

namespace AMDiS {

  using namespace std;

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


31
32
33
34
35
36
37
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
38
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
39

40
41
42
43
    MatMult(data->subSolver->getMatIntCoarse(), x, data->tmp_vec_b);
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarseCoarse(), x, y);
44
45
46
47
48
49
50
51
52
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


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

55
56
    //    F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
57

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

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

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

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

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

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

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

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

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

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

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

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

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

91
92
93
94
    return 0;
  }


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

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

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


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

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

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

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

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


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

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

132
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
133
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
134
135
136
137
138
139
140
141
142
143
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

    VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);


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

    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

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

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


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

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

157
158
159
160
161
162
163
164
165
166
167
168
169
170
    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
171
172


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

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

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

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

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


    // === K_DD ===

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

196

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

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

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

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

209
210

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

    return 0;
  }


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

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

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

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

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


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

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

265
266
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

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

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

281
282
    primalDofMap.init(levelData,
		      feSpaces, meshDistributor->getFeSpaces(), 
283
		      true, true);
284
285
    primalDofMap.setComputeMatIndex(true);
    
286
287
    dualDofMap.init(levelData, 
		    feSpaces, meshDistributor->getFeSpaces(),
288
		    false, false);
289
290
    dualDofMap.setComputeMatIndex(true);

291
292
293
294
295
296
297
298
299
300
301
    if (meshLevel == 0) {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       false, false);
      localDofMap.setComputeMatIndex(true);
    } else {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       true, true);
      localDofMap.setComputeMatIndex(true);
    }
302

Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
305
306
307
308
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
		     true, true);
    lagrangeMap.setComputeMatIndex(true);


309
    if (fetiPreconditioner != FETI_NONE) {
310
311
312
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

313
314
      interiorDofMap.init(levelData, 
			  feSpaces,  meshDistributor->getFeSpaces(),
315
			  false, false);
316
317
      interiorDofMap.setComputeMatIndex(true);
    }
318
319
320
321
322
323
324
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
    double timeCounter = MPI::Wtime();

327
328
329
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
330

331
332
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

333
334
335
336
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
337
    if (fetiPreconditioner != FETI_NONE)
338
339
      interiorDofMap.clear();

340
341
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
342

343
344
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
345
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
346
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
347
348
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
349

350
351
352
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
353
354
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

355
356
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
357

358
359
      createPrimals(feSpace);  

360
      createDuals(feSpace);      
361

362
363
      createIndexB(feSpace);
    }
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
374
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
375

376
377
378
    lagrangeMap.update();


379
380
381
382
383
384
385
386
387
388
389
390
    // === ===

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

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

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

406
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
407
408
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
409
      
410
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
411
412
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
413

414
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
415
416
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
417

418
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
419
420
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
421
    }
422

Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
425
    subdomain->setDofMapping(&localDofMap, &primalDofMap);

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
426
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
427
      timeCounter = MPI::Wtime() - timeCounter;
428
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
429
430
	  timeCounter);
    }
431
432
433
  }


434
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
435
  {
436
    FUNCNAME("PetscSolverFeti::createPrimals()");  
437

438
439
440
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
441
    // Set of DOF indices that are considered to be primal variables.
442
    DofContainerSet& vertices = 
443
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
444
445

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
455
456
457
458
459
460
461
	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);
	}
      }
462
    }
463

464
465
466
467

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

468
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
469
470
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
471
      else
472
  	primalDofMap[feSpace].insertNonRankDof(*it);
473
474
475
  }


476
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
477
478
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
479

480
481
482
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
483
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
484
485
486
487

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

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

501
502
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
503
504
505
    // 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
506
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
507

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
508
    if (meshLevel > 0) {
509
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
510
511
512
513
514
515
516
517
518

      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
519
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
520
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
521
	  else
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
	    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
537
538
539
540
541
542
	   !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());
    }
543

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
544
545
546
547
548
549
550
551
552
553
554
555
556
557
    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);

558
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
559
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
560
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
561
562
563
564
	}
      }
    }

565
566
567
568

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

569
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
570

571
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
572
573
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
574
	if (!isPrimal(feSpace, it.getDofIndex()))
575
576
577
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
578
579
580

    stdMpi.updateSendDataSize();

581
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
582
	 !it.end(); it.nextRank()) {
583
      bool recvFromRank = false;
584
      for (; !it.endDofIter(); it.nextDof()) {
585
	if (!isPrimal(feSpace, it.getDofIndex())) {
586
587
588
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
589
590
591
	    recvFromRank = true;
	    break;
	  }
592
	}
593
      }
594
595

      if (recvFromRank)
596
	stdMpi.recv(it.getRank());
597
    }
598

599
600
    stdMpi.startCommunication();

601
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
602
	 !it.end(); it.nextRank()) {
603
      int i = 0;
604
      for (; !it.endDofIter(); it.nextDof())
605
	if (!isPrimal(feSpace, it.getDofIndex()))
606
607
608
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
609
610
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
611
612
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
613
614
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
615

616
617
618
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

619
    int nRankLagrange = 0;
620
    DofMap& dualMap = dualDofMap[feSpace].getMap();
621
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
622
623
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
624
	int degree = boundaryDofRanks[feSpace][it->first].size();
625
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
626
      } else {
627
	lagrangeMap[feSpace].insertNonRankDof(it->first);
628
629
      }
    }
630
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
631
632
633
  }


634
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
635
  {
636
    FUNCNAME("PetscSolverFeti::createIndexB()");
637

638
    DOFAdmin* admin = feSpace->getAdmin();
639
640
641
642

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

644
    int nLocalInterior = 0;    
645
    for (int i = 0; i < admin->getUsedSize(); i++) {
646
      if (admin->isDofFree(i) == false && 
647
	  isPrimal(feSpace, i) == false &&
648
	  isDual(feSpace, i) == false) {
649
650
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
651

652
653
654
655
656
657
658
659
660
	  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);
661
662
663

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
664
	}
665
      }
666
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
667
    
668
669
    // === And finally, add the global indicies of all dual nodes. ===

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


683
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
684
685
686
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
687
688
    double wtime = MPI::Wtime();

689
690
    // === Create distributed matrix for Lagrange constraints. ===

691
692
693
694
695
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
696
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
697

698
699
700
701
702
703
704
    // === 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
705
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
706
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
707

708
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
709
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
712
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
713
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
714

Thomas Witkowski's avatar
Thomas Witkowski committed
715
	// Copy set of all ranks that contain this dual node.
716
717
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
720
721

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
722
723
724
725
	
	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
726
727
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
728

729
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
730
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
731
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
732
	    index++;	      
733
734
735
736
737
738
739
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
741
742
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
745
    }
746
747
748
  }


749
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
750
  {
751
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
752

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

Thomas Witkowski's avatar
Thomas Witkowski committed
756
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
757
      
758
      VecCreateMPI(mpiCommGlobal, 
759
		   localDofMap.getRankDofs(), 
760
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
761
		   &(schurPrimalData.tmp_vec_b));
762
      VecCreateMPI(mpiCommGlobal, 
763
764
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
765
766
		   &(schurPrimalData.tmp_vec_primal));

767
      MatCreateShell(mpiCommGlobal,
768
769
770
771
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
772
773
774
775
776
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
777
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
778
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
779
780
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
781
782
      KSPSetFromOptions(ksp_schur_primal);
    } else {
783
784
      MSG("Create direct schur primal solver!\n");

785
786
      double wtime = MPI::Wtime();

787
788
789
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

790
791
792
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
793

Thomas Witkowski's avatar
Thomas Witkowski committed
794
      Mat matBPi;
795
796
797
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
798
799
800
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

801

802
803
804
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
805

806
807
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

808
      for (int i = 0; i < nRowsRankB; i++) {
809
	PetscInt row = localDofMap.getStartDofs() + i;
Thomas Witkowski's avatar
Thomas Witkowski committed
810
	MatGetRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
811
812
813
814
815

	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
816
	MatRestoreRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
817
818
      }

819
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
820
		primalDofMap.getLocalDofs())
821
	("Should not happen!\n");
822

823
824
825
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
826
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
827
828
829
830
831
832
833
834

 	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
835
	subdomain->solve(tmpVec, tmpVec);
836

Thomas Witkowski's avatar
Thomas Witkowski committed
837
	PetscScalar *tmpValues;
838
	VecGetArray(tmpVec, &tmpValues);
839
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
840
	  MatSetValue(matBPi, 
841
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
842
843
844
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
845
846
847
848
849
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
850
851
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
852

Thomas Witkowski's avatar
Thomas Witkowski committed
853
      MatDuplicate(subdomain->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
854
855

      Mat matPrimal;
Thomas Witkowski's avatar
Thomas Witkowski committed
856
      MatMatMult(subdomain->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
857
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
858
859

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
860
      MatDestroy(&matBPi);
861

862
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
863
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
864
		      SAME_NONZERO_PATTERN);
865
866
867
868
869
870
      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);
871
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
872

Thomas Witkowski's avatar
Thomas Witkowski committed
873
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
874
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
875
876
877
878
879
880
881
882
883
884
	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
885
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
886
887
888
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
889
    }
890
891
892
893
894
895
896
  }


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

897
    if (schurPrimalSolver == 0) {
898
      schurPrimalData.subSolver = NULL;
899

900
901
902
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
903
904
905
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
906
907
908
  }


909
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
910
911
912
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

913
914
    // === Create FETI-DP solver object. ===

915
    fetiData.mat_lagrange = &mat_lagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
916
    fetiData.subSolver = subdomain;
917
    fetiData.ksp_schur_primal = &ksp_schur_primal;
918

919
    VecCreateMPI(mpiCommGlobal, 
920
		 localDofMap.getRankDofs(),