PetscSolverFeti.cc 52.8 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
234
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

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

245
246
247
248

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

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

257

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

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

266
267
268

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

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

273
274
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	nRankPrimals, nOverallPrimals);
275

276
277
278
279
280

    // === Communicate primal's global index from ranks that own the     ===
    // === primals to ranks that contain this primals but are not owning ===
    // === them.                                                         ===

281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (globalPrimalIndex.count(**dofIt))
	  stdMpi.getSendData(it->first).push_back(globalPrimalIndex[**dofIt]);
    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvFromRank = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false) {
	  recvFromRank = true;
	  break;
	}

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

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int i = 0;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false)
	  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()");
    
329
330
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
331
332
333
334
335
336
337
338
339

    boundaryDofRanks.clear();

    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it) {
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	// If DOF is not primal, i.e., its a dual node
340
	if (globalPrimalIndex.count(**dofIt) == 0) {
341
342
343
344
345
346
	  boundaryDofRanks[**dofIt].insert(mpiRank);
	  boundaryDofRanks[**dofIt].insert(it->first);
	}
      }
    }

347
348
349
350

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

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

    stdMpi.updateSendDataSize();

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

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

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


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

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

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

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

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

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

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

424
425
426
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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

436
437
438
439
440

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

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

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

450
451
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	nRankLagrange, nOverallLagrange);
452
453


454
    // === Communicate dofFirstLagrange to all other ranks. ===
455
456
457
458
459
460
461

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
462
	if (globalPrimalIndex.count(**dofIt) == 0) {
463
464
465
466
467
468
469
470
471
472
473
474
	  TEST_EXIT_DBG(dofFirstLagrange.count(**dofIt))("Should not happen!\n");
	  stdMpi.getSendData(it->first).push_back(dofFirstLagrange[**dofIt]);
	}
      }
    stdMpi.updateSendDataSize();

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

    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int counter = 0;
      for (unsigned int i = 0; i < it->second.size(); i++)
490
	if (globalPrimalIndex.count(*(it->second[i])) == 0)
491
	  dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++];
492
    }     
493
494
495
496
497
498
499
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
500
    localIndexB.clear();
501
    DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin();
502
503
504
505

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

512
513
514

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

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

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

527
528
529

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

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


536
  void PetscSolverFeti::createMatLagrange()
537
538
539
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

540
541
    // === Create distributed matrix for Lagrange constraints. ===

542
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
543
544
545
546
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
547
548
549
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

550
551
552
553
554
555
556
    // === 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.                                                     ===

557
558
559
560
    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");

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


594
  void PetscSolverFeti::createSchurPrimalKsp()
595
596
597
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

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

628
629
630
631
632
633
      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
634
635
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
636

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

644
645
646
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
647

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

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

	VecDestroy(&tmpVec);
      }

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

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
700
      MatDestroy(&matBPi);
701
702
703
704
705
706
707
708
709
710

      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
711

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


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

722
723
724
725
726
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
727

728
729
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
730

731
732
733
734
735
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
736
737
738
739
740
741
742
  }


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

743
744
    // === Create FETI-DP solver object. ===

745
746
747
    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
748
    fetiData.fetiSolver = this;
749
    fetiData.ksp_schur_primal = &ksp_schur_primal;
750

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

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


774
    // === Create FETI-DP preconditioner object. ===
775

776
777
778
779
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
780

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

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

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

836
837
    // === Destroy FETI-DP solver object. ===

838
839
840
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
841
    fetiData.fetiSolver = NULL;
842
    fetiData.ksp_schur_primal = PETSC_NULL;
843

Thomas Witkowski's avatar
Thomas Witkowski committed
844
845
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
846
    VecDestroy(&fetiData.tmp_vec_primal);
847
848
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
849
850


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


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

891
    // === Get local part of the solution for B variables. ===
892
893
894
895
896

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


897
898
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
    
    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);

938
939
940
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
941
942
943
944

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

945
    // === And copy from PETSc local vectors to the DOF vectors. ===
946
947
948
949

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

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


970
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
971
972
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
973

974
    nComponents = mat->getSize();
975
976
977

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

978
979
    updateDofData();

980
981
982
983
984
985
986

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

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
987
    int nRowsInterior = nLocalInterior * nComponents;
988
    int nRowsDual = nLocalDuals * nComponents;
989

Thomas Witkowski's avatar
Thomas Witkowski committed
990
991
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
992
993

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
994
995
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
996
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
997
998

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
999
1000
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
1001
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
1002
1003

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1004
1005
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
1006
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
1007

1008

1009
1010
1011
1012
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1013
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
1014
1015
1016
1017
		      &mat_duals_duals);

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

1030
1031
    
    // === Prepare traverse of sequentially created matrices. ===
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046

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

1047
1048
1049
1050
1051
1052
1053
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068

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

1070
	  bool rowPrimal = globalPrimalIndex.count(*cursor) != 0;
1071
  
1072
1073
	  cols.clear();
	  colsOther.clear();
1074
	  values.clear();	  
1075
	  valuesOther.clear();
1076
1077
1078
1079
1080
1081

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

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

1087
	    bool colPrimal = globalPrimalIndex.count(col(*icursor)) != 0;
1088
1089

	    if (colPrimal) {
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
	      // 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));
1100
	      } else {
1101
1102
1103
1104
1105
1106
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      // Column is not a primal variable.

Thomas Witkowski's avatar
Thomas Witkowski committed
1107
	      TEST_EXIT_DBG(localIndexB.count(col(*icursor)))
1108
		("No global B index for DOF %d!\n", col(*icursor));
1109
	      
Thomas Witkowski's avatar
Thomas Witkowski committed
1110
1111
	      int colIndex = 
		(localIndexB[col(*icursor)] + rStartB) * nComponents + j;
1112
1113
1114
1115
1116
1117
1118

	      if (rowPrimal) {
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(colIndex);
		values.push_back(value(*icursor));
1119
1120
	      }
	    }
1121
1122
1123
1124
1125
1126



	    // === For preconditioner ===

	    if (!rowPrimal && !colPrimal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1127
1128
	      int rowIndex = localIndexB[*cursor];
	      int colIndex = localIndexB[col(*icursor)];
1129
1130
1131
		
	      if (rowIndex < nLocalInterior) {
		if (colIndex < nLocalInterior) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1132
		  int colIndex2 = localIndexB[col(*icursor)] * nComponents + j;
1133
1134
1135
1136

		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1137
1138
		  int colIndex2 = (localIndexB[col(*icursor)] - nLocalInterior) *
		    nComponents + j;
1139
1140
1141
1142
1143
1144

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		}
	      } else {
		if (colIndex < nLocalInterior) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1145
		  int colIndex2 = localIndexB[col(*icursor)] * nComponents + j;
1146
1147
1148
1149

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1150
1151
		  int colIndex2 = (localIndexB[col(*icursor)] - nLocalInterior) *
		    nComponents + j;
1152
1153
1154
1155
1156
1157
1158
1159

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


1160
	  }  // for each nnz in row
1161

1162
1163
1164
	  if (rowPrimal) {
	    TEST_EXIT_DBG(globalPrimalIndex.count(*cursor))
	      ("Should not happen!\n");
1165

1166
1167
1168
	    int rowIndex = globalPrimalIndex[*cursor] * nComponents + i;
	    MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1169

1170
1171
1172
1173
	    if (colsOther.size())
	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1174
	    TEST_EXIT_DBG(localIndexB.count(*cursor))
1175
	      ("Should not happen!\n");
1176

1177
	    int localRowIndex = localIndexB[*cursor] * nComponents + i; 
Thomas Witkowski's avatar
Thomas Witkowski committed
1178
1179
	    for (unsigned int k = 0; k < cols.size(); k++)
	      cols[k] -= rStartB * nComponents;