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
223
224
    createPrimals();

    createDuals();

    createLagrange();

    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
462
463
464
465
466
467
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      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
    stdMpi.updateSendDataSize();

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

    stdMpi.startCommunication();

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


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

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

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

515
516
517

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

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

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

530
531
532

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

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


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

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

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

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

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

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


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

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

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

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

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

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

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

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

	VecDestroy(&tmpVec);
      }

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

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

      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
714

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


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

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

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

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


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

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

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

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

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


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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
847
848
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
849
    VecDestroy(&fetiData.tmp_vec_primal);
850
851
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
852
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
    // === 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;
882
883
    default:
      break;
884
    }
885
886
887
888
889
890
891
892
893
  }


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

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

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


900
901
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
902
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
    
    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);

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

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

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

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

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


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

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

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

981
982
    updateDofData();

983
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();
1077
	  values.clear();	  
1078
	  valuesOther.clear();
1079
1080
1081
1082
1083
1084

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

1085
1086
1087
1088
1089
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

1090
	    bool colPrimal = globalPrimalIndex.count(col(*icursor)) != 0;