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

namespace AMDiS {

  using namespace std;

23
24
25
26
27
28
29
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
30
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
31
32

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

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

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

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

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

    return 0;
  }


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

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


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

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

87
88
89
90
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

91
92
93
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
94
95

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


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

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

103
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
104
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
105
106
107
108
109
110
111
112
113
114
    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);

115
116
117
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139

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


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
140
141


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
153
154
155
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
156
157

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


    // === K_DD ===

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

165

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

168
169
170
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
173
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
174

175
176
    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
177

178
179

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

    return 0;
  }


186
187
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
Thomas Witkowski's avatar
Thomas Witkowski committed
188
      schurPrimalSolver(0)
189
190
191
192
193
194
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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


215
  void PetscSolverFeti::updateDofData()
216
217
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
218
219
220

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
221
   
222
223
224
225
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createPrimals(feSpace);      
      createDuals(feSpace);
226
      createLagrange(feSpace);      
227
228
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
229

230
231
232
233
234
    primalDofMap.setDofComm(meshDistributor->getSendDofs(),
			    meshDistributor->getRecvDofs());
    lagrangeMap.setDofComm(meshDistributor->getSendDofs(), 
			   meshDistributor->getRecvDofs());

Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237
238
    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();
239
240
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
241
242
243
244
245
246
247
248

    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);

      MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
      
      MSG("nRankDuals = %d   nOverallDuals = %d\n",
249
	  dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
250
251

      MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
252
	  lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
253
254
255
256

      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
257
    }
258
259
260
  }


261
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
262
  {
263
    FUNCNAME("PetscSolverFeti::createPrimals()");  
264

265
266
267
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

268
269
270
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
271
    DofContainerSet& vertices = 
272
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
273
274
275
276
    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);
277

278
279
280
281

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

282
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
283
284
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
285
286
      else
  	primalDofMap[feSpace].insert(*it);
287
288
289
  }


290
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
291
292
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
293

294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
	dualDofMap[feSpace].insertRankDof(**it);
  }

  
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

310
311
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
312

313
    boundaryDofRanks[feSpace].clear();
314

315
316
317
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
318
	if (!isPrimal(feSpace, it.getDofIndex())) {
319
320
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
321
322
	}
      }
323
	
324
325
326
327

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

328
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
329
330
331
332

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
333
	if (!isPrimal(feSpace, it.getDofIndex()))
334
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
335
336
337

    stdMpi.updateSendDataSize();

338
339
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
340
      bool recvFromRank = false;
341
      for (; !it.endDofIter(); it.nextDof()) {
342
	if (!isPrimal(feSpace, it.getDofIndex())) {
343
344
345
	  recvFromRank = true;
	  break;
	}
346
      }
347
348

      if (recvFromRank)
349
	stdMpi.recv(it.getRank());
350
    }
351

352
353
    stdMpi.startCommunication();

354
355
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
356
      int i = 0;
357
      for (; !it.endDofIter(); it.nextDof())
358
	if (!isPrimal(feSpace, it.getDofIndex()))
359
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
360
	    stdMpi.getRecvData(it.getRank())[i++];
361
362
363
    }


364
365
366
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

367
    int nRankLagrange = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
368
369
    map<DegreeOfFreedom, MultiIndex>& dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
370
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
371
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
372
	int degree = boundaryDofRanks[feSpace][it->first].size();
373
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
374
375
      } else {
	lagrangeMap[feSpace].insert(it->first);
376
377
      }
    }
378
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
379
380
381
  }


382
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
383
  {
384
    FUNCNAME("PetscSolverFeti::createIndexB()");
385

386
    DOFAdmin* admin = feSpace->getAdmin();
387
388
389
390

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

392
    int nLocalInterior = 0;    
393
    for (int i = 0; i < admin->getUsedSize(); i++) {
394
      if (admin->isDofFree(i) == false && 
395
	  isPrimal(feSpace, i) == false &&
396
	  dualDofMap[feSpace].isSet(i) == false) {
397
398
399
400
401
402
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);

	if (fetiPreconditioner == FETI_DIRICHLET)
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);

	nLocalInterior++;
403
      }
404
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
405
    
406
407
    // === And finally, add the global indicies of all dual nodes. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
408
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualDofMap[feSpace].getMap().begin();
409
	 it != dualDofMap[feSpace].getMap().end(); ++it)
410
      localDofMap[feSpace].insertRankDof(it->first);
411
412
413
  }


414
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
415
416
417
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

418
419
    // === Create distributed matrix for Lagrange constraints. ===

420
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
421
422
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
423
424
425
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

426
427
428
429
430
431
432
    // === 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.                                                     ===

Thomas Witkowski's avatar
Thomas Witkowski committed
433
434
435
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
      map<DegreeOfFreedom, MultiIndex> &dualMap = 
	dualDofMap[feSpaces[k]].getMap();
436

Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	   it != dualMap.end(); ++it) {
439
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
442
443
444
445
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
	int index = lagrangeMap.getMatIndex(k, it->first);
	
	// Copy set of all ranks that contain this dual node.
446
447
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
450
451
452
453
	// Number of ranks that contain this dual node.
	int degree = W.size();
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
454
	      int colIndex = localDofMap.getMatIndex(k, it->first);
455
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
456
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
457
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
458
	    index++;	      
459
460
461
462
463
464
465
466
467
468
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
  }


469
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
470
471
472
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
478
479
480
481
      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, 
482
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
483
484
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
485
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
486
487
488
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
489
490
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
493
494
495
496
497
498
499
500
		     &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 {
501
502
      MSG("Create direct schur primal solver!\n");

503
504
505
506
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

507
508
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
509
510
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
511

Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
514
515
516
517
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
518

519
520
521
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
522

523
524
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

525
526
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
527
528
529
530
531
532
	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
533
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
534
535
536
537
538
      }

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

539
540
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
		primalDofMap.getLocalDofs())
541
	("Should not happen!\n");
542

543
544
545
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
546
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
547
548
549
550
551
552
553
554
555
556

 	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
557
	PetscScalar *tmpValues;
558
	VecGetArray(tmpVec, &tmpValues);
559
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
560
	  MatSetValue(matBPi, 
561
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
562
563
564
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
565
566
567
568
569
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
570
571
572
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
573
574
575
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
576
      MatDestroy(&matBPi);
577
578
579
580
581
582
583
584
585
586

      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
587

588
589
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
590
    }
591
592
593
594
595
596
597
  }


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

598
599
600
601
602
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
603

604
605
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
606

607
608
609
610
611
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
612
613
614
  }


615
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
616
617
618
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

619
620
    // === Create FETI-DP solver object. ===

621
622
623
    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
624
    fetiData.fetiSolver = this;
625
    fetiData.ksp_schur_primal = &ksp_schur_primal;
626

627
    VecCreateMPI(PETSC_COMM_WORLD, 
628
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
631
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
632
		 &(fetiData.tmp_vec_lagrange));
633
    VecCreateMPI(PETSC_COMM_WORLD, 
634
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
635
		 &(fetiData.tmp_vec_primal));
636
637

    MatCreateShell(PETSC_COMM_WORLD,
638
639
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
640
		   &fetiData, &mat_feti);
641
642
643
644
645
646
647
    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);
648
649


650
    // === Create FETI-DP preconditioner object. ===
651

652
653
654
655
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
656

657
658
659
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
660
661
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
662
663
664
665
666
667
668
669
670
671
      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;
      
672
      VecCreateMPI(PETSC_COMM_WORLD, 
673
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
674
		   &(fetiDirichletPreconData.tmp_vec_b));      
675
676
677
678
679
680
      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));
681
682
683
684
685
686
687
688
689
690
691
692
693

      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex> &dualMap = 
	  dualDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	     it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

694
695
696
697
698
699
700
701
702
703
704
      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;

Thomas Witkowski's avatar
Thomas Witkowski committed
705
706
707
708
709
710
711
712
713
714
715
716
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex> &dualMap = 
	  dualDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	     it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

717
      VecCreateMPI(PETSC_COMM_WORLD, 
718
719
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
720
		   &(fetiLumpedPreconData.tmp_vec_b));
721
722
723
724
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
725
726
727
728
729
730

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
731
732
      break;
    default:
733
734
      break;
    }
735
736
737
738
739
740
741
  }
  

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

742
743
    // === Destroy FETI-DP solver object. ===

744
745
746
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
747
    fetiData.fetiSolver = NULL;
748
    fetiData.ksp_schur_primal = PETSC_NULL;
749

Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
752
    VecDestroy(&fetiData.tmp_vec_primal);
753
754
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
755
756


757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
    // === 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;
785
786
    default:
      break;
787
    }
788
789
790
791
792
793
794
795
796
  }


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

797
    // === Get local part of the solution for B variables. ===
798
799
800
801
802

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


803
804
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
805
    
Thomas Witkowski's avatar
Thomas Witkowski committed
806
807
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);

808
    vector<PetscInt> globalIsIndex, localIsIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
809
810
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());
811
812

    {
Thomas Witkowski's avatar
Thomas Witkowski committed
813
814
815
816
817
818
819
820
      int cnt = 0;
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex>& dofMap = 
	  primalDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin();
	     it != dofMap.end(); ++it) {
	  globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
	  localIsIndex.push_back(cnt++);
821
822
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
823
824

      TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs())("Should not happen!\n");
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
    }
    
    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);

850
851
852
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
853
854
855
856

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

857
    // === And copy from PETSc local vectors to the DOF vectors. ===
858

Thomas Witkowski's avatar
Thomas Witkowski committed
859
860
    int cnt = 0;
    for (int i = 0; i < vec.getSize(); i++) {
861
862
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));

Thomas Witkowski's avatar
Thomas Witkowski committed
863
864
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpaces[i]].getMap().begin();
	   it != localDofMap[feSpaces[i]].getMap().end(); ++it) {
865
	int petscIndex = localDofMap.getLocalMatIndex(i, it->first);
866
867
868
	dofVec[it->first] = localSolB[petscIndex];
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
869
870
871
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpaces[i]].getMap().begin();
	   it != primalDofMap[feSpaces[i]].getMap().end(); ++it)
	dofVec[it->first] = localSolPrimal[cnt++];
872
873
874
875
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
876
    VecDestroy(&local_sol_primal);
877
878
879
  }


880
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
881
882
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
883
    
884
    // === Create all sets and indices. ===
885
    
Thomas Witkowski's avatar
Thomas Witkowski committed
886
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
887
    
888
889
890
891
    primalDofMap.init(mpiComm, feSpaces, true, true);
    dualDofMap.init(mpiComm, feSpaces, false, false);
    lagrangeMap.init(mpiComm, feSpaces, true, true);
    localDofMap.init(mpiComm, feSpaces, false, false);
892
893
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.init(mpiComm, feSpaces, false, false);
894

895
    updateDofData();
896
    
897
    // === Create matrices for the FETI-DP method. ===
898
    
899
900
    int nRowsRankB = localDofMap.getRankDofs();
    int nRowsOverallB = localDofMap.getOverallDofs();
901
902
    int nRowsRankPrimal = primalDofMap.getRankDofs();
    int nRowsOverallPrimal = primalDofMap.getOverallDofs();
903
    
Thomas Witkowski's avatar
Thomas Witkowski committed
904
905
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
906
    
907
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
908
909
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
910
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
911
    
912
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
913
914
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
915
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
916
    
917
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
918
919
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
920
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
921
922
    
    
923
    // === Create matrices for FETI-DP preconditioner. ===
924
925
926
    
    if (fetiPreconditioner != FETI_NONE) {
      int nRowsDual = dualDofMap.getRankDofs();
927
928

      MatCreateSeqAIJ(PETSC_COMM_SELF,
929
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
930
931
		      &mat_duals_duals);

932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
      if (fetiPreconditioner == FETI_DIRICHLET) {
	int nRowsInterior = interiorDofMap.getRankDofs();

	MatCreateSeqAIJ(PETSC_COMM_SELF,
			nRowsInterior, nRowsInterior, 30, PETSC_NULL,
			&mat_interior_interior);
	
	MatCreateSeqAIJ(PETSC_COMM_SELF,
			nRowsInterior, nRowsDual, 30, PETSC_NULL,
			&mat_interior_duals);
	
	MatCreateSeqAIJ(PETSC_COMM_SELF,
			nRowsDual, nRowsInterior, 30, PETSC_NULL,
			&mat_duals_interior);
      }
947
    }
948

949
950
    
    // === Prepare traverse of sequentially created matrices. ===
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965

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

966
967
968
969
970
971
972
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

973
974
975
976

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

Thomas Witkowski's avatar
Thomas Witkowski committed
977
    int nComponents = mat->getSize();
978
979
980
981
982
983
984
985
986
987
988
    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