PetscSolverFeti.cc 45 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
91
92
    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];
93
94

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


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

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

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


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

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

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

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

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


    // === K_DD ===

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

162

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

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

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

174
175

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

    return 0;
  }


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

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

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


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

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
218
   
219
220
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
221

222
223
      createPrimals(feSpace);      
      createDuals(feSpace);
224
      createLagrange(feSpace);      
225
226
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291

    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();

    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",
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);

      MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
    }

    // === Communicate lagrangeMap to all other ranks. ===

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

      StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
      
      for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	   !it.end(); it.nextRank()) {
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex())) {
	    DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()].local;
	    stdMpi.getSendData(it.getRank()).push_back(d);
	  }
      }
      
      stdMpi.updateSendDataSize();
      
      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	   !it.end(); it.nextRank()) {
	bool recvData = false;
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex())) {
	    recvData = true;
	    break;
	  }
	
	if (recvData)
	  stdMpi.recv(it.getRank());
      }
      
      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	   !it.end(); it.nextRank()) {
	int counter = 0;
	for (; !it.endDofIter(); it.nextDof()) {
	  if (!isPrimal(feSpace, it.getDofIndex())) {
	    DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
	    lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
	  }
	}
      }
    }
292
293
294
  }


295
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
296
 {
297
    FUNCNAME("PetscSolverFeti::createPrimals()");  
298

299
300
301
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

302
303
304
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
305
    DofContainerSet& vertices = 
306
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
307
308
309
310
    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);
311

312
313
314
315

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

316
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
317
318
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
319
320
      else
  	primalDofMap[feSpace].insert(*it);
321

Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
322
323
324
    primalDofMap[feSpace].setOverlap(true);
    primalDofMap[feSpace].setDofComm(meshDistributor->getSendDofs(),
  				     meshDistributor->getRecvDofs());
325

Thomas Witkowski's avatar
Thomas Witkowski committed
326
    primalDofMap[feSpace].update();
327
328
329
  }


330
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
331
332
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
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 (!isPrimal(feSpace, it.getDofIndex())) {
343
344
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
345
346
	}
      }
347
	
348
349
350
351

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

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

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
357
	if (!isPrimal(feSpace, it.getDofIndex()))
358
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
359
360
361

    stdMpi.updateSendDataSize();

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

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

376
377
    stdMpi.startCommunication();

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

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

    DofContainer allBoundaryDofs;
390
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
391
392

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

Thomas Witkowski's avatar
Thomas Witkowski committed
397
    dualDofMap[feSpace].update();
398
399
400
  }

  
401
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
402
403
404
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

405
406
407
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

408
    int nRankLagrange = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
409
410
    map<DegreeOfFreedom, MultiIndex>& dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
411
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
412
	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
413
	int degree = boundaryDofRanks[it->first].size();
414
415
416
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
417
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
418
    lagrangeMap[feSpace].update();
419
420
421
  }


422
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
423
  {
424
    FUNCNAME("PetscSolverFeti::createIndexB()");
425

426
    DOFAdmin* admin = feSpace->getAdmin();
427
428
429
430

    // === 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.                                 ===
431
432
433

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
434
      if (admin->isDofFree(i) == false && 
435
	  isPrimal(feSpace, i) == false &&
436
	  dualDofMap[feSpace].isSet(i) == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
437
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
438
439
	nLocalInterior++;
      }
440
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
441
    
Thomas Witkowski's avatar
Thomas Witkowski committed
442
    TEST_EXIT_DBG(nLocalInterior + primalDofMap[feSpace].size() + 
443
		  dualDofMap[feSpace].size() == 
444
445
446
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualDofMap[feSpace].getMap().begin();
450
	 it != dualDofMap[feSpace].getMap().end(); ++it)
451
452
      localDofMap[feSpace].insertRankDof(it->first);

Thomas Witkowski's avatar
Thomas Witkowski committed
453
    localDofMap[feSpace].update();
454
455
456
  }


457
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
458
459
460
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

461
462
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

463
464
    // === Create distributed matrix for Lagrange constraints. ===

465
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
466
467
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
468
469
470
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

471
472
473
474
475
476
477
    // === 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
478
479
    map<DegreeOfFreedom, MultiIndex> &dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
480
481
482
483
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");

      // Global index of the first Lagrange constriant for this node.
Thomas Witkowski's avatar
Thomas Witkowski committed
484
      int index = lagrangeMap[feSpace][it->first].local;
485
486
487
488
489
490
491
492
493
494
495
496
      // Copy set of all ranks that contain this dual node.
      vector<int> W(boundaryDofRanks[it->first].begin(), 
		    boundaryDofRanks[it->first].end());
      // 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) {
	    // Set the constraint for all components of the system.
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
Thomas Witkowski's avatar
Thomas Witkowski committed
497
498
499
500
	      //	      int colIndex = it->second * nComponents + k;
	      int colIndex = 
		(localDofMap[feSpace].rStartDofs +
		 localDofMap[feSpace][it->first].local) * nComponents + k;
501
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
502
503
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
504
505
	    }
	  }
506
507

	  index++;
508
509
510
511
512
513
514
515
516
	}
      }
    }

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


517
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
518
519
520
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

521
522
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
526
527
528
529
530
531
      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, 
532
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
533
534
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
535
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
536
537
538
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
539
540
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
541
542
543
544
545
546
547
548
549
550
		     &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 {
551
552
      MSG("Create direct schur primal solver!\n");

553
554
555
556
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

557
558
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
559
560
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
561

Thomas Witkowski's avatar
Thomas Witkowski committed
562
563
564
565
566
567
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
568

569
570
571
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
572

573
574
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

575
576
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
577
578
579
580
581
582
	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
583
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
584
585
586
587
588
589
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
590
591
		primalDofMap[feSpace].size() * nComponents)
	("Should not happen!\n");
592
593
594
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
595
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
596
597
598
599
600
601
602
603
604
605

 	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
606
	PetscScalar *tmpValues;
607
	VecGetArray(tmpVec, &tmpValues);
608
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	  MatSetValue(matBPi, 
610
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
613
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
614
615
616
617
618
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
619
620
621
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
622
623
624
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
625
      MatDestroy(&matBPi);
626
627
628
629
630
631
632
633
634
635

      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
636

637
638
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
639
    }
640
641
642
643
644
645
646
  }


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

647
648
649
650
651
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
652

653
654
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
655

656
657
658
659
660
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
661
662
663
  }


664
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
665
666
667
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

668
669
    // === Create FETI-DP solver object. ===

670
671
672
    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
673
    fetiData.fetiSolver = this;
674
    fetiData.ksp_schur_primal = &ksp_schur_primal;
675

676
    VecCreateMPI(PETSC_COMM_WORLD, 
677
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
678
679
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
680
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
681
		 &(fetiData.tmp_vec_lagrange));
682
    VecCreateMPI(PETSC_COMM_WORLD, 
683
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
684
		 &(fetiData.tmp_vec_primal));
685
686

    MatCreateShell(PETSC_COMM_WORLD,
687
688
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
689
		   &fetiData, &mat_feti);
690
691
692
693
694
695
696
    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);
697
698


699
    // === Create FETI-DP preconditioner object. ===
700

701
702
703
704
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
705

706
707
708
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
709
710
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
711
712
713
714
715
716
717
718
719
720
      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;
      
721
      VecCreateMPI(PETSC_COMM_WORLD, 
722
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
723
		   &(fetiDirichletPreconData.tmp_vec_b));      
724
725
726
727
728
729
      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));
730
731
732
733
734
735
736
737
738
739
740
741
      
      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;

742
      VecCreateMPI(PETSC_COMM_WORLD, 
743
744
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
745
		   &(fetiLumpedPreconData.tmp_vec_b));
746
747
748
749
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
750
751
752
753
754
755

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
756
757
      break;
    default:
758
759
      break;
    }
760
761
762
763
764
765
766
  }
  

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

767
768
    // === Destroy FETI-DP solver object. ===

769
770
771
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
772
    fetiData.fetiSolver = NULL;
773
    fetiData.ksp_schur_primal = PETSC_NULL;
774

Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
777
    VecDestroy(&fetiData.tmp_vec_primal);
778
779
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
780
781


782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
    // === 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;
810
811
    default:
      break;
812
    }
813
814
815
816
817
818
819
820
821
  }


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

822
823
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

824
    // === Get local part of the solution for B variables. ===
825
826
827
828
829

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


830
831
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
832
833
    
    vector<PetscInt> globalIsIndex, localIsIndex;
834
835
    globalIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
    localIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
836
837
838

    {
      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
839
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
840
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
841
	for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
842
	  globalIsIndex.push_back(it->second.local * nComponents + i);
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
	  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);

871
872
873
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
874
875
876
877

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

878
    // === And copy from PETSc local vectors to the DOF vectors. ===
879
880
881
882

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

Thomas Witkowski's avatar
Thomas Witkowski committed
883
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpace].getMap().begin();
884
	   it != localDofMap[feSpace].getMap().end(); ++it) {
885
	int petscIndex = localDofMap.mapLocal(it->first, i);
886
887
888
889
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
890
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
891
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
892
893
894
895
896
897
898
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
899
    VecDestroy(&local_sol_primal);
900
901
902
  }


903
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
904
905
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
906

907
908
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);

909
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
910
    nComponents = mat->getSize();
911
912
913

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

Thomas Witkowski's avatar
Thomas Witkowski committed
914
915
916
917
    primalDofMap.init(mpiComm, feSpaces, true);
    dualDofMap.init(mpiComm, feSpaces, false);
    lagrangeMap.init(mpiComm, feSpaces, true);
    localDofMap.init(mpiComm, feSpaces, false);
918

919
920
    updateDofData();

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

923
924
    int nRowsRankB = localDofMap.getRankDofs();
    int nRowsOverallB = localDofMap.getOverallDofs();
925
926
927
    int nRowsRankPrimal = primalDofMap.getRankDofs();
    int nRowsOverallPrimal = primalDofMap.getOverallDofs();
    int nRowsDual = dualDofMap.getRankDofs();
928
929
    int nRowsInterior = nRowsRankB - nRowsDual;

Thomas Witkowski's avatar
Thomas Witkowski committed
930
931
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
932
933

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
934
935
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
936
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
937
938

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
939
940
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
941
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
942
943

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
944
945
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
946
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
947

948

949
950
951
952
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
953
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
954
955
956
957
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
958
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
959
960
961
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
962
		      nRowsInterior, nRowsDual, 30, PETSC_NULL,
963
964
965
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
966
		      nRowsDual, nRowsInterior, 30, PETSC_NULL,
967
968
		      &mat_duals_interior);
    }
969

970
971
    
    // === Prepare traverse of sequentially created matrices. ===
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986

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

987
988
989
990
991
992
993
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008

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

1010
	  bool rowPrimal = primalDofMap[feSpace].isSet(*cursor);
1011
  
1012
1013
	  cols.clear();
	  colsOther.clear();