PetscSolverFeti.cc 45.2 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
    VecDuplicate(f_b, &(schurPrimalData.tmp_vec_b));
    VecDuplicate(f_primal, &(schurPrimalData.tmp_vec_primal));
600
601
602
603

    MatCreateShell(PETSC_COMM_WORLD,
		   nRankPrimals * nComponents, nRankPrimals * nComponents,
		   nOverallPrimals * nComponents, nOverallPrimals * nComponents,
604
		   &schurPrimalData, 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
		   &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()");

620
621
622
623
    schurPrimalData.mat_primal_primal = PETSC_NULL;
    schurPrimalData.mat_primal_b = PETSC_NULL;
    schurPrimalData.mat_b_primal = PETSC_NULL;
    schurPrimalData.ksp_b = PETSC_NULL;
624

625
626
    VecDestroy(&schurPrimalData.tmp_vec_b);
    VecDestroy(&schurPrimalData.tmp_vec_primal);
627

628
629
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
630
631
632
633
634
635
636
  }


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

637
638
    // === Create FETI-DP solver object. ===

639
640
641
642
643
644
    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;
645

646
647
    VecDuplicate(f_b, &(fetiData.tmp_vec_b0));
    VecDuplicate(f_b, &(fetiData.tmp_vec_b1));
648
    VecDuplicate(f_primal, &(fetiData.tmp_vec_primal));
649
650

    MatCreateShell(PETSC_COMM_WORLD,
651
652
		   nRankLagrange * nComponents, nRankLagrange * nComponents,
		   nOverallLagrange * nComponents, nOverallLagrange * nComponents,
653
		   &fetiData, &mat_feti);
654
655
656
657
658
659
660
    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);
661
662


663
    // === Create FETI-DP preconditioner object. ===
664

665
666
667
668
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
669

670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
    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;
      
      VecDuplicate(f_b, &(fetiDirichletPreconData.tmp_vec_b));      
      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;

      VecDuplicate(f_b, &(fetiLumpedPreconData.tmp_vec_b));
      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;
    }
711
712
713
714
715
716
717
  }
  

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

718
719
    // === Destroy FETI-DP solver object. ===

720
721
722
723
724
725
    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;
726

727
728
    VecDestroy(&fetiData.tmp_vec_b0);
    VecDestroy(&fetiData.tmp_vec_b1);
729
    VecDestroy(&fetiData.tmp_vec_primal);
730
731
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
732
733


734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
    // === 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;
    }
763
764
765
766
767
768
769
770
771
  }


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

772
    // === Get local part of the solution for B variables. ===
773
774
775
776
777

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


778
779
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
    
    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);

819
820
821
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
822
823
824
825
826

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


827
    // === And copy from PETSc local vectors to the DOF vectors. ===
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849

    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);
850
    VecDestroy(&local_sol_primal);
851
852
853
  }


854
855
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat, 
					SystemVector *vec)
856
857
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
858

859
860
861
862
    nComponents = vec->getSize();

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

863
864
    updateDofData();

865
866
867
868
869
870
871

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

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
872
    int nRowsInterior = nLocalInterior * nComponents;
873
    int nRowsDual = nLocalDuals * nComponents;
874
875

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
876
877
		    nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_b_b);
878
879

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
880
881
882
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
		    10, PETSC_NULL, 10, PETSC_NULL, &mat_primal_primal);
883
884

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
885
886
887
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_b_primal);
888
889

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
890
891
892
893
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_primal_b);

894

895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsDual, nRowsDual, 100, PETSC_NULL,
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsInterior, nRowsInterior, 100, PETSC_NULL,
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsInterior, nRowsDual, 100, PETSC_NULL,
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsDual, nRowsInterior, 100, PETSC_NULL,
		      &mat_duals_interior);
    }
915

916
917
    
    // === Prepare traverse of sequentially created matrices. ===
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932

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

933
934
935
936
937
938
939
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

940
941
942
943
944
945
946
947
948
949
950
951
952
953
954

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

956
	  bool rowPrimal = primals.count(*cursor) != 0;
957
  
958
959
	  cols.clear();
	  colsOther.clear();
960
	  values.clear();	  
961
	  valuesOther.clear();
962
963
964
965
966
967

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

968
969
970
971
972
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

973
974
975
	    bool colPrimal = primals.count(col(*icursor)) != 0;

	    if (colPrimal) {
976
977
978
979
980
981
982
983
984
985
	      // 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));
986
	      } else {
987
988
989
990
991
992
993
994
		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));
995
	      
996
997
998
999
1000
1001
1002
1003
	      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));
1004
1005
	      }
	    }
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
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



	    // === 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));
		}
	      }		
	    }


1047
	  }
1048

1049
1050
1051
	  if (rowPrimal) {
	    TEST_EXIT_DBG(globalPrimalIndex.count(*cursor))
	      ("Should not happen!\n");
1052

1053
1054
1055
	    int rowIndex = globalPrimalIndex[*cursor] * nComponents + i;
	    MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1056

1057
1058
1059
1060
1061
1062
	    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");
1063

1064
1065
1066
	    int rowIndex = globalIndexB[*cursor] * nComponents + i;
	    MatSetValues(mat_b_b, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1067

1068
1069
1070
1071
	    if (colsOther.size())
	      MatSetValues(mat_b_primal, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  }
1072
1073
1074
1075

	  // === For preconditioner ===

	  if (!rowPrimal) {
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
	    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;
1105

1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
	    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;
	    }	  
1120
	  }
1121
1122
1123
1124
	} 
      }
    }
    
1125

1126
    // === Start global assembly procedure. ===
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139

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

1140

1141
    // === Start global assembly procedure for preconditioner matrices. ===
1142

1143
1144
1145
1146
    if (fetiPreconditioner != FETI_NONE) {
      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
    }
1147

1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
    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);
    }
1158
1159


1160
    // === Create and fill PETSc's right hand side vectors. ===
1161

1162
1163
1164
    VecCreate(PETSC_COMM_WORLD, &f_b);
    VecSetSizes(f_b, nRankB * nComponents, nOverallB * nComponents);
    VecSetType(f_b, VECMPI);
1165

1166
1167
    VecCreate(PETSC_COMM_WORLD, &f_primal);
    VecSetSizes(f_primal, nRankPrimals * nComponents, 
1168
		nOverallPrimals * nComponents);
1169
    VecSetType(f_primal, VECMPI);
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
    
    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;
1181
	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
1182
1183
1184
1185
1186
1187
	} else {
	  TEST_EXIT_DBG(globalIndexB.count(index))
	    ("Should not happen!\n");

	  index = globalIndexB[index] * nComponents + i;
	  double value = *dofIt;
1188
	  VecSetValues(f_b, 1, &index, &value, ADD_VALUES);
1189
1190
1191
1192
	}      
      }
    }

1193
1194
    VecAssemblyBegin(f_b);
    VecAssemblyEnd(f_b);
1195

1196
1197
    VecAssemblyBegin(f_primal);
    VecAssemblyEnd(f_primal);
1198
1199


1200
    // === Create and fill PETSc matrix for Lagrange constraints. ===
1201

1202
    createMatLagrange();
1203
1204

    
1205
1206
1207
1208
1209
1210
    // === Create PETSc solver for the Schur complement on primal variables. ===
    
    createSchurPrimalKsp();


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

    createFetiKsp();
1213
1214
1215
  }


1216
  void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
1217
  {
1218
    FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
1219

1220
1221
1222
    // Create transpose of Lagrange matrix.
    Mat mat_lagrange_transpose;
    MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);
1223
1224


1225
    // === Create nested matrix which will contain the overall FETI system. ===
1226

1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
    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;
1238

1239
    MatCreateNest(PETSC_COMM_WORLD, 3, PETSC_NULL, 3, PETSC_NULL, &(nestedA[0][0]), &A);
1240

1241
1242
1243
    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  
1244
1245


1246
1247
1248
    int nRankNest = (nRankB + nRankPrimals + nRankLagrange) * nComponents;
    int nOverallNest = (nOverallB + nOverallPrimals + nOverallLagrange) * nComponents;
    int rStartNest = (rStartB + rStartPrimals + rStartLagrange) * nComponents;
1249

1250
1251
    {
      // === Test some matrix sizes. ===
1252

1253
1254
1255
1256
1257
      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");
1258

1259
1260
1261
      MatGetOwnershipRange(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == rStartNest)("Should not happen!\n");
    }
1262

1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
    // === Create rhs and solution vectors for the overall FETI system. ===

    Vec f;
    VecCreate(PETSC_COMM_WORLD, &f);
    VecSetSizes(f, nRankNest, nOverallNest);
    VecSetType(f, VECMPI);

    Vec b;
    VecDuplicate(f, &b);

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

1276
1277
    PetscScalar *local_f_b;
    VecGetArray(f_b, &local_f_b);
1278

1279
1280
    PetscScalar *local_f_primal;
    VecGetArray(f_primal, &local_f_primal);
1281

1282
1283
1284
1285
1286
1287
1288
    {
      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");
    }
1289

1290
1291
    PetscScalar *local_f;
    VecGetArray(f, &local_f);
1292

1293
1294
1295
1296
1297
    int index = 0;
    for (int i = 0; i < nRankB * nComponents; i++)
      local_f[index++] = local_f_b[i];
    for (int i = 0; i < nRankPrimals * nComponents; i++)
      local_f[index++] = local_f_primal[i];
1298

1299
1300
1301
    VecRestoreArray(f, &local_f);  
    VecRestoreArray(f_b, &local_f_b);
    VecRestoreArray(f_primal, &local_f_primal);
1302

1303
1304
    
    // === Create solver and solve the overall FETI system. ===
1305

1306
1307
1308
1309
    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, A, A, SAME_NONZERO_PATTERN);
    KSPSetFromOptions(ksp);
1310
1311


1312
    KSPSolve(ksp, f, b);