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

namespace AMDiS {

  using namespace std;

22
#ifdef HAVE_PETSC_DEV 
Thomas Witkowski's avatar
Thomas Witkowski committed
23
  FetiStatisticsData PetscSolverFeti::fetiStatistics;
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

Thomas Witkowski's avatar
Thomas Witkowski committed
30
    double first = MPI::Wtime();
31
32
    void *ctx;
    MatShellGetContext(mat, &ctx);
33
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
34
35

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
36
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
37
38
39
40
    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);

Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
43
    PetscSolverFeti::fetiStatistics.nSchurPrimalApply++;
    PetscSolverFeti::fetiStatistics.timeSchurPrimalApply += MPI::Wtime() - first;

44
45
46
47
48
49
50
    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
51
52
    //    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)
53
54
55

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

58
59
    // tmp_vec_b0 = inv(K_BB) trans(L) x
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
Thomas Witkowski's avatar
Thomas Witkowski committed
60
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b0, data->tmp_vec_b0);
61

62
63
    // 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);
Thomas Witkowski's avatar
Thomas Witkowski committed
64
65

    double first = MPI::Wtime();
66
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
67
68
69
70
    PetscSolverFeti::fetiStatistics.nSchurPrimalSolve++;
    PetscSolverFeti::fetiStatistics.timeSchurPrimalSolve +=
      MPI::Wtime() - first;

71
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b1);
72

73
74
75
    // 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);
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
    PetscSolverFeti::fetiStatistics.nFetiApply++;

79
80
81
82
    return 0;
  }


83
  // y = PC * x
84
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
85
  {
86
    // Get data for the preconditioner
87
88
    void *ctx;
    PCShellGetContext(pc, &ctx);
89
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
90

91
    // Multiply with scaled Lagrange constraint matrix.
92
93
94
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


95
    // === Restriction of the B nodes to the boundary nodes. ===
96

97
98
99
100
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
101

102
103
104
105
106
107
    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];
108
109

    VecRestoreArray(data->tmp_vec_b, &local_b);
110
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
111
112


113
    // === K_DD - K_DI inv(K_II) K_ID ===
114

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

117
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
118
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
    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);
153
154


155
    // === Restriction of the B nodes to the boundary nodes. ===
156

157
158
159
160
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
161

162
    PetscScalar *local_b, *local_duals;
163
    VecGetArray(data->tmp_vec_b, &local_b);
164
    VecGetArray(data->tmp_vec_duals0, &local_duals);
165

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

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


    // === K_DD ===

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

177

178
    // === Prolongation from local dual nodes to B nodes.
179

180
181
182
183
184
    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];
185

186
187
    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
188

189
190

    // Multiply with scaled Lagrange constraint matrix.
191
192
193
194
195
196
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


197
198
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
      nComponents(-1),
      schurPrimalSolver(0)
201
202
203
204
205
206
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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


227
  void PetscSolverFeti::updateDofData()
228
229
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
230
231
232

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
233
   
234
235
236
237
238
239
240
    createPrimals();

    createDuals();

    createLagrange();

    createIndexB();
241
242
243
  }


244
  void PetscSolverFeti::createPrimals()
245
  {
246
    FUNCNAME("PetscSolverFeti::createPrimals()");  
247

248
249
250
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

251
252
253
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
254
255
256
257
258
259
    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);
260

261
262
263
264

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

265
    globalPrimalIndex.clear();
266
267
268
269
    nRankPrimals = 0;
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it)) {
	globalPrimalIndex[*it] = nRankPrimals;
270
271
272
	nRankPrimals++;
      }

273

274
275
276
    // === Get overall number of primals and rank's displacement in the ===
    // === numbering of the primals.                                    ===

277
    nOverallPrimals = 0;
278
    rStartPrimals = 0;
279
280
281
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankPrimals, rStartPrimals, nOverallPrimals);

282
283
284

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

285
286
287
288
    for (DofMapping::iterator it = globalPrimalIndex.begin();
	 it != globalPrimalIndex.end(); ++it)
      it->second += rStartPrimals;

289
290
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	nRankPrimals, nOverallPrimals);
291

292
293
294
295
296

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

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
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
    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());
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
    
345
346
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
347
348
349
350
351
352
353
354
355

    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
356
	if (globalPrimalIndex.count(**dofIt) == 0) {
357
358
359
360
361
362
	  boundaryDofRanks[**dofIt].insert(mpiRank);
	  boundaryDofRanks[**dofIt].insert(it->first);
	}
      }
    }

363
364
365
366

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

367
368
369
370
371
    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)
372
	if (globalPrimalIndex.count(**dofIt) == 0)
373
374
375
376
377
378
379
380
381
382
	  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)
383
	if (globalPrimalIndex.count(**dofIt) == 0) {
384
385
386
387
388
389
390
391
392
393
394
395
396
397
	  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)	
398
	if (globalPrimalIndex.count(**dofIt) == 0)
399
400
401
402
403
404
405
406
407
408
	  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();
409
    nRankB = nRankAllDofs - globalPrimalIndex.size();
410
411
412
413
414
415
416
417
418
419
420
    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) {
421
      if (globalPrimalIndex.count(**it) == 0) {
422
423
424
425
426
427
428
429
430
	duals.insert(**it);
	globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals;
	nRankDuals++;
      }
    }

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

431
432
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
	nRankDuals, nOverallDuals);
433
434
435
436
437
438
439
  }

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

440
441
442
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

443
444
445
446
447
448
449
450
451
    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;
      }
    }

452
453
454
455
456

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

457
    nOverallLagrange = 0;
458
    rStartLagrange = 0;
459
460
461
462
463
464
465
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankLagrange, rStartLagrange, nOverallLagrange);

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

466
467
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	nRankLagrange, nOverallLagrange);
468
469


470
    // === Communicate dofFirstLagrange to all other ranks. ===
471
472
473
474
475
476
477

    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) {
478
	if (globalPrimalIndex.count(**dofIt) == 0) {
479
480
481
482
483
484
485
486
487
488
489
490
	  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)
491
	if (globalPrimalIndex.count(**dofIt) == 0) {
492
493
494
495
496
497
498
499
500
501
502
503
504
505
	  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++)
506
	if (globalPrimalIndex.count(*(it->second[i])) == 0)
507
	  dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++];
508
    }     
509
510
511
512
513
514
515
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
516
    localIndexB.clear();
517
    DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin();
518
519
520
521

    // === 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.                                 ===
522
523
    
    for (int i = 0; i < admin->getUsedSize(); i++)
524
525
      if (admin->isDofFree(i) == false && globalPrimalIndex.count(i) == 0)
	if (duals.count(i) == 0 && globalPrimalIndex.count(i) == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
526
	  localIndexB[i] = -1;
527

528
529
530

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

531
    nLocalInterior = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
532
533
534
    for (DofMapping::iterator it = localIndexB.begin(); 
	 it != localIndexB.end(); ++it) {
      it->second = nLocalInterior;
535
      nLocalInterior++;
536
    }
537
    nLocalDuals = duals.size();
538

539
    TEST_EXIT_DBG(nLocalInterior + globalPrimalIndex.size() + duals.size() == 
540
541
542
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

543
544
545

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

546
547
    for (DofIndexSet::iterator it = duals.begin();
	 it != duals.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
548
      localIndexB[*it] = globalDualIndex[*it] - rStartB;
549
550
551
  }


552
  void PetscSolverFeti::createMatLagrange()
553
554
555
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

556
557
    // === Create distributed matrix for Lagrange constraints. ===

558
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
559
560
561
562
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
563
564
565
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

566
567
568
569
570
571
572
    // === 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.                                                     ===

573
574
575
576
    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");

577
      // Global index of the first Lagrange constriant for this node.
578
      int index = dofFirstLagrange[*it];
579
      // Copy set of all ranks that contain this dual node.
580
      vector<int> W(boundaryDofRanks[*it].begin(), boundaryDofRanks[*it].end());
581
      // Number of ranks that contain this dual node.
582
583
584
585
586
587
588
      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++) {
589
590
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
	    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);
  }


610
  void PetscSolverFeti::createSchurPrimalKsp()
611
612
613
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
      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, 
		   nRankPrimals * nComponents, nOverallPrimals * nComponents,
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
		     nRankPrimals * nComponents, nRankPrimals * nComponents,
		     nOverallPrimals * nComponents, nOverallPrimals * nComponents,
		     &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 {
642
643
      MSG("Create direct schur primal solver!\n");

644
645
646
647
648
649
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

      int nRowsRankPrimal = nRankPrimals * nComponents;
      int nRowsOverallPrimal = nOverallPrimals * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
650
651
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
652

Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
655
656
657
658
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
659

660
661
662
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
663

664
665
666
667
668
669
670
671
672
673
674
      int nLocalPrimals = globalPrimalIndex.size() * nComponents;
      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
675
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
		globalPrimalIndex.size() * nComponents)("Should not happen!\n");
      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
697
	PetscScalar *tmpValues;
698
699
	VecGetArray(tmpVec, &tmpValues);
	for (int i  = 0; i < nRankB * nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
702
703
704
	  MatSetValue(matBPi, 
		      rStartB * nComponents + i,
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
705
706
707
708
709
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
712
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
713
714
715
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
716
      MatDestroy(&matBPi);
717
718
719
720
721
722
723
724
725
726

      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
727

728
729
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
730
    }
731
732
733
734
735
736
737
  }


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

738
739
740
741
742
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
743

744
745
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
746

747
748
749
750
751
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
752
753
754
755
756
757
758
  }


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

759
760
    // === Create FETI-DP solver object. ===

761
762
763
    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
764
    fetiData.fetiSolver = this;
765
    fetiData.ksp_schur_primal = &ksp_schur_primal;
766

767
768
769
770
771
772
773
774
775
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
		 &(fetiData.tmp_vec_b0));
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
		 &(fetiData.tmp_vec_b1));
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankPrimals * nComponents, nOverallPrimals * nComponents,
		 &(fetiData.tmp_vec_primal));
776
777

    MatCreateShell(PETSC_COMM_WORLD,
778
779
		   nRankLagrange * nComponents, nRankLagrange * nComponents,
		   nOverallLagrange * nComponents, nOverallLagrange * nComponents,
780
		   &fetiData, &mat_feti);
781
782
783
784
785
786
787
    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);
788
789


790
    // === Create FETI-DP preconditioner object. ===
791

792
793
794
795
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
796

797
798
799
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
800
801
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
802
803
804
805
806
807
808
809
810
811
      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;
      
812
813
814
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiDirichletPreconData.tmp_vec_b));      
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
      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;

830
831
832
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiLumpedPreconData.tmp_vec_b));
833
834
835
836
837
838
839
840
      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);
      
841
842
      break;
    default:
843
844
      break;
    }
845
846
847
848
849
850
851
  }
  

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

852
853
    // === Destroy FETI-DP solver object. ===

854
855
856
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
857
    fetiData.fetiSolver = NULL;
858
    fetiData.ksp_schur_primal = PETSC_NULL;
859

860
861
    VecDestroy(&fetiData.tmp_vec_b0);
    VecDestroy(&fetiData.tmp_vec_b1);
862
    VecDestroy(&fetiData.tmp_vec_primal);
863
864
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
865
866


867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
    // === 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;
895
896
    default:
      break;
897
    }
898
899
900
901
902
903
904
905
906
  }


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

907
908
909
910
911
912
    PetscScalar *ttt;
    VecGetArray(vec_sol_primal, &ttt);
    for (int i = 0; i < nRankPrimals; i++)
      MSG("PRIM VAL = %f\n", ttt[i * 3 + 1]);
    VecRestoreArray(vec_sol_primal, &ttt);

913
    // === Get local part of the solution for B variables. ===
914
915
916
917
918

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


919
920
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
921
922
923
924
925
926
927
928
929
930
    
    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++) {
931
932
	  MSG("COPY GLOBAL PRIM %d of COMP %d TO LOCAL TMP %d\n",
	      it->second, i, counter);
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
	  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);

962
963
964
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
965
966
967
968

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

969
970
971
    for (int i = 0; i < globalPrimalIndex.size(); i++)
      MSG("TEST VAL %f\n", localSolPrimal[i * 3 + 1]);

972

973
    // === And copy from PETSc local vectors to the DOF vectors. ===
974
975
976
977

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

Thomas Witkowski's avatar
Thomas Witkowski committed
978
979
980
      for (DofMapping::iterator it = localIndexB.begin();
	   it != localIndexB.end(); ++it) {
	int petscIndex = it->second * nComponents + i;
981
982
983
984
985
986
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
      for (DofMapping::iterator it = globalPrimalIndex.begin();
	   it != globalPrimalIndex.end(); ++it) {
987
988
989
990
	if (i == 1)
	  MSG("AND SET THE VAL OF DOF %d TO VAL %f\n", 
	      it->first,
	      localSolPrimal[counter * nComponents + i]);
991
992
993
994
995
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

996
    MSG("VAL BOUNDS: %f %f\n", vec.getDOFVector(1)->min(), vec.getDOFVector(1)->max());
997
998
999

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
1000
    VecDestroy(&local_sol_primal);
1001
1002
1003
  }


1004
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
1005
1006
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
1007

1008
    nComponents = mat->getSize();
1009
1010
1011

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

1012
1013
    updateDofData();

1014
1015
1016
1017
1018
1019
1020

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

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
1021
    int nRowsInterior = nLocalInterior * nComponents;
1022
    int nRowsDual = nLocalDuals * nComponents;
1023

Thomas Witkowski's avatar
Thomas Witkowski committed
1024
1025
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
1026
1027

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1028
1029
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
1030
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
1031
1032

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1033
1034
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
1035
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
1036
1037

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1038
1039
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
1040
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
1041

1042

1043
1044
1045
1046
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1047
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
1048
1049
1050
1051
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1052
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
1053
1054
1055
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1056
		      nRowsInterior, nRowsDual, 30, PETSC_NULL,
1057
1058
1059
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1060
		      nRowsDual, nRowsInterior, 30, PETSC_NULL,
1061
1062
		      &mat_duals_interior);
    }
1063

1064
1065
    
    // === Prepare traverse of sequentially created matrices. ===
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080

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

1081
1082
1083
1084
1085
1086
1087
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102

    // === 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) {
1103

1104
	  bool rowPrimal = globalPrimalIndex.count(*cursor) != 0;
1105
  
1106
1107
	  cols.clear();
	  colsOther.clear();
1108
	  values.clear();	  
1109
	  valuesOther.clear();
1110
1111
1112
1113
1114
1115

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

1116
1117
1118
1119
1120
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

1121
	    bool colPrimal = globalPrimalIndex.count(col(*icursor)) != 0;
1122
1123

	    if (colPrimal) {
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
	      // 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));
1134
	      } else {
1135
1136
1137
1138
1139
1140
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      // Column is not a primal variable.

Thomas Witkowski's avatar
Thomas Witkowski committed
1141
	      TEST_EXIT_DBG(localIndexB.count(col(*icursor)))
1142
		("No global B index for DOF %d!\n", col(*icursor));
1143
	      
Thomas Witkowski's avatar
Thomas Witkowski committed
1144
1145
	      int colIndex = 
		(localIndexB[col(*icursor)] + rStartB) * nComponents + j;
1146
1147
1148
1149
1150
1151
1152

	      if (rowPrimal) {
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(colIndex);
		values.push_back(value(*icursor));
1153
1154
	      }
	    }
1155
1156
1157
1158
1159
1160



	    // === For preconditioner ===

	    if (!rowPrimal && !colPrimal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1161
1162
	      int rowIndex = localIndexB[*cursor];
	      int colIndex = localIndexB[col(*icursor)];
1163
1164
1165
		
	      if (rowIndex < nLocalInterior) {
		if (colIndex < nLocalInterior) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1166
		  int colIndex2 = localIndexB[col(*icursor)] * nComponents + j;
1167
1168
1169
1170

		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		} else {