PetscSolverFeti.cc 50.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"
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
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
293
		     true, true);
294
    lagrangeMap.setComputeMatIndex(true);
295
296
297
298
299
300
301
302
303
304
305
306

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

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

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


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

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

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

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

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

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

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

    if (meshLevel > 0)
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

352
353
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
354

355
356
      createPrimals(feSpace);  

357
      createDuals(feSpace);      
358

359
360
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
363
364

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
365
    if (fetiPreconditioner != FETI_NONE)
366
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
367

368
369
370
371
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
372

373
374
375
    lagrangeMap.update();


376
377
378
379
380
381
382
383
384
385
386
387
    // === ===

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

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

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

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

411
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
412
413
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
414

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

Thomas Witkowski's avatar
Thomas Witkowski committed
420
421
422
423
    subdomain->setDofMapping(&localDofMap, &primalDofMap);

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


430
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
431
  {
432
    FUNCNAME("PetscSolverFeti::createPrimals()");  
433

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

437
438
    /// Set of DOF indices that are considered to be primal variables.
    
439
    DofContainerSet& vertices = 
440
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
441
442

    DofIndexSet primals;
443
    for (DofContainerSet::iterator it = vertices.begin(); 
444
445
446
	 it != vertices.end(); ++it) {
      WorldVector<double> c;
      feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
447

448
      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
452
453
454
455
456
457
458
      } else {
	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);
	}
      }
459
    }
460

461
462
463
464

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

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


473
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
474
475
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
476

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
  }

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

498
499
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
500
501
502
503
    // Stores for all rank own communication DOFs, if the counterpart is
    // a rank owmed DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFS that fulfill this requirenment.
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
504

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

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

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

555
556
557
558
559
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex()))) {
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
 	  }

560
561
562
563
	}
      }
    }

564
565
566
567

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

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

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

    stdMpi.updateSendDataSize();

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

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

598
599
    stdMpi.startCommunication();

600
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
601
	 !it.end(); it.nextRank()) {
602
      int i = 0;
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
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
610
611
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
612
613
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
614

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

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


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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

    if (printTimings)
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
743
744
745
  }


746
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
747
  {
748
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
749

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

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

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

782
783
      double wtime = MPI::Wtime();

784
785
786
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

787
788
789
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
790

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

Thomas Witkowski's avatar
Thomas Witkowski committed
798
      Mat matPrimal;
799

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

804
805
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

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

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

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

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

 	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
833
	subdomain->solve(tmpVec, tmpVec);
834

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
848
849
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
850

Thomas Witkowski's avatar
Thomas Witkowski committed
851
852
      MatDuplicate(subdomain->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
      MatMatMult(subdomain->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
853
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
854
855

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
856
      MatDestroy(&matBPi);
857

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

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


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

891
    if (schurPrimalSolver == 0) {
892
      schurPrimalData.subSolver = NULL;
893

894
895
896
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
897
898
899
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
900
901
902
  }


903
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
904
905
906
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

907
908
    // === Create FETI-DP solver object. ===

909
    fetiData.mat_lagrange = &mat_lagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
910
    fetiData.subSolver = subdomain;
911
    fetiData.ksp_schur_primal = &ksp_schur_primal;
912

913
    VecCreateMPI(mpiCommGlobal, 
914
		 localDofMap.getRankDofs(), 
915
		 nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
916
		 &(fetiData.tmp_vec_b));
917
    VecCreateMPI(mpiCommGlobal,
918
919
		 lagrangeMap.getRankDofs(), 
		 lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
920
		 &(fetiData.tmp_vec_lagrange));
921
    VecCreateMPI(mpiCommGlobal, 
922
923
		 primalDofMap.getRankDofs(), 
		 primalDofMap.getOverallDofs(),
924
		 &(fetiData.tmp_vec_primal));
925

926
    MatCreateShell(mpiCommGlobal,
927
928
929
930
		   lagrangeMap.getRankDofs(), 
		   lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(),