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->getDofComm());
    else
352
353
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

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

357
358
      createPrimals(feSpace);  

359
      createDuals(feSpace);      
360

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

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

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

375
376
377
    lagrangeMap.update();


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

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

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

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

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

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

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

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

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


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

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

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

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

450
      double e = 1e-8;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
451
      if (meshLevel == 0) {
452
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
453
454
455
456
457
458
459
460
      } 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);
	}
      }
461
    }
462

463
464
465
466

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

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


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

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

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

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

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

500
501
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
502
503
504
505
    // 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;
506

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

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

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

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

562
563
564
565
	}
      }
    }

566
567
568
569

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

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

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

    stdMpi.updateSendDataSize();

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

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

600
601
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
616

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

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


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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
800
      Mat matPrimal;
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
854
      MatDuplicate(subdomain->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
      MatMatMult(subdomain->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
855
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
856
857

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
858
      MatDestroy(&matBPi);
859

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

Thomas Witkowski's avatar
Thomas Witkowski committed
871
872
873
874
875
876
877
878
879
880
881
882
883
884
      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
885
    }
886
887
888
889
890
891
892
  }


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

893
    if (schurPrimalSolver == 0) {
894
      schurPrimalData.subSolver = NULL;
895

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


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

909
910
    // === Create FETI-DP solver object. ===

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

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