PetscSolverFeti.cc 45.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
//
// 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"

namespace AMDiS {

  using namespace std;


#ifdef HAVE_PETSC_DEV 
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
33
34
35
36
37
38
39
40
41
42
43
44
45

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);

    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)
  {
46
47
    //    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)
48
49
50

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

53
54
55
    // tmp_vec_b0 = inv(K_BB) trans(L) x
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b0, data->tmp_vec_b0);
56

57
58
    // tmp_vec_b1 = inv(K_BB) K_BPi  inv(S_PiPi) K_PiB tmp_vec_b0
    MatMult(*(data->mat_primal_b), data->tmp_vec_b0, data->tmp_vec_primal);
59
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
60
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b1);
61

62
63
64
    // y = L (tmp_vec_b0 + tmp_vec_b1)
    VecAXPBY(data->tmp_vec_b0, 1.0, 1.0, data->tmp_vec_b1);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y);
65
66
67
68
69

    return 0;
  }


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

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


82
    // === Restriction of the B nodes to the boundary nodes. ===
83

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

89
90
91
92
93
94
    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];
95
96

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


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

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

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

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

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


    // === K_DD ===

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

164

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

167
168
169
170
171
    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];
172

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

176
177

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

    return 0;
  }


184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
      nComponents(-1)
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
      fetiPreconditioner = FETI_LUMPED;
    } else {
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", preconditionerName.c_str());
    }
  }


204
  void PetscSolverFeti::updateDofData()
205
206
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
207
208
209

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
210
   
211
212
213
214
215
216
217
    createPrimals();

    createDuals();

    createLagrange();

    createIndexB();
218
219
220
  }


221
  void PetscSolverFeti::createPrimals()
222
  {
223
    FUNCNAME("PetscSolverFeti::createPrimals()");  
224

225
226
227
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

228
229
230
231
232
233
234
    primals.clear();
    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);
235

236
237
238
239

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

240
    globalPrimalIndex.clear();
241
242
243
244
    nRankPrimals = 0;
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it)) {
	globalPrimalIndex[*it] = nRankPrimals;
245
246
247
	nRankPrimals++;
      }

248

249
250
251
    // === Get overall number of primals and rank's displacement in the ===
    // === numbering of the primals.                                    ===

252
    nOverallPrimals = 0;
253
    rStartPrimals = 0;
254
255
256
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankPrimals, rStartPrimals, nOverallPrimals);

257
258
259

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

260
261
262
263
    for (DofMapping::iterator it = globalPrimalIndex.begin();
	 it != globalPrimalIndex.end(); ++it)
      it->second += rStartPrimals;

264
265
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	nRankPrimals, nOverallPrimals);
266

267
268
269
270
271

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

272
273
274
275
276
277
278
279
280
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
    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());


    TEST_EXIT_DBG(nOverallPrimals > 0)
      ("There are zero primal nodes in domain!\n");
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
    
324
325
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341

    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
	if (primals.count(**dofIt) == 0) {
	  boundaryDofRanks[**dofIt].insert(mpiRank);
	  boundaryDofRanks[**dofIt].insert(it->first);
	}
      }
    }

342
343
344
345

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

346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
    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)
	if (primals.count(**dofIt) == 0)
	  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)
	if (primals.count(**dofIt) == 0) {
	  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) == 0)
	  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();
    nRankB = nRankAllDofs - primals.size();
    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) {
      if (primals.count(**it) == 0) {
	duals.insert(**it);
	globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals;
	nRankDuals++;
      }
    }

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

410
411
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
	nRankDuals, nOverallDuals);
412
413
414
415
416
417
418
  }

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

419
420
421
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

422
423
424
425
426
427
428
429
430
    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;
      }
    }

431
432
433
434
435

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

436
    nOverallLagrange = 0;
437
    rStartLagrange = 0;
438
439
440
441
442
443
444
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankLagrange, rStartLagrange, nOverallLagrange);

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

445
446
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	nRankLagrange, nOverallLagrange);
447
448


449
    // === Communicate dofFirstLagrange to all other ranks. ===
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486

    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 (primals.count(**dofIt) == 0) {
	  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)
	if (primals.count(**dofIt) == 0) {
	  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++)
	if (primals.count(*(it->second[i])) == 0)
	  dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++];
487
    }     
488
489
490
491
492
493
494
495
496
  }


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

    globalIndexB.clear();
    DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin();
497
498
499
500

    // === 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.                                 ===
501
502
503
504
505
506
    
    for (int i = 0; i < admin->getUsedSize(); i++)
      if (admin->isDofFree(i) == false && primals.count(i) == 0)
	if (duals.count(i) == 0 && primals.count(i) == 0)
	  globalIndexB[i] = -1;

507
508
509

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

510
    nLocalInterior = 0;
511
512
    for (DofMapping::iterator it = globalIndexB.begin(); 
	 it != globalIndexB.end(); ++it) {
513
514
      it->second = nLocalInterior + rStartB;
      nLocalInterior++;
515
    }
516
    nLocalDuals = duals.size();
517

518
    TEST_EXIT_DBG(nLocalInterior + primals.size() + duals.size() == 
519
520
521
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

522
523
524

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

525
526
527
    for (DofIndexSet::iterator it = duals.begin();
	 it != duals.end(); ++it)
      globalIndexB[*it] = globalDualIndex[*it];
528
529
530
  }


531
  void PetscSolverFeti::createMatLagrange()
532
533
534
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

535
536
    // === Create distributed matrix for Lagrange constraints. ===

537
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
538
539
540
541
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
542
543
544
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

545
546
547
548
549
550
551
    // === 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.                                                     ===

552
553
554
555
    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");

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


589
  void PetscSolverFeti::createSchurPrimalKsp()
590
591
592
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

593
594
595
596
    schurPrimalData.mat_primal_primal = &mat_primal_primal;
    schurPrimalData.mat_primal_b = &mat_primal_b;
    schurPrimalData.mat_b_primal = &mat_b_primal;
    schurPrimalData.ksp_b = &ksp_b;
597

598
599
600
601
602
603
604
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
		 &(schurPrimalData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankPrimals * nComponents, nOverallPrimals * nComponents,
		 &(schurPrimalData.tmp_vec_primal));

605
606
607
608

    MatCreateShell(PETSC_COMM_WORLD,
		   nRankPrimals * nComponents, nRankPrimals * nComponents,
		   nOverallPrimals * nComponents, nOverallPrimals * nComponents,
609
		   &schurPrimalData, 
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
		   &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);
  }


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

625
626
627
628
    schurPrimalData.mat_primal_primal = PETSC_NULL;
    schurPrimalData.mat_primal_b = PETSC_NULL;
    schurPrimalData.mat_b_primal = PETSC_NULL;
    schurPrimalData.ksp_b = PETSC_NULL;
629

630
631
    VecDestroy(&schurPrimalData.tmp_vec_b);
    VecDestroy(&schurPrimalData.tmp_vec_primal);
632

633
634
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
635
636
637
638
639
640
641
  }


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

642
643
    // === Create FETI-DP solver object. ===

644
645
646
647
648
649
    fetiData.mat_primal_primal = &mat_primal_primal;
    fetiData.mat_primal_b = &mat_primal_b;
    fetiData.mat_b_primal = &mat_b_primal;
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.ksp_b = &ksp_b;
    fetiData.ksp_schur_primal = &ksp_schur_primal;
650

651
652
653
654
655
656
657
658
659
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
		 &(fetiData.tmp_vec_b0));
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
		 &(fetiData.tmp_vec_b1));
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankPrimals * nComponents, nOverallPrimals * nComponents,
		 &(fetiData.tmp_vec_primal));
660
661

    MatCreateShell(PETSC_COMM_WORLD,
662
663
		   nRankLagrange * nComponents, nRankLagrange * nComponents,
		   nOverallLagrange * nComponents, nOverallLagrange * nComponents,
664
		   &fetiData, &mat_feti);
665
666
667
668
669
670
671
    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);
672
673


674
    // === Create FETI-DP preconditioner object. ===
675

676
677
678
679
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
680

681
682
683
684
685
686
687
688
689
690
691
692
693
694
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, SAME_NONZERO_PATTERN);
      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;
      
695
696
697
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiDirichletPreconData.tmp_vec_b));      
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
      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;

713
714
715
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiLumpedPreconData.tmp_vec_b));
716
717
718
719
720
721
722
723
724
725
      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);
      
      break;
    }
726
727
728
729
730
731
732
  }
  

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

733
734
    // === Destroy FETI-DP solver object. ===

735
736
737
738
739
740
    fetiData.mat_primal_primal = PETSC_NULL;
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
    fetiData.ksp_b = PETSC_NULL;
    fetiData.ksp_schur_primal = PETSC_NULL;
741

742
743
    VecDestroy(&fetiData.tmp_vec_b0);
    VecDestroy(&fetiData.tmp_vec_b1);
744
    VecDestroy(&fetiData.tmp_vec_primal);
745
746
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
747
748


749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
    // === 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;
    }
778
779
780
781
782
783
784
785
786
  }


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

787
    // === Get local part of the solution for B variables. ===
788
789
790
791
792

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


793
794
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
    
    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);

834
835
836
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
837
838
839
840
841

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


842
    // === And copy from PETSc local vectors to the DOF vectors. ===
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864

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

      for (DofMapping::iterator it = globalIndexB.begin();
	   it != globalIndexB.end(); ++it) {
	int petscIndex = (it->second - rStartB) * nComponents + i;
	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);
865
    VecDestroy(&local_sol_primal);
866
867
868
  }


869
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
870
871
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
872

873
    nComponents = mat->getSize();
874
875
876

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

877
878
    updateDofData();

879
880
881
882
883
884
885

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

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
886
    int nRowsInterior = nLocalInterior * nComponents;
887
    int nRowsDual = nLocalDuals * nComponents;
888
889

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
890
		    nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB,
891
		    30, PETSC_NULL, 0, PETSC_NULL, &mat_b_b);
892
893

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
894
895
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
896
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
897
898

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
899
900
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
901
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
902
903

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
904
905
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
906
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_b);
907

908

909
910
911
912
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
913
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
914
915
916
917
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
918
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
919
920
921
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
922
		      nRowsInterior, nRowsDual, 30, PETSC_NULL,
923
924
925
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
926
		      nRowsDual, nRowsInterior, 30, PETSC_NULL,
927
928
		      &mat_duals_interior);
    }
929

930
931
    
    // === Prepare traverse of sequentially created matrices. ===
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946

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

947
948
949
950
951
952
953
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

954
955
956
957
958
959
960
961
962
963
964
965
966
967
968

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

970
	  bool rowPrimal = primals.count(*cursor) != 0;
971
  
972
973
	  cols.clear();
	  colsOther.clear();
974
	  values.clear();	  
975
	  valuesOther.clear();
976
977
978
979
980
981

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

982
983
984
985
986
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

987
988
989
	    bool colPrimal = primals.count(col(*icursor)) != 0;

	    if (colPrimal) {
990
991
992
993
994
995
996
997
998
999
	      // 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));
1000
	      } else {
1001
1002
1003
1004
1005
1006
1007
1008
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      // Column is not a primal variable.

	      TEST_EXIT_DBG(globalIndexB.count(col(*icursor)))
		("No global B index for DOF %d!\n", col(*icursor));
1009
	      
1010
1011
1012
1013
1014
1015
1016
1017
	      int colIndex = globalIndexB[col(*icursor)] * nComponents + j;

	      if (rowPrimal) {
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(colIndex);
		values.push_back(value(*icursor));
1018
1019
	      }
	    }
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060



	    // === For preconditioner ===

	    if (!rowPrimal && !colPrimal) {
	      int rowIndex = globalIndexB[*cursor] - rStartB;
	      int colIndex = globalIndexB[col(*icursor)] - rStartB;
		
	      if (rowIndex < nLocalInterior) {
		if (colIndex < nLocalInterior) {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB) * nComponents + j;

		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		} else {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		}
	      } else {
		if (colIndex < nLocalInterior) {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB) * nComponents + j;

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		} else {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;

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


1061
	  }
1062

1063
1064
1065
	  if (rowPrimal) {
	    TEST_EXIT_DBG(globalPrimalIndex.count(*cursor))
	      ("Should not happen!\n");
1066

1067
1068
1069
	    int rowIndex = globalPrimalIndex[*cursor] * nComponents + i;
	    MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1070

1071
1072
1073
1074
1075
1076
	    if (colsOther.size())
	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  } else {
	    TEST_EXIT_DBG(globalIndexB.count(*cursor))
	      ("Should not happen!\n");
1077

1078
1079
1080
	    int rowIndex = globalIndexB[*cursor] * nComponents + i;
	    MatSetValues(mat_b_b, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1081

1082
1083
1084
1085
	    if (colsOther.size())
	      MatSetValues(mat_b_primal, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  }
1086
1087
1088
1089

	  // === For preconditioner ===

	  if (!rowPrimal) {
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
	    switch (fetiPreconditioner) {
	    case FETI_DIRICHLET:
	      {
		int rowIndex = globalIndexB[*cursor] - rStartB;
		
		if (rowIndex < nLocalInterior) {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB) * nComponents + i;
		  
		  MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
		  
		  if (colsLocalOther.size()) 
		    MatSetValues(mat_interior_duals, 1, &rowIndex2, colsLocalOther.size(),
				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
		} else {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
		  
		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
		  
		  if (colsLocalOther.size()) 
		    MatSetValues(mat_duals_interior, 1, &rowIndex2, colsLocalOther.size(),
				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
		  
		}
	      }
	      break;
1119

1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
	    case FETI_LUMPED:
	      {
		int rowIndex = globalIndexB[*cursor] - rStartB;
		
		if (rowIndex >= nLocalInterior) {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
		  
		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);		
		}
	      }		
	      break;
	    }	  
1134
	  }
1135
1136
1137
	} 
      }
    }
1138

1139
    // === Start global assembly procedure. ===
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152

    MatAssemblyBegin(mat_b_b, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_b_b, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_primal_primal, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_primal_primal, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_b_primal, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_b_primal, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_primal_b, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_primal_b, MAT_FINAL_ASSEMBLY);

1153

1154
    // === Start global assembly procedure for preconditioner matrices. ===
1155

1156
1157
1158
1159
    if (fetiPreconditioner != FETI_NONE) {
      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
    }
1160

1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
      
      MatAssemblyBegin(mat_interior_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_interior_duals, MAT_FINAL_ASSEMBLY);
      
      MatAssemblyBegin(mat_duals_interior, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_interior, MAT_FINAL_ASSEMBLY);
    }
1171
1172


1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
    // === Create and fill PETSc matrix for Lagrange constraints. ===

    createMatLagrange();

    
    // === Create PETSc solver for the Schur complement on primal variables. ===
    
    createSchurPrimalKsp();


    // === Create PETSc solver for the FETI-DP operator. ===

    createFetiKsp();

    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(PETSC_COMM_WORLD, &ksp_b);
    KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_b, "solver_b_");
    KSPSetFromOptions(ksp_b);
  }

1195

1196
1197
1198
1199
1200
  void PetscSolverFeti::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverFeti::fillPetscRhs()");

    int nComponents = vec->getSize();
1201

1202
1203
1204
1205
1206
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents, &f_b);
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankPrimals * nComponents, 
		 nOverallPrimals * nComponents, &f_primal);
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
    
    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
      for (dofIt.reset(); !dofIt.end(); ++dofIt) {
	int index = dofIt.getDOFIndex();
	if (primals.count(index)) {
	  TEST_EXIT_DBG(globalPrimalIndex.count(index))
	    ("Should not happen!\n");

	  index = globalPrimalIndex[index] * nComponents + i;
	  double value = *dofIt;
1218
	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
1219
1220
1221
1222
1223
1224
	} else {
	  TEST_EXIT_DBG(globalIndexB.count(index))
	    ("Should not happen!\n");

	  index = globalIndexB[index] * nComponents + i;
	  double value = *dofIt;
1225
	  VecSetValues(f_b, 1, &index, &value, ADD_VALUES);
1226
1227
1228
1229
	}      
      }
    }

1230
1231
    VecAssemblyBegin(f_b);
    VecAssemblyEnd(f_b);
1232

1233
1234
    VecAssemblyBegin(f_primal);
    VecAssemblyEnd(f_primal);
1235
1236
1237
  }


1238
  void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
1239
  {
1240
    FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
1241

1242
    Mat  mat_lagrange_transpose;
1243
    MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);
1244

1245
    // === Create nested matrix which will contain the overall FETI system. ===
1246

1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
    Mat A;
    Mat nestedA[3][3];
    nestedA[0][0] = mat_b_b;
    nestedA[0][1] = mat_b_primal;
    nestedA[0][2] = mat_lagrange_transpose;
    nestedA[1][0] = mat_primal_b;
    nestedA[1][1] = mat_primal_primal;
    nestedA[1][2] = PETSC_NULL;
    nestedA[2][0] = mat_lagrange;
    nestedA[2][1] = PETSC_NULL;
    nestedA[2][2] = PETSC_NULL;
1258

1259
    MatCreateNest(PETSC_COMM_WORLD, 3, PETSC_NULL, 3, PETSC_NULL, &(nestedA[0][0]), &A);
1260

1261
1262
1263
    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  
1264
1265


1266
1267
1268
    int nRankNest = (nRankB + nRankPrimals + nRankLagrange) * nComponents;
    int nOverallNest = (nOverallB + nOverallPrimals + nOverallLagrange) * nComponents;
    int rStartNest = (rStartB + rStartPrimals + rStartLagrange) * nComponents;
1269

1270
1271
    {
      // === Test some matrix sizes. ===
1272

1273
1274
1275
1276
1277
      int matRow, matCol;
      MatGetLocalSize(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == nRankNest)("Should not happen!\n");
      mpi::globalAdd(matRow);
      TEST_EXIT_DBG(matRow == nOverallNest)("Should not happen!\n");
1278

1279
1280
1281
      MatGetOwnershipRange(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == rStartNest)("Should not happen!\n");
    }
1282

1283
1284
    // === Create rhs and solution vectors for the overall FETI system. ===

1285
1286
    Vec f, b;
    VecCreateMPI(PETSC_COMM_WORLD, nRankNest, nOverallNest, &f);
1287
1288
1289
1290
    VecDuplicate(f, &b);

    
    // === Fill rhs vector by coping the primal and non primal PETSc vectors. ===
1291

1292
1293
    PetscScalar *local_f_b;
    VecGetArray(f_b, &local_f_b);
1294

1295
1296
    PetscScalar *local_f_primal;
    VecGetArray(f_primal, &local_f_primal);
1297

1298
1299
1300
1301
1302
1303
1304
    {
      int size;
      VecGetLocalSize(f_b, &size);
      TEST_EXIT_DBG(size == nRankB * nComponents)("Should not happen!\n");
      VecGetLocalSize(f_primal, &size);
      TEST_EXIT_DBG(size == nRankPrimals * nComponents)("Should not happen!\n");
    }
1305

1306
1307
    PetscScalar *local_f;
    VecGetArray(f, &local_f);