PetscSolverFeti.cc 45.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
//
// 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"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

  using namespace std;


#ifdef HAVE_PETSC_DEV 
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
33
34
35
36
37
38
39
40
41
42
43
44
45

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);

    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)
  {
46
47
    //    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)
48
49
50

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

53
54
55
    // tmp_vec_b0 = inv(K_BB) trans(L) x
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b0, data->tmp_vec_b0);
56

57
58
    // tmp_vec_b1 = inv(K_BB) K_BPi  inv(S_PiPi) K_PiB tmp_vec_b0
    MatMult(*(data->mat_primal_b), data->tmp_vec_b0, data->tmp_vec_primal);
59
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
60
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b1);
61

62
63
64
    // y = L (tmp_vec_b0 + tmp_vec_b1)
    VecAXPBY(data->tmp_vec_b0, 1.0, 1.0, data->tmp_vec_b1);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y);
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
93
94
    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];
95
96

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


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

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

104
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
105
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
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
138
139
    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);
140
141


142
    // === Restriction of the B nodes to the boundary nodes. ===
143

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

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

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

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


    // === K_DD ===

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

164

165
    // === Prolongation from local dual nodes to B nodes.
166

167
168
169
170
171
    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];
172

173
174
    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
175

176
177

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

    return 0;
  }


184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
      nComponents(-1)
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
      fetiPreconditioner = FETI_LUMPED;
    } else {
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", preconditionerName.c_str());
    }
  }


204
  void PetscSolverFeti::updateDofData()
205
206
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
207
208
209

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
210
   
211
212
213
214
215
216
217
    createPrimals();

    createDuals();

    createLagrange();

    createIndexB();
218
219
220
  }


221
  void PetscSolverFeti::createPrimals()
222
  {
223
    FUNCNAME("PetscSolverFeti::createPrimals()");  
224

225
226
227
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

228
229
230
231
232
233
234
    primals.clear();
    DofContainerSet& vertices = 
      meshDistributor->getBoundaryDofInfo().geoDofs[VERTEX];
    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);
235

236
237
238
239

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

240
    globalPrimalIndex.clear();
241
242
243
244
    nRankPrimals = 0;
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it)) {
	globalPrimalIndex[*it] = nRankPrimals;
245
246
247
	nRankPrimals++;
      }

248

249
250
251
    // === Get overall number of primals and rank's displacement in the ===
    // === numbering of the primals.                                    ===

252
    nOverallPrimals = 0;
253
    rStartPrimals = 0;
254
255
256
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankPrimals, rStartPrimals, nOverallPrimals);

257
258
259

    // === Create global primal index for all primals. ===

260
261
262
263
    for (DofMapping::iterator it = globalPrimalIndex.begin();
	 it != globalPrimalIndex.end(); ++it)
      it->second += rStartPrimals;

264
265
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	nRankPrimals, nOverallPrimals);
266

267
268
269
270
271

    // === Communicate primal's global index from ranks that own the     ===
    // === primals to ranks that contain this primals but are not owning ===
    // === them.                                                         ===

272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (globalPrimalIndex.count(**dofIt))
	  stdMpi.getSendData(it->first).push_back(globalPrimalIndex[**dofIt]);
    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvFromRank = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false) {
	  recvFromRank = true;
	  break;
	}

      if (recvFromRank) 
	stdMpi.recv(it->first);
    }
    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int i = 0;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false)
	  globalPrimalIndex[**dofIt] = stdMpi.getRecvData(it->first)[i++];
      }
    }

    TEST_EXIT_DBG(primals.size() == globalPrimalIndex.size())
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
       primals.size(), globalPrimalIndex.size());


    TEST_EXIT_DBG(nOverallPrimals > 0)
      ("There are zero primal nodes in domain!\n");
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
    
324
325
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341

    boundaryDofRanks.clear();

    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it) {
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	// If DOF is not primal, i.e., its a dual node
	if (primals.count(**dofIt) == 0) {
	  boundaryDofRanks[**dofIt].insert(mpiRank);
	  boundaryDofRanks[**dofIt].insert(it->first);
	}
      }
    }

342
343
344
345

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

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
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) == 0)
	  stdMpi.getSendData(it->first).push_back(boundaryDofRanks[**dofIt]);

    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvFromRank = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) == 0) {
	  recvFromRank = true;
	  break;
	}

      if (recvFromRank)
	stdMpi.recv(it->first);
    }
    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int i = 0;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)	
	if (primals.count(**dofIt) == 0)
	  boundaryDofRanks[**dofIt] = stdMpi.getRecvData(it->first)[i++];	      
    }


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

    duals.clear();
    globalDualIndex.clear();

    int nRankAllDofs = meshDistributor->getFeSpace()->getAdmin()->getUsedDofs();
    nRankB = nRankAllDofs - primals.size();
    nOverallB = 0;
    rStartB = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankB, rStartB, nOverallB);
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(allBoundaryDofs);
    int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size();

    int nRankDuals = 0;
    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it) {
      if (primals.count(**it) == 0) {
	duals.insert(**it);
	globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals;
	nRankDuals++;
      }
    }

    int nOverallDuals = nRankDuals;
    mpi::globalAdd(nOverallDuals);

410
411
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
	nRankDuals, nOverallDuals);
412
413
414
415
416
417
418
  }

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

419
420
421
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

422
423
424
425
426
427
428
429
430
    nRankLagrange = 0;
    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
      if (meshDistributor->getIsRankDof(*it)) {
	dofFirstLagrange[*it] = nRankLagrange;
	int degree = boundaryDofRanks[*it].size();
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }

431
432
433
434
435

    // === Get the overall number of Lagrange constraints and create the ===
    // === mapping dofFirstLagrange, that defines for each dual boundary ===
    // === node the first Lagrange constraint global index.              ===

436
    nOverallLagrange = 0;
437
    rStartLagrange = 0;
438
439
440
441
442
443
444
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankLagrange, rStartLagrange, nOverallLagrange);

    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it))
	dofFirstLagrange[*it] += rStartLagrange;

445
446
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	nRankLagrange, nOverallLagrange);
447
448


449
    // === Communicate dofFirstLagrange to all other ranks. ===
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (primals.count(**dofIt) == 0) {
	  TEST_EXIT_DBG(dofFirstLagrange.count(**dofIt))("Should not happen!\n");
	  stdMpi.getSendData(it->first).push_back(dofFirstLagrange[**dofIt]);
	}
      }
    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvData = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) == 0) {
	  recvData = true;
	  break;
	}
	  
      if (recvData)
	stdMpi.recv(it->first);
    }

    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int counter = 0;
      for (unsigned int i = 0; i < it->second.size(); i++)
	if (primals.count(*(it->second[i])) == 0)
	  dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++];
487
    }     
488
489
490
491
492
493
494
495
496
  }


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

    globalIndexB.clear();
    DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin();
497
498
499
500

    // === 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.                                 ===
501
502
503
504
505
506
    
    for (int i = 0; i < admin->getUsedSize(); i++)
      if (admin->isDofFree(i) == false && primals.count(i) == 0)
	if (duals.count(i) == 0 && primals.count(i) == 0)
	  globalIndexB[i] = -1;

507
508
509

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

510
    nLocalInterior = 0;
511
512
    for (DofMapping::iterator it = globalIndexB.begin(); 
	 it != globalIndexB.end(); ++it) {
513
514
      it->second = nLocalInterior + rStartB;
      nLocalInterior++;
515
    }
516
    nLocalDuals = duals.size();
517

518
    TEST_EXIT_DBG(nLocalInterior + primals.size() + duals.size() == 
519
520
521
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

522
523
524

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

525
526
527
    for (DofIndexSet::iterator it = duals.begin();
	 it != duals.end(); ++it)
      globalIndexB[*it] = globalDualIndex[*it];
528
529
530
  }


531
  void PetscSolverFeti::createMatLagrange()
532
533
534
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

535
536
    // === Create distributed matrix for Lagrange constraints. ===

537
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
538
539
540
541
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
542
543
544
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

545
546
547
548
549
550
551
    // === 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.                                                     ===

552
553
554
555
    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
      TEST_EXIT_DBG(dofFirstLagrange.count(*it))("Should not happen!\n");
      TEST_EXIT_DBG(boundaryDofRanks.count(*it))("Should not happen!\n");

556
      // Global index of the first Lagrange constriant for this node.
557
      int index = dofFirstLagrange[*it];
558
      // Copy set of all ranks that contain this dual node.
559
      vector<int> W(boundaryDofRanks[*it].begin(), boundaryDofRanks[*it].end());
560
      // Number of ranks that contain this dual node.
561
562
563
564
565
566
567
      int degree = W.size();

      TEST_EXIT_DBG(globalDualIndex.count(*it))("Should not happen!\n");
      int dualCol = globalDualIndex[*it];

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
568
569
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
	      int colIndex = dualCol * nComponents + k;
	      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);
  }


589
  void PetscSolverFeti::createSchurPrimalKsp()
590
591
592
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

593
594
595
596
    schurPrimalData.mat_primal_primal = &mat_primal_primal;
    schurPrimalData.mat_primal_b = &mat_primal_b;
    schurPrimalData.mat_b_primal = &mat_b_primal;
    schurPrimalData.ksp_b = &ksp_b;
597

598
599
    VecDuplicate(f_b, &(schurPrimalData.tmp_vec_b));
    VecDuplicate(f_primal, &(schurPrimalData.tmp_vec_primal));
600
601
602
603

    MatCreateShell(PETSC_COMM_WORLD,
		   nRankPrimals * nComponents, nRankPrimals * nComponents,
		   nOverallPrimals * nComponents, nOverallPrimals * nComponents,
604
		   &schurPrimalData, 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
		   &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);
  }


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

620
621
622
623
    schurPrimalData.mat_primal_primal = PETSC_NULL;
    schurPrimalData.mat_primal_b = PETSC_NULL;
    schurPrimalData.mat_b_primal = PETSC_NULL;
    schurPrimalData.ksp_b = PETSC_NULL;
624

625
626
    VecDestroy(&schurPrimalData.tmp_vec_b);
    VecDestroy(&schurPrimalData.tmp_vec_primal);
627

628
629
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
630
631
632
633
634
635
636
  }


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

637
638
    // === Create FETI-DP solver object. ===

639
640
641
642
643
644
    fetiData.mat_primal_primal = &mat_primal_primal;
    fetiData.mat_primal_b = &mat_primal_b;
    fetiData.mat_b_primal = &mat_b_primal;
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.ksp_b = &ksp_b;
    fetiData.ksp_schur_primal = &ksp_schur_primal;
645

646
647
    VecDuplicate(f_b, &(fetiData.tmp_vec_b0));
    VecDuplicate(f_b, &(fetiData.tmp_vec_b1));
648
    VecDuplicate(f_primal, &(fetiData.tmp_vec_primal));
649
650

    MatCreateShell(PETSC_COMM_WORLD,
651
652
		   nRankLagrange * nComponents, nRankLagrange * nComponents,
		   nOverallLagrange * nComponents, nOverallLagrange * nComponents,
653
		   &fetiData, &mat_feti);
654
655
656
657
658
659
660
    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);
661
662


663
    // === Create FETI-DP preconditioner object. ===
664

665
666
667
668
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
669

670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, SAME_NONZERO_PATTERN);
      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;
      
      VecDuplicate(f_b, &(fetiDirichletPreconData.tmp_vec_b));      
      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;

      VecDuplicate(f_b, &(fetiLumpedPreconData.tmp_vec_b));
      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);
      
      break;
    }
711
712
713
714
715
716
717
  }
  

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

718
719
    // === Destroy FETI-DP solver object. ===

720
721
722
723
724
725
    fetiData.mat_primal_primal = PETSC_NULL;
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
    fetiData.ksp_b = PETSC_NULL;
    fetiData.ksp_schur_primal = PETSC_NULL;
726

727
728
    VecDestroy(&fetiData.tmp_vec_b0);
    VecDestroy(&fetiData.tmp_vec_b1);
729
    VecDestroy(&fetiData.tmp_vec_primal);
730
731
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
732
733


734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
    // === 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;
    }
763
764
765
766
767
768
769
770
771
  }


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

772
    // === Get local part of the solution for B variables. ===
773
774
775
776
777

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


778
779
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
    
    vector<PetscInt> globalIsIndex, localIsIndex;
    globalIsIndex.reserve(globalPrimalIndex.size() * nComponents);
    localIsIndex.reserve(globalPrimalIndex.size() * nComponents);

    {
      int counter = 0;
      for (DofMapping::iterator it = globalPrimalIndex.begin();
	   it != globalPrimalIndex.end(); ++it) {
	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);

819
820
821
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
822
823
824
825
826

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


827
    // === And copy from PETSc local vectors to the DOF vectors. ===
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849

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

      for (DofMapping::iterator it = globalIndexB.begin();
	   it != globalIndexB.end(); ++it) {
	int petscIndex = (it->second - rStartB) * nComponents + i;
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
      for (DofMapping::iterator it = globalPrimalIndex.begin();
	   it != globalPrimalIndex.end(); ++it) {
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }



    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
850
    VecDestroy(&local_sol_primal);
851
852
853
  }


854
855
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat, 
					SystemVector *vec)
856
857
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
858

859
860
861
862
    nComponents = vec->getSize();

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

863
864
    updateDofData();

865
866
867
868
869
870
871

    // === Create matrices for the FETI-DP method. ===

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
872
    int nRowsInterior = nLocalInterior * nComponents;
873
    int nRowsDual = nLocalDuals * nComponents;
874
875

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
876
		    nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB,
877
		    30, PETSC_NULL, 0, PETSC_NULL, &mat_b_b);
878
879

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
880
881
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
882
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
883
884

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
885
886
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
887
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
888
889

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
890
891
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
892
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_b);
893

894

895
896
897
898
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
899
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
900
901
902
903
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
904
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
905
906
907
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
908
		      nRowsInterior, nRowsDual, 30, PETSC_NULL,
909
910
911
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
912
		      nRowsDual, nRowsInterior, 30, PETSC_NULL,
913
914
		      &mat_duals_interior);
    }
915

916
917
    
    // === Prepare traverse of sequentially created matrices. ===
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

    vector<int> cols, colsOther;
    vector<double> values, valuesOther;
    cols.reserve(300);
    colsOther.reserve(300);
    values.reserve(300);
    valuesOther.reserve(300);

933
934
935
936
937
938
939
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

940
941
942
943
944
945
946
947
948
949
950
951
952
953
954

    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	if (!(*mat)[i][j])
	  continue;

	traits::col<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
	traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
	
	// Traverse all rows.
	for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()), 
	       cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
955

956
	  bool rowPrimal = primals.count(*cursor) != 0;
957
  
958
959
	  cols.clear();
	  colsOther.clear();
960
	  values.clear();	  
961
	  valuesOther.clear();
962
963
964
965
966
967

	  colsLocal.clear();
	  colsLocalOther.clear();
	  valuesLocal.clear();
	  valuesLocalOther.clear();

968
969
970
971
972
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

973
974
975
	    bool colPrimal = primals.count(col(*icursor)) != 0;

	    if (colPrimal) {
976
977
978
979
980
981
982
983
984
985
	      // Column is a primal variable.

	      TEST_EXIT_DBG(globalPrimalIndex.count(col(*icursor)))
		("No global primal index for DOF %d!\n", col(*icursor));
	      
	      int colIndex = globalPrimalIndex[col(*icursor)] * nComponents + j;
	      
	      if (rowPrimal) {
		cols.push_back(colIndex);
		values.push_back(value(*icursor));
986
	      } else {
987
988
989
990
991
992
993
994
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      // Column is not a primal variable.

	      TEST_EXIT_DBG(globalIndexB.count(col(*icursor)))
		("No global B index for DOF %d!\n", col(*icursor));
995
	      
996
997
998
999
1000
	      int colIndex = globalIndexB[col(*icursor)] * nComponents + j;

	      if (rowPrimal) {
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));