PetscSolverFeti.cc 53.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
//
// 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.


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

namespace AMDiS {

  using namespace std;

23
24
25
26
27
28
29
  // 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);
30
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
31
32

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
33
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
34
35
36
37
38
39
40
41
42
43
44
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(*(data->mat_primal_primal), x, y);
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
54
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
55

Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
57
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
60
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
61

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

    return 0;
  }


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

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


80
    // === Restriction of the B nodes to the boundary nodes. ===
81

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

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

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];
93
94

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


98
    // === K_DD - K_DI inv(K_II) K_ID ===
99

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

102
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
103
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
    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);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];

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


140
    // === Restriction of the B nodes to the boundary nodes. ===
141

142
143
144
145
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
146

147
    PetscScalar *local_b, *local_duals;
148
    VecGetArray(data->tmp_vec_b, &local_b);
149
    VecGetArray(data->tmp_vec_duals0, &local_duals);
150

151
152
    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];
153
154

    VecRestoreArray(data->tmp_vec_b, &local_b);
155
156
157
158
159
160
161
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

162

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

165
166
167
168
169
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];
170

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

174
175

    // Multiply with scaled Lagrange constraint matrix.
176
177
178
179
180
181
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


182
183
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
Thomas Witkowski's avatar
Thomas Witkowski committed
184
185
      nComponents(-1),
      schurPrimalSolver(0)
186
187
188
189
190
191
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
192
      MSG("Create FETI-DP solver with no preconditioner!\n");
193
194
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
195
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
196
197
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
198
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
199
200
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
203
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
204
205
206
207
208

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


212
  void PetscSolverFeti::updateDofData()
213
214
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
215
216
217

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
218
   
219
220
221
222
223
    createPrimals();

    createDuals();

    createLagrange();
224
    
225
    createIndexB();
226
227
228
  }


229
  void PetscSolverFeti::createPrimals()
230
  {
231
    FUNCNAME("PetscSolverFeti::createPrimals()");  
232

233
234
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

235
236
237
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

238
239
240
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
241
    DofContainerSet& vertices = 
242
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
243
244
245
246
    TEST_EXIT_DBG(vertices.size())("No primal vertices on this rank!\n");
    for (DofContainerSet::iterator it = vertices.begin(); 
	 it != vertices.end(); ++it)
      primals.insert(**it);
247

248
249
250
251

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

252
    primalDofMap.addFeSpace(feSpace);
253
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
254
255
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
256

257
    primalDofMap[feSpace].update();
258

259
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
260
261
	primalDofMap[feSpace].nRankDofs, 
	primalDofMap[feSpace].nOverallDofs);
262

263
264
265
266

    // === Communicate primal's global index from ranks that own the     ===
    // === primals to ranks that contain this primals but are not owning ===
    // === them.                                                         ===
267
268
    
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
269

270
    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
271
272
273
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
274
275
276
277
	if (primalDofMap[feSpace].isSet(it.getDofIndex())) {
	  DegreeOfFreedom d = primalDofMap[feSpace][it.getDofIndex()];
	  stdMpi.getSendData(it.getRank()).push_back(d);
	}
278

279
280
    stdMpi.updateSendDataSize();

281
282
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
283
      bool recvFromRank = false;
284
285
286
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
287
288
289
	  recvFromRank = true;
	  break;
	}
290
      }
291
292

      if (recvFromRank) 
293
	stdMpi.recv(it.getRank());
294
    }
295

296
297
    stdMpi.startCommunication();

298
299
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
300
      int i = 0;
301
302
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
303
304
305
306
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[i++];
	  primalDofMap[feSpace].insert(it.getDofIndex(), d);
	}
307
308
309
      }
    }

310
    TEST_EXIT_DBG(primals.size() == primalDofMap[feSpace].size())
311
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
312
       primals.size(), primalDofMap[feSpace].size());
313
314
315
316
317
318
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
319
320

    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
321
    
322
323
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
324
325
326

    boundaryDofRanks.clear();

327
328
329
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
330
	// If DOF is not primal, i.e., its a dual node
331
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
332
333
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
334
335
	}
      }
336
	
337
338
339
340

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

341
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
342
343
344
345

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
346
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
347
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
348
349
350

    stdMpi.updateSendDataSize();

351
352
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
353
      bool recvFromRank = false;
354
      for (; !it.endDofIter(); it.nextDof()) {
355
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
356
357
358
	  recvFromRank = true;
	  break;
	}
359
      }
360
361

      if (recvFromRank)
362
	stdMpi.recv(it.getRank());
363
    }
364

365
366
    stdMpi.startCommunication();

367
368
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
369
      int i = 0;
370
      for (; !it.endDofIter(); it.nextDof())
371
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
372
373
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
374
375
376
377
    }

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

378
    dualDofMap.addFeSpace(feSpace);
379

380
    int nRankAllDofs = feSpace->getAdmin()->getUsedDofs();
381
    nRankB = nRankAllDofs - primalDofMap[feSpace].size();
382
383
384
385
386
    nOverallB = 0;
    rStartB = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankB, rStartB, nOverallB);
    DofContainer allBoundaryDofs;
387
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
388
389
390
    int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size();

    for (DofContainer::iterator it = allBoundaryDofs.begin();
391
392
393
	 it != allBoundaryDofs.end(); ++it)
      if (primalDofMap[feSpace].isSet(**it) == false)
	dualDofMap[feSpace].insertRankDof(**it);
394

395
396
    dualDofMap[feSpace].update(false);
    dualDofMap[feSpace].addOffset(rStartB + nRankInteriorDofs);
397

398
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
399
400
	dualDofMap[feSpace].nRankDofs, 
	dualDofMap[feSpace].nOverallDofs);
401
402
403
404
405
406
407
  }

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

408
409
410
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

411
412
413
414
415
416
417
418
419
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
    dofFirstLagrange.addFeSpace(feSpace);

    int nRankLagrange = 0;
    DofMapping& dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
	dofFirstLagrange[feSpace].insert(it->first, nRankLagrange);
	int degree = boundaryDofRanks[it->first].size();
420
421
422
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
423
424
425
    dofFirstLagrange[feSpace].nRankDofs = nRankLagrange;
    dofFirstLagrange[feSpace].update();
    
426
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
427
428
	dofFirstLagrange[feSpace].nRankDofs, 
	dofFirstLagrange[feSpace].nOverallDofs);
429
430


431
    // === Communicate dofFirstLagrange to all other ranks. ===
432
433

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
434

435
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
436
	 !it.end(); it.nextRank()) {
437
      for (; !it.endDofIter(); it.nextDof())
438
439
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
	  TEST_EXIT_DBG(dofFirstLagrange[feSpace].isSet(it.getDofIndex()))
440
	    ("Should not happen!\n");
441
442
	  DegreeOfFreedom d = dofFirstLagrange[feSpace][it.getDofIndex()];
	  stdMpi.getSendData(it.getRank()).push_back(d);
443
	}
444
    }
445

446
447
    stdMpi.updateSendDataSize();

448
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
449
	 !it.end(); it.nextRank()) {
450
      bool recvData = false;
451
      for (; !it.endDofIter(); it.nextDof())
452
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
453
454
455
456
457
	  recvData = true;
	  break;
	}
	  
      if (recvData)
458
	stdMpi.recv(it.getRank());
459
460
461
462
    }

    stdMpi.startCommunication();

463
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
464
	 !it.end(); it.nextRank()) {
465
      int counter = 0;
466
467
468
469
470
471
      for (; !it.endDofIter(); it.nextDof()) {
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
	  dofFirstLagrange[feSpace].insert(it.getDofIndex(), d); 
	}
      }
472
    }
473
474
475
476
477
478
479
  }


  void PetscSolverFeti::createIndexB()
  {
    FUNCNAME("PetscSolverFeti::createIndeB()");

480
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
481
    localIndexB.clear();
482
    DOFAdmin* admin = feSpace->getAdmin();
483
484
485
486

    // === 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.                                 ===
487
488
    
    for (int i = 0; i < admin->getUsedSize(); i++)
489
490
491
492
      if (admin->isDofFree(i) == false && 
	  primalDofMap[feSpace].isSet(i) == false && 
	  dualDofMap[feSpace].isSet(i) == false)
	localIndexB[i] = -1;
493

494
495
496

    // === Get correct index for all interior nodes. ===

497
    nLocalInterior = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
498
499
500
    for (DofMapping::iterator it = localIndexB.begin(); 
	 it != localIndexB.end(); ++it) {
      it->second = nLocalInterior;
501
      nLocalInterior++;
502
    }
503
    nLocalDuals = dualDofMap[feSpace].size();
504

505
506
507
    TEST_EXIT_DBG(nLocalInterior + 
		  primalDofMap[feSpace].size() + 
		  dualDofMap[feSpace].size() == 
508
509
510
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

511
512
513

    // === And finally, add the global indicies of all dual nodes. ===

514
515
516
    for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      localIndexB[it->first] = it->second - rStartB;
517
518
519
  }


520
  void PetscSolverFeti::createMatLagrange()
521
522
523
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

524
525
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

526
527
    // === Create distributed matrix for Lagrange constraints. ===

528
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
529
		    dofFirstLagrange[feSpace].nRankDofs * nComponents, 
530
		    nRankB * nComponents,
531
		    dofFirstLagrange[feSpace].nOverallDofs * nComponents, 
532
		    nOverallB * nComponents,
533
534
535
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

536
537
538
539
540
541
542
    // === 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.                                                     ===

543
544
545
546
547
548
    DofMapping &dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      TEST_EXIT_DBG(dofFirstLagrange[feSpace].isSet(it->first))
	("Should not happen!\n");
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");
549

550
      // Global index of the first Lagrange constriant for this node.
551
      int index = dofFirstLagrange[feSpace][it->first];
552
      // Copy set of all ranks that contain this dual node.
553
554
      vector<int> W(boundaryDofRanks[it->first].begin(), 
		    boundaryDofRanks[it->first].end());
555
      // Number of ranks that contain this dual node.
556
557
558
559
      int degree = W.size();

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
560
561
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
562
563
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
564
	      int colIndex = it->second * nComponents + k;
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
	    }
	  }

	  index++;
	}
      }
    }

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


581
  void PetscSolverFeti::createSchurPrimalKsp()
582
583
584
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

585
586
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
590
591
592
593
594
595
596
597
598
      schurPrimalData.mat_primal_primal = &mat_primal_primal;
      schurPrimalData.mat_primal_b = &mat_primal_b;
      schurPrimalData.mat_b_primal = &mat_b_primal;
      schurPrimalData.fetiSolver = this;
      
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
599
		   primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
600
601
602
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
603
604
		     primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nRankDofs * nComponents,
		     primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
605
606
607
608
609
610
611
612
613
614
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
    } else {
615
616
      MSG("Create direct schur primal solver!\n");

617
618
619
620
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

621
622
      int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents;
      int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
623
624
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
625

Thomas Witkowski's avatar
Thomas Witkowski committed
626
627
628
629
630
631
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
632

633
634
635
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
636

637
638
639
640
641
642
643
644
645
646
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

      for (int i = 0; i < nRankB * nComponents; i++) {
	PetscInt row = rStartB * nComponents + i;
	MatGetRow(mat_b_primal, row, &nCols, &cols, &values);

	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
647
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
648
649
650
651
652
653
      }

      int maxLocalPrimal = mat_b_primal_cols.size();
      mpi::globalMax(maxLocalPrimal);

      TEST_EXIT(mat_b_primal_cols.size() ==
654
		primalDofMap[feSpace].size() * nComponents)("Should not happen!\n");
655
656
657
658
659
660
661
662
663
664
665
666
667
668
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
	VecCreateSeq(PETSC_COMM_SELF, nRankB * nComponents, &tmpVec);

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

       	KSPSolve(ksp_b, tmpVec, tmpVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
669
	PetscScalar *tmpValues;
670
671
	VecGetArray(tmpVec, &tmpValues);
	for (int i  = 0; i < nRankB * nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
672
673
674
675
676
	  MatSetValue(matBPi, 
		      rStartB * nComponents + i,
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
677
678
679
680
681
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
684
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
685
686
687
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
688
      MatDestroy(&matBPi);
689
690
691
692
693
694
695
696
697
698

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

      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_primal_primal, 
		      mat_primal_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
699

700
701
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
702
    }
703
704
705
706
707
708
709
  }


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

710
711
712
713
714
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
715

716
717
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
718

719
720
721
722
723
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
724
725
726
727
728
729
730
  }


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

731
732
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

733
734
    // === Create FETI-DP solver object. ===

735
736
737
    fetiData.mat_primal_b = &mat_primal_b;
    fetiData.mat_b_primal = &mat_b_primal;
    fetiData.mat_lagrange = &mat_lagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
738
    fetiData.fetiSolver = this;
739
    fetiData.ksp_schur_primal = &ksp_schur_primal;
740

741
742
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
745
746
		 dofFirstLagrange[feSpace].nRankDofs * nComponents, 
		 dofFirstLagrange[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
747
		 &(fetiData.tmp_vec_lagrange));
748
    VecCreateMPI(PETSC_COMM_WORLD, 
749
		 primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
750
		 &(fetiData.tmp_vec_primal));
751
752

    MatCreateShell(PETSC_COMM_WORLD,
753
754
755
756
		   dofFirstLagrange[feSpace].nRankDofs * nComponents, 
		   dofFirstLagrange[feSpace].nRankDofs * nComponents,
		   dofFirstLagrange[feSpace].nOverallDofs * nComponents, 
		   dofFirstLagrange[feSpace].nOverallDofs * nComponents,
757
		   &fetiData, &mat_feti);
758
759
760
761
762
763
764
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


    KSPCreate(PETSC_COMM_WORLD, &ksp_feti);
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_feti, "solver_feti_");
    KSPSetFromOptions(ksp_feti);
765
766


767
    // === Create FETI-DP preconditioner object. ===
768

769
770
771
772
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
773

774
775
776
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
777
778
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
779
780
781
782
783
784
785
786
787
788
      KSPSetOptionsPrefix(ksp_interior, "solver_interior_");
      KSPSetFromOptions(ksp_interior);
            
      fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
      fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
      fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
      fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
      fetiDirichletPreconData.ksp_interior = &ksp_interior;
      
789
790
791
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiDirichletPreconData.tmp_vec_b));      
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_interior));
      
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
      
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;

807
808
809
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiLumpedPreconData.tmp_vec_b));
810
811
812
813
814
815
816
817
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals1));

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
818
819
      break;
    default:
820
821
      break;
    }
822
823
824
825
826
827
828
  }
  

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

829
830
    // === Destroy FETI-DP solver object. ===

831
832
833
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
834
    fetiData.fetiSolver = NULL;
835
    fetiData.ksp_schur_primal = PETSC_NULL;
836

Thomas Witkowski's avatar
Thomas Witkowski committed
837
838
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
839
    VecDestroy(&fetiData.tmp_vec_primal);
840
841
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
842
843


844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
    // === Destroy FETI-DP preconditioner object. ===

    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPDestroy(&ksp_interior);

      fetiDirichletPreconData.mat_lagrange_scaled = NULL;
      fetiDirichletPreconData.mat_interior_interior = NULL;
      fetiDirichletPreconData.mat_duals_duals = NULL;
      fetiDirichletPreconData.mat_interior_duals = NULL;
      fetiDirichletPreconData.mat_duals_interior = NULL;
      fetiDirichletPreconData.ksp_interior = NULL;
      
      VecDestroy(&fetiDirichletPreconData.tmp_vec_b);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals0);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals1);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_interior);
      MatDestroy(&mat_lagrange_scaled);
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = NULL;
      fetiLumpedPreconData.mat_duals_duals = NULL;

      VecDestroy(&fetiLumpedPreconData.tmp_vec_b);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1);
      break;
872
873
    default:
      break;
874
    }
875
876
877
878
879
880
881
882
883
  }


  void PetscSolverFeti::recoverSolution(Vec &vec_sol_b,
					Vec &vec_sol_primal,
					SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::recoverSolution()");

884
885
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

886
    // === Get local part of the solution for B variables. ===
887
888
889
890
891

    PetscScalar *localSolB;
    VecGetArray(vec_sol_b, &localSolB);


892
893
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
894
895
    
    vector<PetscInt> globalIsIndex, localIsIndex;
896
897
    globalIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
    localIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
898
899
900

    {
      int counter = 0;
901
902
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
	for (int i = 0; i < nComponents; i++) {
	  globalIsIndex.push_back(it->second * nComponents + i);
	  localIsIndex.push_back(counter++);
	}
      }
    }
    
    IS globalIs, localIs;
    ISCreateGeneral(PETSC_COMM_SELF, 
		    globalIsIndex.size(), 
		    &(globalIsIndex[0]),
		    PETSC_USE_POINTER,
		    &globalIs);

    ISCreateGeneral(PETSC_COMM_SELF, 
		    localIsIndex.size(), 
		    &(localIsIndex[0]),
		    PETSC_USE_POINTER,
		    &localIs);

    Vec local_sol_primal;
    VecCreateSeq(PETSC_COMM_SELF, localIsIndex.size(), &local_sol_primal);

    VecScatter primalScatter;
    VecScatterCreate(vec_sol_primal, globalIs, local_sol_primal, localIs, &primalScatter);
    VecScatterBegin(primalScatter, vec_sol_primal, local_sol_primal, 
		    INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(primalScatter, vec_sol_primal, local_sol_primal, 
		  INSERT_VALUES, SCATTER_FORWARD);

933
934
935
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
936
937
938
939

    PetscScalar *localSolPrimal;
    VecGetArray(local_sol_primal, &localSolPrimal);

940
    // === And copy from PETSc local vectors to the DOF vectors. ===
941
942
943
944

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));

Thomas Witkowski's avatar
Thomas Witkowski committed
945
946
947
      for (DofMapping::iterator it = localIndexB.begin();
	   it != localIndexB.end(); ++it) {
	int petscIndex = it->second * nComponents + i;
948
949
950
951
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
952
953
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
954
955
956
957
958
959
960
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
961
    VecDestroy(&local_sol_primal);
962
963
964
  }


965
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
966
967
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
968

969
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
970
    nComponents = mat->getSize();
971
972
973

    // === Create all sets and indices. ===

974
975
976
    primalDofMap.setMpiComm(mpiComm);
    dualDofMap.setMpiComm(mpiComm);
    dofFirstLagrange.setMpiComm(mpiComm);
977
978
    updateDofData();

979
980
981
982
    // === Create matrices for the FETI-DP method. ===

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
983
984
    int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents;
    int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents;
985
    int nRowsInterior = nLocalInterior * nComponents;
986
    int nRowsDual = nLocalDuals * nComponents;
987

Thomas Witkowski's avatar
Thomas Witkowski committed
988
989
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
990
991

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
992
993
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,