PetscSolverFeti.cc 47.2 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/SubDomainSolver.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),
192
      subDomainSolver(NULL),
193
194
195
      meshLevel(0),
      rStartInterior(0),
      nGlobalOverallInterior(0)
196
197
198
199
200
201
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

    Parameters::get("parallel->multi level test", multiLevelTest);
221
222
    if (multiLevelTest)
      meshLevel = 1;
223
224
225
  }


226
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
227
  {
228
229
    FUNCNAME("PetscSolverFeti::initialize()");

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

233
234
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

235
236
237
238
239
240
241
242
243
    if (subDomainSolver == NULL) {
      if (meshLevel == 0) {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);
      } else {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, 
			      levelData.getMpiComm(meshLevel - 1),
			      levelData.getMpiComm(meshLevel));
244
	subDomainSolver->setLevel(meshLevel);
245
246
      }
    }
247

248
249
    primalDofMap.init(levelData,
		      feSpaces, meshDistributor->getFeSpaces(), 
250
		      true, true);
251
252
    primalDofMap.setComputeMatIndex(true);
    
253
254
    dualDofMap.init(levelData, 
		    feSpaces, meshDistributor->getFeSpaces(),
255
		    false, false);
256
257
    dualDofMap.setComputeMatIndex(true);

258
259
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
260
		     true, true);
261
    lagrangeMap.setComputeMatIndex(true);
262
263
264
265
266
267
268
269
270
271
272
273

    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);
    }
274
275

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

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


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

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

295
296
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

297
298
299
300
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
301
302
303
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.clear();

304
305
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
306

307
308
309
310
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);    
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
311
312
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
313
314
315
316

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

317
318
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
319

320
321
      createPrimals(feSpace);  

322
      createDuals(feSpace);      
323

324
325
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
328
329

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
330
    if (fetiPreconditioner != FETI_NONE)
331
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
332

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
333

334
335
336
337
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
338
    
339
340
341
    lagrangeMap.update();


342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      MSG("RANK %d FROM %d\n",
	  levelData.getMpiComm(1).Get_rank(),
	  levelData.getMpiComm(1).Get_size());

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

      mpi::getDofNumbering(mpiComm, groupRowsInterior,
			   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);

      MSG("COMM TEST FETI-DP: %d %d %d %d %d\n",
	  levelData.getMpiComm(1).Get_size(),
	  localDofMap.getRankDofs(),
	  localDofMap.getOverallDofs(),
	  nGlobalOverallInterior, rStartInterior);
    }

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

379
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
380
381
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
382
      
383
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
384
385
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
386

387
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
388
389
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
390

391
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
392
393
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
394
    }
395
396
397


    // If multi level test, inform sub domain solver about coarse space.
398
    subDomainSolver->setDofMapping(&primalDofMap, &localDofMap);
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
459
460
461
462
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
	} 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
479
480
481
482
483
484
485
486
487
      StdMpi<vector<int> > stdMpi(mpiComm);

      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
513


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

	  if (meshLevel == 0 ||
	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
531
532
533
534
	}
      }
    }

535
536
537
538

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

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

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

    stdMpi.updateSendDataSize();

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

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

569
570
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
583

584
585
586
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


602
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
603
  {
604
    FUNCNAME("PetscSolverFeti::createIndexB()");
605

606
    DOFAdmin* admin = feSpace->getAdmin();
607
608
609
610

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

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

620
621
622
623
624
625
626
627
628
	  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);
629
630
631

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

638
639
640
641
642
643
644
645
646
647
    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);
      }
648
649
650
  }


651
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
652
653
654
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

655
656
    // === Create distributed matrix for Lagrange constraints. ===

657
    MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
658
659
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
660
661
662
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

663
664
665
666
667
668
669
    // === 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
670
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
671
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
672

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

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

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

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

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
  }


708
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
709
  {
710
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
711

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

715
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
716
      
717
      VecCreateMPI(mpiComm, 
718
		   localDofMap.getRankDofs(), 
719
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
720
		   &(schurPrimalData.tmp_vec_b));
721
      VecCreateMPI(mpiComm, 
722
723
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
724
725
		   &(schurPrimalData.tmp_vec_primal));

726
      MatCreateShell(mpiComm,
727
728
729
730
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
733
734
735
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
736
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
737
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
738
739
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
740
741
      KSPSetFromOptions(ksp_schur_primal);
    } else {
742
743
      MSG("Create direct schur primal solver!\n");

744
745
      double wtime = MPI::Wtime();

746
747
748
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

749
750
751
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
752

Thomas Witkowski's avatar
Thomas Witkowski committed
753
      Mat matBPi;
754
      MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
755
		      nRowsRankB, nRowsRankPrimal, 
756
		      nGlobalOverallInterior, nRowsOverallPrimal,
Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
759

760
761
762
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
763

764
765
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

766
      for (int i = 0; i < nRowsRankB; i++) {
767
	PetscInt row = localDofMap.getStartDofs() + i;
768
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
769
770
771
772
773

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

774
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
775
776
      }

777
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
778
		primalDofMap.getLocalDofs())
779
	("Should not happen!\n");
780

781
782
783
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
784
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
785
786
787
788
789
790
791
792

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

793
	subDomainSolver->solve(tmpVec, tmpVec);
794

Thomas Witkowski's avatar
Thomas Witkowski committed
795
	PetscScalar *tmpValues;
796
	VecGetArray(tmpVec, &tmpValues);
797
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
798
	  MatSetValue(matBPi, 
799
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
800
801
802
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
803
804
805
806
807
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
808
809
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
810
811

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
812
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
813
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
814
815

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
816
      MatDestroy(&matBPi);
817
818

      MatInfo minfo;
819
      MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
820
821
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

822
      KSPCreate(mpiComm, &ksp_schur_primal);
823
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
824
		      SAME_NONZERO_PATTERN);
825
826
827
828
829
830
      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);
831
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
832

833
834
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
835
    }
836
837
838
839
840
841
842
  }


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

843
    if (schurPrimalSolver == 0) {
844
      schurPrimalData.subSolver = NULL;
845

846
847
848
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
849
850
851
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
852
853
854
  }


855
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
856
857
858
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

859
860
    // === Create FETI-DP solver object. ===

861
    fetiData.mat_lagrange = &mat_lagrange;
862
    fetiData.subSolver = subDomainSolver;
863
    fetiData.ksp_schur_primal = &ksp_schur_primal;
864

865
    VecCreateMPI(mpiComm, 
866
		 localDofMap.getRankDofs(), 
867
		 nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
868
		 &(fetiData.tmp_vec_b));
869
    VecCreateMPI(mpiComm,
870
871
		 lagrangeMap.getRankDofs(), 
		 lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
872
		 &(fetiData.tmp_vec_lagrange));
873
    VecCreateMPI(mpiComm, 
874
875
		 primalDofMap.getRankDofs(), 
		 primalDofMap.getOverallDofs(),
876
		 &(fetiData.tmp_vec_primal));
877

878
    MatCreateShell(mpiComm,
879
880
881
882
		   lagrangeMap.getRankDofs(), 
		   lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), 
		   lagrangeMap.getOverallDofs(),
883
		   &fetiData, &mat_feti);
884
885
886
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);