PetscSolverFeti.cc 48.9 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
31
  // 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);
32
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
33

34
35
36
37
    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);
38
39
40
41
42
43
44
45
46
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
47
48
    //    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)
49
50
51

    void *ctx;
    MatShellGetContext(mat, &ctx);
52
    FetiData* data = static_cast<FetiData*>(ctx);
53

Thomas Witkowski's avatar
Thomas Witkowski committed
54
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
55
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
57

58
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
59
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
60
    MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
61
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
63

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
65
66
67
68
69

    return 0;
  }


70
  // y = PC * x
71
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
72
  {
73
    // Get data for the preconditioner
74
75
    void *ctx;
    PCShellGetContext(pc, &ctx);
76
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
77

78
    // Multiply with scaled Lagrange constraint matrix.
79
80
81
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


82
    // === Restriction of the B nodes to the boundary nodes. ===
83

84
85
86
87
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
88

89
90
91
92
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

93
94
95
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
96
97

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


101
    // === K_DD - K_DI inv(K_II) K_ID ===
102

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

105
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
106
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
107
108
109
110
111
112
113
114
115
116
    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);

117
118
119
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

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

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


144
    // === Restriction of the B nodes to the boundary nodes. ===
145

146
147
148
149
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
150

151
    PetscScalar *local_b, *local_duals;
152
    VecGetArray(data->tmp_vec_b, &local_b);
153
    VecGetArray(data->tmp_vec_duals0, &local_duals);
154

Thomas Witkowski's avatar
Thomas Witkowski committed
155
156
157
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
158
159

    VecRestoreArray(data->tmp_vec_b, &local_b);
160
161
162
163
164
165
166
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

167

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

170
171
172
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

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

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

180
181

    // Multiply with scaled Lagrange constraint matrix.
182
183
184
185
186
187
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


188
189
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
190
      schurPrimalSolver(0),
191
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
192
      subdomain(NULL),
193
194
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
      nGlobalOverallInterior(0),
      printTimings(false)
197
198
199
200
201
202
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
203
      MSG("Create FETI-DP solver with no preconditioner!\n");
204
205
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
206
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
207
208
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
209
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
210
211
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
212
213
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
214
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
215
216
217
218
219

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

    Parameters::get("parallel->multi level test", multiLevelTest);
222
223
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
224
225

    Parameters::get("parallel->print timings", printTimings);
226
227
228
  }


229
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
230
  {
231
232
    FUNCNAME("PetscSolverFeti::initialize()");

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

236
237
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

Thomas Witkowski's avatar
Thomas Witkowski committed
238
239
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
240

241
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
242
	subdomain->setMeshDistributor(meshDistributor, 
243
					    mpiCommGlobal, mpiCommLocal);
244
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
245
	subdomain->setMeshDistributor(meshDistributor, 
246
247
					    levelData.getMpiComm(meshLevel - 1),
					    levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
248
	subdomain->setLevel(meshLevel);
249
250
      }
    }
251

252
253
    primalDofMap.init(levelData,
		      feSpaces, meshDistributor->getFeSpaces(), 
254
		      true, true);
255
256
    primalDofMap.setComputeMatIndex(true);
    
257
258
    dualDofMap.init(levelData, 
		    feSpaces, meshDistributor->getFeSpaces(),
259
		    false, false);
260
261
    dualDofMap.setComputeMatIndex(true);

262
263
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
264
		     true, true);
265
    lagrangeMap.setComputeMatIndex(true);
266
267
268
269
270
271
272
273
274
275
276
277

    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);
    }
278
279

    if (fetiPreconditioner == FETI_DIRICHLET) {
280
281
282
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

283
284
      interiorDofMap.init(levelData, 
			  feSpaces,  meshDistributor->getFeSpaces(),
285
			  false, false);
286
287
      interiorDofMap.setComputeMatIndex(true);
    }
288
289
290
291
292
293
294
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
295
296
    double timeCounter = MPI::Wtime();

297
298
299
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
300

301
302
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

303
304
305
306
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
307
308
309
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.clear();

310
311
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
312

313
314
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
315
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
316
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
317
318
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
319
320
321
322

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

323
324
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
325

326
327
      createPrimals(feSpace);  

328
      createDuals(feSpace);      
329

330
331
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
334
335

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
336
    if (fetiPreconditioner != FETI_NONE)
337
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
338

339
340
341
342
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
343

344
345
346
    lagrangeMap.update();


347
348
349
350
351
352
353
354
355
356
357
358
    // === ===

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

359
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
360
361
362
363
364
365
366
367
368
			   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);
    }

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

374
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
375
376
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
377
      
378
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
379
380
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
381

382
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
383
384
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
385

386
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
387
388
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
389
    }
390
391

    // If multi level test, inform sub domain solver about coarse space.
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
395
396
397
398
    subdomain->setDofMapping(&localDofMap, &primalDofMap);

    if (printTimings) {
      timeCounter = MPI::Wtime() - timeCounter;
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)", 
	  timeCounter);
    }
399
400
401
  }


402
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
403
  {
404
    FUNCNAME("PetscSolverFeti::createPrimals()");  
405

406
407
408
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

409
410
    /// Set of DOF indices that are considered to be primal variables.
    
411
    DofContainerSet& vertices = 
412
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
413
414

    DofIndexSet primals;
415
    for (DofContainerSet::iterator it = vertices.begin(); 
416
417
418
	 it != vertices.end(); ++it) {
      WorldVector<double> c;
      feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
419

420
      double e = 1e-8;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
421
      if (meshLevel == 0) {
422
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
423
424
425
426
427
428
429
430
      } 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);
	}
      }
431
    }
432

433
434
435
436

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

437
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
438
439
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
440
      else
441
  	primalDofMap[feSpace].insertNonRankDof(*it);
442
443
444
  }


445
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
446
447
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
448

449
450
451
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
452
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
453
454
455
456

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
457
458
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
459
460
461
462
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
463
464
465
466
467
468
469
  }

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

470
471
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
472
473
474
475
    // 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;
476

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
477
    if (meshLevel > 0) {
478
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
479
480
481
482
483
484
485
486
487

      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
488
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
489
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
490
	  else
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
	    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
506
507
508
509
510
511
	   !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());
    }
512

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
513
514
515
516
517
518
519
520
521
522
523
524
525
526
    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);

527
528
529
530
531
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex()))) {
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
 	  }

532
533
534
535
	}
      }
    }

536
537
538
539

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

540
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
541

542
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
543
544
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
545
	if (!isPrimal(feSpace, it.getDofIndex()))
546
547
548
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
549
550
551

    stdMpi.updateSendDataSize();

552
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
553
	 !it.end(); it.nextRank()) {
554
      bool recvFromRank = false;
555
      for (; !it.endDofIter(); it.nextDof()) {
556
	if (!isPrimal(feSpace, it.getDofIndex())) {
557
558
559
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
560
561
562
	    recvFromRank = true;
	    break;
	  }
563
	}
564
      }
565
566

      if (recvFromRank)
567
	stdMpi.recv(it.getRank());
568
    }
569

570
571
    stdMpi.startCommunication();

572
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
573
	 !it.end(); it.nextRank()) {
574
      int i = 0;
575
      for (; !it.endDofIter(); it.nextDof())
576
	if (!isPrimal(feSpace, it.getDofIndex()))
577
578
579
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
580
581
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
582
583
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
584
585
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
586

587
588
589
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

590
    int nRankLagrange = 0;
591
    DofMap& dualMap = dualDofMap[feSpace].getMap();
592
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
593
594
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
595
	int degree = boundaryDofRanks[feSpace][it->first].size();
596
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
597
      } else {
598
	lagrangeMap[feSpace].insertNonRankDof(it->first);
599
600
      }
    }
601
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
602
603
604
  }


605
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
606
  {
607
    FUNCNAME("PetscSolverFeti::createIndexB()");
608

609
    DOFAdmin* admin = feSpace->getAdmin();
610
611
612
613

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

615
    int nLocalInterior = 0;    
616
    for (int i = 0; i < admin->getUsedSize(); i++) {
617
      if (admin->isDofFree(i) == false && 
618
	  isPrimal(feSpace, i) == false &&
619
	  isDual(feSpace, i) == false) {
620
621
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
622

623
624
625
626
627
628
629
630
631
	  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);
632
633
634

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
635
	}
636
      }
637
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
638
    
639
640
    // === And finally, add the global indicies of all dual nodes. ===

641
642
643
644
645
646
647
648
649
650
    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);
      }
651
652
653
  }


654
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
655
656
657
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
658
659
    double wtime = MPI::Wtime();

660
661
    // === Create distributed matrix for Lagrange constraints. ===

662
    MatCreateMPIAIJ(mpiCommGlobal,
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
663
664
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
665
666
667
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

668
669
670
671
672
673
674
    // === 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
675
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
676
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
677

678
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
679
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
682
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
683
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
684

Thomas Witkowski's avatar
Thomas Witkowski committed
685
	// Copy set of all ranks that contain this dual node.
686
687
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
690
691

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
694
695
	
	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
696
697
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
698

699
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
700
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
701
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
702
	    index++;	      
703
704
705
706
707
708
709
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
712
713

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


717
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
718
  {
719
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
720

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

Thomas Witkowski's avatar
Thomas Witkowski committed
724
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
725
      
726
      VecCreateMPI(mpiCommGlobal, 
727
		   localDofMap.getRankDofs(), 
728
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
729
		   &(schurPrimalData.tmp_vec_b));
730
      VecCreateMPI(mpiCommGlobal, 
731
732
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
733
734
		   &(schurPrimalData.tmp_vec_primal));

735
      MatCreateShell(mpiCommGlobal,
736
737
738
739
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
740
741
742
743
744
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
745
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
746
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
747
748
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
749
750
      KSPSetFromOptions(ksp_schur_primal);
    } else {
751
752
      MSG("Create direct schur primal solver!\n");

753
754
      double wtime = MPI::Wtime();

755
756
757
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

758
759
760
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
761

Thomas Witkowski's avatar
Thomas Witkowski committed
762
      Mat matBPi;
763
      MatCreateMPIAIJ(mpiCommGlobal,
Thomas Witkowski's avatar
Thomas Witkowski committed
764
		      nRowsRankB, nRowsRankPrimal, 
765
		      nGlobalOverallInterior, nRowsOverallPrimal,
Thomas Witkowski's avatar
Thomas Witkowski committed
766
767
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
768

769
770
771
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
772

773
774
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

775
      for (int i = 0; i < nRowsRankB; i++) {
776
	PetscInt row = localDofMap.getStartDofs() + i;
Thomas Witkowski's avatar
Thomas Witkowski committed
777
	MatGetRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
778
779
780
781
782

	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
783
	MatRestoreRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
784
785
      }

786
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
787
		primalDofMap.getLocalDofs())
788
	("Should not happen!\n");
789

790
791
792
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
793
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
794
795
796
797
798
799
800
801

 	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
802
	subdomain->solve(tmpVec, tmpVec);
803

Thomas Witkowski's avatar
Thomas Witkowski committed
804
	PetscScalar *tmpValues;
805
	VecGetArray(tmpVec, &tmpValues);
806
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
807
	  MatSetValue(matBPi, 
808
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
809
810
811
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
812
813
814
815
816
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
817
818
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
819

Thomas Witkowski's avatar
Thomas Witkowski committed
820
821
      MatDuplicate(subdomain->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
      MatMatMult(subdomain->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
822
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
823
824

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
825
      MatDestroy(&matBPi);
826

827
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
828
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
829
		      SAME_NONZERO_PATTERN);
830
831
832
833
834
835
      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);
836
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
837

Thomas Witkowski's avatar
Thomas Witkowski committed
838
839
840
841
842
843
844
845
846
847
848
849
850
851
      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
852
    }
853
854
855
856
857
858
859
  }


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

860
    if (schurPrimalSolver == 0) {
861
      schurPrimalData.subSolver = NULL;
862

863
864
865
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
866
867
868
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
869
870
871
  }


872
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
873
874
875
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

876
877
    // === Create FETI-DP solver object. ===

878
    fetiData.mat_lagrange = &mat_lagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
879
    fetiData.subSolver = subdomain