PetscSolverFeti.cc 53.4 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
    for (it_type it = meshDistributor->getSendDofs().begin();
	 it != meshDistributor->getSendDofs().end(); ++it)
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt)
290
291
292
293
	if (globalPrimalIndex.count(**dofIt))
	  stdMpi.getSendData(it->first).push_back(globalPrimalIndex[**dofIt]);
    stdMpi.updateSendDataSize();

294
295
    for (it_type it = meshDistributor->getRecvDofs().begin();
	 it != meshDistributor->getRecvDofs().end(); ++it) {
296
      bool recvFromRank = false;
297
298
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt)
299
	if (primals.count(**dofIt) && 
300
	    meshDistributor->getIsRankDof(feSpace, **dofIt) == false) {
301
302
303
304
305
306
307
308
309
	  recvFromRank = true;
	  break;
	}

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

310
311
    for (it_type it = meshDistributor->getRecvDofs().begin();
	 it != meshDistributor->getRecvDofs().end(); ++it) {
312
      int i = 0;
313
314
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt) {
315
	if (primals.count(**dofIt) && 
316
	    meshDistributor->getIsRankDof(feSpace, **dofIt) == false)
317
318
319
320
321
322
323
324
325
326
327
328
329
	  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()");
330
331

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

    boundaryDofRanks.clear();

338
339
340
341
342
343
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

    for (it_type it = meshDistributor->getSendDofs().begin();
	 it != meshDistributor->getSendDofs().end(); ++it) {
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt) {
344
	// If DOF is not primal, i.e., its a dual node
345
	if (globalPrimalIndex.count(**dofIt) == 0) {
346
347
348
349
350
351
	  boundaryDofRanks[**dofIt].insert(mpiRank);
	  boundaryDofRanks[**dofIt].insert(it->first);
	}
      }
    }

352
353
354
355

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

356
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
357
358
359
360
    for (it_type it = meshDistributor->getSendDofs().begin();
	 it != meshDistributor->getSendDofs().end(); ++it)
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt)
361
	if (globalPrimalIndex.count(**dofIt) == 0)
362
363
364
365
	  stdMpi.getSendData(it->first).push_back(boundaryDofRanks[**dofIt]);

    stdMpi.updateSendDataSize();

366
367
    for (it_type it = meshDistributor->getRecvDofs().begin();
	 it != meshDistributor->getRecvDofs().end(); ++it) {
368
      bool recvFromRank = false;
369
370
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt)
371
	if (globalPrimalIndex.count(**dofIt) == 0) {
372
373
374
375
376
377
378
379
380
	  recvFromRank = true;
	  break;
	}

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

381
382
    for (it_type it = meshDistributor->getRecvDofs().begin();
	 it != meshDistributor->getRecvDofs().end(); ++it) {
383
      int i = 0;
384
385
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt)	
386
	if (globalPrimalIndex.count(**dofIt) == 0)
387
	  boundaryDofRanks[**dofIt] = stdMpi.getRecvData(it->first)[i++];
388
389
390
391
392
393
394
395
    }


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

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

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

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

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

419
420
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
	nRankDuals, nOverallDuals);
421
422
423
424
425
426
427
  }

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

428
429
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

430
431
432
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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

442
443
444
445
446

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

447
    nOverallLagrange = 0;
448
    rStartLagrange = 0;
449
450
451
452
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankLagrange, rStartLagrange, nOverallLagrange);

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

456
457
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	nRankLagrange, nOverallLagrange);
458
459


460
    // === Communicate dofFirstLagrange to all other ranks. ===
461
462

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
463
464
465
466
467
468
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

    for (it_type it = meshDistributor->getSendDofs().begin();
	 it != meshDistributor->getSendDofs().end(); ++it)
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt) {
469
	if (globalPrimalIndex.count(**dofIt) == 0) {
470
471
472
473
474
475
	  TEST_EXIT_DBG(dofFirstLagrange.count(**dofIt))("Should not happen!\n");
	  stdMpi.getSendData(it->first).push_back(dofFirstLagrange[**dofIt]);
	}
      }
    stdMpi.updateSendDataSize();

476
477
    for (it_type it = meshDistributor->getRecvDofs().begin();
	 it != meshDistributor->getRecvDofs().end(); ++it) {
478
      bool recvData = false;
479
480
      for (DofContainer::iterator dofIt = it->second[feSpace].begin();
	   dofIt != it->second[feSpace].end(); ++dofIt)
481
	if (globalPrimalIndex.count(**dofIt) == 0) {
482
483
484
485
486
487
488
489
490
491
	  recvData = true;
	  break;
	}
	  
      if (recvData)
	stdMpi.recv(it->first);
    }

    stdMpi.startCommunication();

492
493
    for (it_type it = meshDistributor->getRecvDofs().begin();
	 it != meshDistributor->getRecvDofs().end(); ++it) {
494
      int counter = 0;
495
496
497
498
      DofContainer &dc = it->second[feSpace];
      for (unsigned int i = 0; i < dc.size(); i++)
	if (globalPrimalIndex.count(*(dc[i])) == 0)
	  dofFirstLagrange[*(dc[i])] = stdMpi.getRecvData(it->first)[counter++];
499
    }     
500
501
502
503
504
505
506
  }


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

507
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
508
    localIndexB.clear();
509
    DOFAdmin* admin = feSpace->getAdmin();
510
511
512
513

    // === 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.                                 ===
514
515
    
    for (int i = 0; i < admin->getUsedSize(); i++)
516
517
      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
518
	  localIndexB[i] = -1;
519

520
521
522

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

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

531
    TEST_EXIT_DBG(nLocalInterior + globalPrimalIndex.size() + duals.size() == 
532
533
534
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

535
536
537

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

538
539
    for (DofIndexSet::iterator it = duals.begin();
	 it != duals.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
540
      localIndexB[*it] = globalDualIndex[*it] - rStartB;
541
542
543
  }


544
  void PetscSolverFeti::createMatLagrange()
545
546
547
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

548
549
    // === Create distributed matrix for Lagrange constraints. ===

550
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
551
552
553
554
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
555
556
557
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

558
559
560
561
562
563
564
    // === 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.                                                     ===

565
566
567
568
    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");

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


602
  void PetscSolverFeti::createSchurPrimalKsp()
603
604
605
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

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

636
637
638
639
640
641
      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
642
643
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
644

Thomas Witkowski's avatar
Thomas Witkowski committed
645
646
647
648
649
650
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
651

652
653
654
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
655

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

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
702
703
704
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
705
706
707
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
708
      MatDestroy(&matBPi);
709
710
711
712
713
714
715
716
717
718

      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
719

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


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

730
731
732
733
734
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
735

736
737
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
738

739
740
741
742
743
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
744
745
746
747
748
749
750
  }


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

751
752
    // === Create FETI-DP solver object. ===

753
754
755
    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
756
    fetiData.fetiSolver = this;
757
    fetiData.ksp_schur_primal = &ksp_schur_primal;
758

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

    MatCreateShell(PETSC_COMM_WORLD,
770
771
		   nRankLagrange * nComponents, nRankLagrange * nComponents,
		   nOverallLagrange * nComponents, nOverallLagrange * nComponents,
772
		   &fetiData, &mat_feti);
773
774
775
776
777
778
779
    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);
780
781


782
    // === Create FETI-DP preconditioner object. ===
783

784
785
786
787
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
788

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

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

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

844
845
    // === Destroy FETI-DP solver object. ===

846
847
848
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
849
    fetiData.fetiSolver = NULL;
850
    fetiData.ksp_schur_primal = PETSC_NULL;
851

Thomas Witkowski's avatar
Thomas Witkowski committed
852
853
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
854
    VecDestroy(&fetiData.tmp_vec_primal);
855
856
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
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
883
884
885
886
    // === 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;
887
888
    default:
      break;
889
    }
890
891
892
893
894
895
896
897
898
  }


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

899
    // === Get local part of the solution for B variables. ===
900
901
902
903
904

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


905
906
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
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
942
943
944
945
    
    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);

946
947
948
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
949
950
951
952

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

953
    // === And copy from PETSc local vectors to the DOF vectors. ===
954
955
956
957

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

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


978
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
979
980
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
981

982
    nComponents = mat->getSize();
983
984
985

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

986
987
    updateDofData();

988
989
990
991
992
993
994

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

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
995
    int nRowsInterior = nLocalInterior * nComponents;
996
    int nRowsDual = nLocalDuals * nComponents;
997

Thomas Witkowski's avatar
Thomas Witkowski committed
998
999
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
1000
1001

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

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1007
1008
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
1009
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
1010
1011

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1012
1013
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
1014
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
1015

1016

1017
1018
1019
1020
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1021
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
1022
1023
1024
1025
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1026
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
1027
1028
1029
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1030
		      nRowsInterior, nRowsDual, 30, PETSC_NULL,
1031
1032
1033
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1034
		      nRowsDual, nRowsInterior, 30, PETSC_NULL,
1035
1036
		      &mat_duals_interior);
    }
1037

1038
1039
    
    // === Prepare traverse of sequentially created matrices. ===
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054

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

1055
1056
1057
1058
1059
1060
1061
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076

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

1078
	  bool rowPrimal = globalPrimalIndex.count(*cursor) != 0;
1079
  
1080
1081
	  cols.clear();
	  colsOther.clear();
1082
	  values.clear();	  
1083
	  valuesOther.clear();
1084
1085
1086
1087
1088
1089

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

1090
1091
1092
1093
1094
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

1095
	    bool colPrimal = globalPrimalIndex.count(col(*icursor)) != 0;
1096
1097

	    if (colPrimal) {
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
	      // 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));
1108
	      } else {
1109
1110
1111
1112
1113
1114
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      // Column is not a primal variable.

Thomas Witkowski's avatar
Thomas Witkowski committed
1115
	      TEST_EXIT_DBG(localIndexB.count(col(*icursor)))
1116
		("No global B index for DOF %d!\n", col(*icursor));
1117
	      
Thomas Witkowski's avatar
Thomas Witkowski committed
1118
1119
	      int colIndex = 
		(localIndexB[col(*icursor)] + rStartB) * nComponents + j;
1120
1121
1122
1123
1124
1125
1126

	      if (rowPrimal) {
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(colIndex);
		values.push_back(value(*icursor));
1127
1128
	      }
	    }
1129
1130
1131
1132
1133
1134



	    // === For preconditioner ===

	    if (!rowPrimal && !colPrimal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1135
1136
	      int rowIndex = localIndexB[*cursor];
	      int colIndex = localIndexB[col(*icursor)];
1137
1138
1139
		
	      if (rowIndex < nLocalInterior) {
		if (colIndex < nLocalInterior) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1140
		  int colIndex2 = localIndexB[col(*icursor)] * nComponents + j;
1141
1142
1143
1144

		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1145
1146
		  int colIndex2 = (localIndexB[col(*icursor)] - nLocalInterior) *
		    nComponents + j;
1147
1148
1149
1150
1151
1152

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		}
	      } else {
		if (colIndex < nLocalInterior) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1153
		  int colIndex2 = localIndexB[col(*icursor)] * nComponents + j;
1154
1155
1156
1157

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1158
1159
		  int colIndex2 = (localIndexB[col(*icursor)] - nLocalInterior) *
		    nComponents + j;
1160
1161
1162
1163
1164
1165
1166
1167

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


1168
	  }  // for each nnz in row
1169

1170
1171
1172
	  if (rowPrimal) {
	    TEST_EXIT_DBG(globalPrimalIndex.count(*cursor))
	      ("Should not happen!\n");
1173