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

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
32
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
33
34
35
36
37
38
39
40
41
42
43
    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)
  {
44
45
    //    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)
46
47
48

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

Thomas Witkowski's avatar
Thomas Witkowski committed
51
52
53
    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
54

Thomas Witkowski's avatar
Thomas Witkowski committed
55
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
56
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
57
58
59
    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);
60

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

    return 0;
  }


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

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


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

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

86
87
88
89
90
91
    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];
92
93

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


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

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

101
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
102
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
103
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
    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);
137
138


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

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

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

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

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


    // === K_DD ===

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

161

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

164
165
166
167
168
    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];
169

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

173
174

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

    return 0;
  }


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

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

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


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

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

    createDuals();

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


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

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

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

237
238
239
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
240
    DofContainerSet& vertices = 
241
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
242
243
244
245
    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);
246

247
248
249
250

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

251
    globalPrimalIndex.clear();
252
253
    nRankPrimals = 0;
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
254
      if (meshDistributor->getIsRankDof(feSpace, *it)) {
255
	globalPrimalIndex[*it] = nRankPrimals;
256
257
258
	nRankPrimals++;
      }

259

260
261
262
    // === Get overall number of primals and rank's displacement in the ===
    // === numbering of the primals.                                    ===

263
    nOverallPrimals = 0;
264
    rStartPrimals = 0;
265
266
267
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankPrimals, rStartPrimals, nOverallPrimals);

268
269
270

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

271
272
273
274
    for (DofMapping::iterator it = globalPrimalIndex.begin();
	 it != globalPrimalIndex.end(); ++it)
      it->second += rStartPrimals;

275
276
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	nRankPrimals, nOverallPrimals);
277

278
279
280
281

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

285
    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
286
287
288
289
290
291
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
	if (globalPrimalIndex.count(it.getDofIndex()))
	  stdMpi.getSendData(it.getRank()).push_back(globalPrimalIndex[it.getDofIndex()]);

292
293
    stdMpi.updateSendDataSize();

294
295
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
296
      bool recvFromRank = false;
297
298
299
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
300
301
302
	  recvFromRank = true;
	  break;
	}
303
      }
304
305

      if (recvFromRank) 
306
	stdMpi.recv(it.getRank());
307
    }
308

309
310
    stdMpi.startCommunication();

311
312
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
313
      int i = 0;
314
315
316
317
318
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false)
	  globalPrimalIndex[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
319
320
321
322
323
324
325
326
327
328
329
330
      }
    }

    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()");
331
332

    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
333
    
334
335
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
336
337
338

    boundaryDofRanks.clear();

339
340
341
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
342
	// If DOF is not primal, i.e., its a dual node
343
344
345
	if (globalPrimalIndex.count(it.getDofIndex()) == 0) {
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
346
347
	}
      }
348
	
349
350
351
352

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

353
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
354
355
356
357
358
359

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
	if (globalPrimalIndex.count(it.getDofIndex()) == 0)
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
360
361
362

    stdMpi.updateSendDataSize();

363
364
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
365
      bool recvFromRank = false;
366
367
      for (; !it.endDofIter(); it.nextDof()) {
	if (globalPrimalIndex.count(it.getDofIndex()) == 0) {
368
369
370
	  recvFromRank = true;
	  break;
	}
371
      }
372
373

      if (recvFromRank)
374
	stdMpi.recv(it.getRank());
375
    }
376

377
378
    stdMpi.startCommunication();

379
380
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
381
      int i = 0;
382
383
384
385
      for (; !it.endDofIter(); it.nextDof())
	if (globalPrimalIndex.count(it.getDofIndex()) == 0)
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
386
387
388
389
390
391
392
    }

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

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

393
    int nRankAllDofs = feSpace->getAdmin()->getUsedDofs();
394
    nRankB = nRankAllDofs - globalPrimalIndex.size();
395
396
397
398
399
    nOverallB = 0;
    rStartB = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankB, rStartB, nOverallB);
    DofContainer allBoundaryDofs;
400
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
401
402
403
404
405
    int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size();

    int nRankDuals = 0;
    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it) {
406
      if (globalPrimalIndex.count(**it) == 0) {
407
408
409
410
411
412
413
414
415
	duals.insert(**it);
	globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals;
	nRankDuals++;
      }
    }

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

416
417
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
	nRankDuals, nOverallDuals);
418
419
420
421
422
423
424
  }

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

425
426
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

427
428
429
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

430
431
    nRankLagrange = 0;
    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
432
      if (meshDistributor->getIsRankDof(feSpace, *it)) {
433
434
435
436
437
438
	dofFirstLagrange[*it] = nRankLagrange;
	int degree = boundaryDofRanks[*it].size();
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }

439
440
441
442
443

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

444
    nOverallLagrange = 0;
445
    rStartLagrange = 0;
446
447
448
449
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankLagrange, rStartLagrange, nOverallLagrange);

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

453
454
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	nRankLagrange, nOverallLagrange);
455
456


457
    // === Communicate dofFirstLagrange to all other ranks. ===
458
459

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

461
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
462
	 !it.end(); it.nextRank()) {
463
464
465
466
467
      for (; !it.endDofIter(); it.nextDof())
	if (globalPrimalIndex.count(it.getDofIndex()) == 0) {
	  TEST_EXIT_DBG(dofFirstLagrange.count(it.getDofIndex()))
	    ("Should not happen!\n");
	  stdMpi.getSendData(it.getRank()).push_back(dofFirstLagrange[it.getDofIndex()]);
468
	}
469
    }
470

471
472
    stdMpi.updateSendDataSize();

473
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
474
	 !it.end(); it.nextRank()) {
475
      bool recvData = false;
476
477
      for (; !it.endDofIter(); it.nextDof())
	if (globalPrimalIndex.count(it.getDofIndex()) == 0) {
478
479
480
481
482
	  recvData = true;
	  break;
	}
	  
      if (recvData)
483
	stdMpi.recv(it.getRank());
484
485
486
487
    }

    stdMpi.startCommunication();

488
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
489
	 !it.end(); it.nextRank()) {
490
      int counter = 0;
491
492
493
494
495
      for (; !it.endDofIter(); it.nextDof())
	if (globalPrimalIndex.count(it.getDofIndex()) == 0)
	  dofFirstLagrange[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[counter++];
    }
496
497
498
499
500
501
502
  }


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

503
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
504
    localIndexB.clear();
505
    DOFAdmin* admin = feSpace->getAdmin();
506
507
508
509

    // === 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.                                 ===
510
511
    
    for (int i = 0; i < admin->getUsedSize(); i++)
512
513
      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
514
	  localIndexB[i] = -1;
515

516
517
518

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

519
    nLocalInterior = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
522
    for (DofMapping::iterator it = localIndexB.begin(); 
	 it != localIndexB.end(); ++it) {
      it->second = nLocalInterior;
523
      nLocalInterior++;
524
    }
525
    nLocalDuals = duals.size();
526

527
    TEST_EXIT_DBG(nLocalInterior + globalPrimalIndex.size() + duals.size() == 
528
529
530
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

531
532
533

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

534
535
    for (DofIndexSet::iterator it = duals.begin();
	 it != duals.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
536
      localIndexB[*it] = globalDualIndex[*it] - rStartB;
537
538
539
  }


540
  void PetscSolverFeti::createMatLagrange()
541
542
543
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

544
545
    // === Create distributed matrix for Lagrange constraints. ===

546
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
547
548
549
550
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
551
552
553
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

554
555
556
557
558
559
560
    // === 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.                                                     ===

561
562
563
564
    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");

565
      // Global index of the first Lagrange constriant for this node.
566
      int index = dofFirstLagrange[*it];
567
      // Copy set of all ranks that contain this dual node.
568
      vector<int> W(boundaryDofRanks[*it].begin(), boundaryDofRanks[*it].end());
569
      // Number of ranks that contain this dual node.
570
571
572
573
574
575
576
      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++) {
577
578
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
	    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);
  }


598
  void PetscSolverFeti::createSchurPrimalKsp()
599
600
601
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
      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 {
630
631
      MSG("Create direct schur primal solver!\n");

632
633
634
635
636
637
      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
638
639
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
640

Thomas Witkowski's avatar
Thomas Witkowski committed
641
642
643
644
645
646
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
647

648
649
650
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
651

652
653
654
655
656
657
658
659
660
661
662
      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
663
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
      }

      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
685
	PetscScalar *tmpValues;
686
687
	VecGetArray(tmpVec, &tmpValues);
	for (int i  = 0; i < nRankB * nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
690
691
692
	  MatSetValue(matBPi, 
		      rStartB * nComponents + i,
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
693
694
695
696
697
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
698
699
700
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
701
702
703
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
704
      MatDestroy(&matBPi);
705
706
707
708
709
710
711
712
713
714

      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
715

716
717
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
718
    }
719
720
721
722
723
724
725
  }


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

726
727
728
729
730
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
731

732
733
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
734

735
736
737
738
739
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
740
741
742
743
744
745
746
  }


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

747
748
    // === Create FETI-DP solver object. ===

749
750
751
    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
752
    fetiData.fetiSolver = this;
753
    fetiData.ksp_schur_primal = &ksp_schur_primal;
754

755
756
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
759
760
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
		 nRankLagrange * nComponents, nOverallLagrange * nComponents,
		 &(fetiData.tmp_vec_lagrange));
761
762
763
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankPrimals * nComponents, nOverallPrimals * nComponents,
		 &(fetiData.tmp_vec_primal));
764
765

    MatCreateShell(PETSC_COMM_WORLD,
766
767
		   nRankLagrange * nComponents, nRankLagrange * nComponents,
		   nOverallLagrange * nComponents, nOverallLagrange * nComponents,
768
		   &fetiData, &mat_feti);
769
770
771
772
773
774
775
    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);
776
777


778
    // === Create FETI-DP preconditioner object. ===
779

780
781
782
783
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
784

785
786
787
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
790
791
792
793
794
795
796
797
798
799
      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;
      
800
801
802
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiDirichletPreconData.tmp_vec_b));      
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
      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;

818
819
820
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiLumpedPreconData.tmp_vec_b));
821
822
823
824
825
826
827
828
      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);
      
829
830
      break;
    default:
831
832
      break;
    }
833
834
835
836
837
838
839
  }
  

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

840
841
    // === Destroy FETI-DP solver object. ===

842
843
844
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
845
    fetiData.fetiSolver = NULL;
846
    fetiData.ksp_schur_primal = PETSC_NULL;
847

Thomas Witkowski's avatar
Thomas Witkowski committed
848
849
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
850
    VecDestroy(&fetiData.tmp_vec_primal);
851
852
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
853
854


855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
    // === 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;
883
884
    default:
      break;
885
    }
886
887
888
889
890
891
892
893
894
  }


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

895
    // === Get local part of the solution for B variables. ===
896
897
898
899
900

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


901
902
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
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
933
934
935
936
937
938
939
940
941
    
    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);

942
943
944
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
945
946
947
948

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

949
    // === And copy from PETSc local vectors to the DOF vectors. ===
950
951
952
953

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

Thomas Witkowski's avatar
Thomas Witkowski committed
954
955
956
      for (DofMapping::iterator it = localIndexB.begin();
	   it != localIndexB.end(); ++it) {
	int petscIndex = it->second * nComponents + i;
957
958
959
960
961
962
963
964
965
966
967
968
969
	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);
970
    VecDestroy(&local_sol_primal);
971
972
973
  }


974
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
975
976
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
977

978
    nComponents = mat->getSize();
979
980
981

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

982
983
    updateDofData();

984
985
986
987
988
989
    // === Create matrices for the FETI-DP method. ===

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
990
    int nRowsInterior = nLocalInterior * nComponents;
991
    int nRowsDual = nLocalDuals * nComponents;
992

Thomas Witkowski's avatar
Thomas Witkowski committed
993
994
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
995
996

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
997
998
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
999
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
1000
1001

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1002
1003
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
1004
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
1005
1006

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1007
1008
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
1009
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
1010

1011

1012
1013
1014
1015
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1016
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
1017
1018
1019
1020
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1021
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
1022
1023
1024
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1025
		      nRowsInterior, nRowsDual, 30, PETSC_NULL,
1026
1027
1028
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1029
		      nRowsDual, nRowsInterior, 30, PETSC_NULL,
1030
1031
		      &mat_duals_interior);
    }
1032

1033
1034
    
    // === Prepare traverse of sequentially created matrices. ===
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049

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

1050
1051
1052
1053
1054
1055
1056
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071

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

1073
	  bool rowPrimal = globalPrimalIndex.count(*cursor) != 0;
1074
  
1075
1076
	  cols.clear();
	  colsOther.clear