PetscSolverFeti.cc 33.5 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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
  // 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);
    PetscSchurPrimalData* data = static_cast<PetscSchurPrimalData*>(ctx);

    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)
  {
    // F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)

    void *ctx;
    MatShellGetContext(mat, &ctx);
    PetscFetiData* data = static_cast<PetscFetiData*>(ctx);

    // y = L inv(K_BB) trans(L) x
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);

    // tmp_vec_primal = inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);

    // tmp_vec_lagrange = L inv(K_BB) K_BPi tmp_vec_primal
    //                  = L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);

    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);

    return 0;
  }


73
  void PetscSolverFeti::updateDofData()
74
75
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
76
77
78
79
80
81

    TEST_EXIT(meshDistributor->getMesh()->getDim() == 2)
      ("Works for 2D problems only!");

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
82
   
83
84
85
86
87
88
89
    createPrimals();

    createDuals();

    createLagrange();

    createIndexB();
90
91
92
  }


93
  void PetscSolverFeti::createPrimals()
94
  {
95
    FUNCNAME("PetscSolverFeti::createPrimals()");  
96

97
98
99
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

100
101
102
103
104
105
106
    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);
107

108
109
110
111

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

112
    globalPrimalIndex.clear();
113
114
115
116
    nRankPrimals = 0;
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it)) {
	globalPrimalIndex[*it] = nRankPrimals;
117
118
119
	nRankPrimals++;
      }

120

121
122
123
    // === Get overall number of primals and rank's displacement in the ===
    // === numbering of the primals.                                    ===

124
    nOverallPrimals = 0;
125
    rStartPrimals = 0;
126
127
128
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankPrimals, rStartPrimals, nOverallPrimals);

129
130
131

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

132
133
134
135
136
137
    for (DofMapping::iterator it = globalPrimalIndex.begin();
	 it != globalPrimalIndex.end(); ++it)
      it->second += rStartPrimals;

    MSG_DBG("nRankPrimals = %d   nOverallPrimals = %d\n",
	    nRankPrimals, nOverallPrimals);
138

139
140
141
142
143

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

144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
    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()");
    
196
197
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213

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

214
215
216
217

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

218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
    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);

    MSG_DBG("nRankDuals = %d   nOverallDuals = %d\n",
	    nRankDuals, nOverallDuals);
  }

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

291
292
293
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

294
295
296
297
298
299
300
301
302
    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;
      }
    }

303
304
305
306
307

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

308
    nOverallLagrange = 0;
309
    rStartLagrange = 0;
310
311
312
313
314
315
316
317
318
319
320
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankLagrange, rStartLagrange, nOverallLagrange);

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

    MSG_DBG("nRankLagrange = %d  nOverallLagrange = %d\n",
	    nRankLagrange, nOverallLagrange);


321
    // === Communicate dofFirstLagrange to all other ranks. ===
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358

    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++];
359
    }     
360
361
362
363
364
365
366
367
368
  }


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

    globalIndexB.clear();
    DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin();
369
370
371
372

    // === 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.                                 ===
373
374
375
376
377
378
    
    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;

379
380
381

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

382
383
384
385
386
387
388
389
390
391
392
    int nInterior = 0;
    for (DofMapping::iterator it = globalIndexB.begin(); 
	 it != globalIndexB.end(); ++it) {
      it->second = nInterior + rStartB;
      nInterior++;
    }

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

393
394
395

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

396
397
398
    for (DofIndexSet::iterator it = duals.begin();
	 it != duals.end(); ++it)
      globalIndexB[*it] = globalDualIndex[*it];
399
400
401
  }


402
  void PetscSolverFeti::createMatLagrange()
403
404
405
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

406
407
    // === Create distributed matrix for Lagrange constraints. ===

408
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
409
410
411
412
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
413
414
415
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

416
417
418
419
420
421
422
    // === 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.                                                     ===

423
424
425
426
    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");

427
      // Global index of the first Lagrange constriant for this node.
428
      int index = dofFirstLagrange[*it];
429
      // Copy set of all ranks that contain this dual node.
430
      vector<int> W(boundaryDofRanks[*it].begin(), boundaryDofRanks[*it].end());
431
      // Number of ranks that contain this dual node.
432
433
434
435
436
437
438
      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++) {
439
440
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
	    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);
  }


460
  void PetscSolverFeti::createSchurPrimalKsp()
461
462
463
464
465
466
467
468
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

    petscSchurPrimalData.mat_primal_primal = &mat_primal_primal;
    petscSchurPrimalData.mat_primal_b = &mat_primal_b;
    petscSchurPrimalData.mat_b_primal = &mat_b_primal;
    petscSchurPrimalData.ksp_b = &ksp_b;

469
470
    VecDuplicate(f_b, &(petscSchurPrimalData.tmp_vec_b));
    VecDuplicate(f_primal, &(petscSchurPrimalData.tmp_vec_primal));
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515

    MatCreateShell(PETSC_COMM_WORLD,
		   nRankPrimals * nComponents, nRankPrimals * nComponents,
		   nOverallPrimals * nComponents, nOverallPrimals * nComponents,
		   &petscSchurPrimalData, 
		   &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()");

    petscSchurPrimalData.mat_primal_primal = PETSC_NULL;
    petscSchurPrimalData.mat_primal_b = PETSC_NULL;
    petscSchurPrimalData.mat_b_primal = PETSC_NULL;
    petscSchurPrimalData.ksp_b = PETSC_NULL;

    VecDestroy(petscSchurPrimalData.tmp_vec_b);
    VecDestroy(petscSchurPrimalData.tmp_vec_primal);

    MatDestroy(mat_schur_primal);
    KSPDestroy(ksp_schur_primal);
  }


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

    petscFetiData.mat_primal_primal = &mat_primal_primal;
    petscFetiData.mat_primal_b = &mat_primal_b;
    petscFetiData.mat_b_primal = &mat_b_primal;
    petscFetiData.mat_lagrange = &mat_lagrange;
    petscFetiData.ksp_b = &ksp_b;
    petscFetiData.ksp_schur_primal = &ksp_schur_primal;


516
517
    VecDuplicate(f_b, &(petscFetiData.tmp_vec_b));
    VecDuplicate(f_primal, &(petscFetiData.tmp_vec_primal));
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
    MatGetVecs(mat_lagrange, PETSC_NULL, &(petscFetiData.tmp_vec_lagrange));


    MatCreateShell(PETSC_COMM_WORLD,
		   nRankLagrange, nRankLagrange,
		   nOverallLagrange, nOverallLagrange,
		   &petscFetiData, &mat_feti);
    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);
  }
  

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

    petscFetiData.mat_primal_primal = PETSC_NULL;
    petscFetiData.mat_primal_b = PETSC_NULL;
    petscFetiData.mat_b_primal = PETSC_NULL;
    petscFetiData.mat_lagrange = PETSC_NULL;
    petscFetiData.ksp_b = PETSC_NULL;
    petscFetiData.ksp_schur_primal = PETSC_NULL;

    VecDestroy(petscFetiData.tmp_vec_b);
    VecDestroy(petscFetiData.tmp_vec_primal);
    VecDestroy(petscFetiData.tmp_vec_lagrange);

    MatDestroy(mat_feti);
    KSPDestroy(ksp_feti);
  }


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

561
    // === Get local part of the solution for B variables. ===
562
563
564
565
566

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


567
568
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
    
    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);

    ISDestroy(globalIs);
    ISDestroy(localIs);    
    VecScatterDestroy(primalScatter);    

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


616
    // === And copy from PETSc local vectors to the DOF vectors. ===
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642

    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);
    VecDestroy(local_sol_primal);
  }


643
644
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat, 
					SystemVector *vec)
645
646
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
647

648
649
650
651
    nComponents = vec->getSize();

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

652
653
    updateDofData();

654
655
656
657
658
659
660

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

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
661
662

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
663
664
		    nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_b_b);
665
666

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
667
668
669
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
		    10, PETSC_NULL, 10, PETSC_NULL, &mat_primal_primal);
670
671

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
672
673
674
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_b_primal);
675
676

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
677
678
679
680
681
682
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_primal_b);

    
    // === Prepare traverse of sequentially created matrices. ===
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697

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

698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735

    // === 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) {
	  bool rowPrimal = primals.count(*cursor) != 0;
	  
	  cols.clear();
	  values.clear();
	  
	  colsOther.clear();
	  valuesOther.clear();
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

	    if (primals.count(col(*icursor)) != 0) {
	      // 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));
736
	      } else {
737
738
739
740
741
742
743
744
		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));
745
	      
746
747
748
749
750
751
752
753
	      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));
754
755
	      }
	    }
756
	  }
757

758
759
760
	  if (rowPrimal) {
	    TEST_EXIT_DBG(globalPrimalIndex.count(*cursor))
	      ("Should not happen!\n");
761

762
763
764
	    int rowIndex = globalPrimalIndex[*cursor] * nComponents + i;
	    MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
765

766
767
768
769
770
771
	    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");
772

773
774
775
	    int rowIndex = globalIndexB[*cursor] * nComponents + i;
	    MatSetValues(mat_b_b, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
776

777
778
779
780
781
782
783
784
	    if (colsOther.size())
	      MatSetValues(mat_b_primal, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  }
	} 
      }
    }
    
785

786
    // === Start global assembly procedure. ===
787
788
789
790
791
792
793
794
795
796
797
798
799
800

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

801
    // === Create and fill PETSc's right hand side vectors. ===
802

803
804
805
    VecCreate(PETSC_COMM_WORLD, &f_b);
    VecSetSizes(f_b, nRankB * nComponents, nOverallB * nComponents);
    VecSetType(f_b, VECMPI);
806

807
808
    VecCreate(PETSC_COMM_WORLD, &f_primal);
    VecSetSizes(f_primal, nRankPrimals * nComponents, 
809
		nOverallPrimals * nComponents);
810
    VecSetType(f_primal, VECMPI);
811
812
813
814
815
816
817
818
819
820
821
    
    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;
822
	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
823
824
825
826
827
828
	} else {
	  TEST_EXIT_DBG(globalIndexB.count(index))
	    ("Should not happen!\n");

	  index = globalIndexB[index] * nComponents + i;
	  double value = *dofIt;
829
	  VecSetValues(f_b, 1, &index, &value, ADD_VALUES);
830
831
832
833
	}      
      }
    }

834
835
    VecAssemblyBegin(f_b);
    VecAssemblyEnd(f_b);
836

837
838
    VecAssemblyBegin(f_primal);
    VecAssemblyEnd(f_primal);
839
840


841
    // === Create and fill PETSc matrix for Lagrange constraints. ===
842

843
    createMatLagrange();
844
845

    
846
847
848
849
850
851
    // === Create PETSc solver for the Schur complement on primal variables. ===
    
    createSchurPrimalKsp();


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

    createFetiKsp();
854
855
856
  }


857
  void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
858
  {
859
    FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
860

861
862
863
    // Create transpose of Lagrange matrix.
    Mat mat_lagrange_transpose;
    MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);
864
865


866
    // === Create nested matrix which will contain the overall FETI system. ===
867

868
869
870
871
872
873
874
875
876
877
878
    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;
879

880
    MatCreateNest(PETSC_COMM_WORLD, 3, PETSC_NULL, 3, PETSC_NULL, &(nestedA[0][0]), &A);
881

882
883
884
    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  
885
886


887
888
889
    int nRankNest = (nRankB + nRankPrimals) * nComponents + nRankLagrange;    
    int nOverallNest = (nOverallB + nOverallPrimals) * nComponents + nOverallLagrange;
    int rStartNest = (rStartB + rStartPrimals) * nComponents + rStartLagrange;
890

891
892
    {
      // === Test some matrix sizes. ===
893

894
895
896
897
898
      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");
899

900
901
902
      MatGetOwnershipRange(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == rStartNest)("Should not happen!\n");
    }
903

904
905
906
907
908
909
910
911
912
913
914
915
    // === 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. ===
916

917
918
    PetscScalar *local_f_b;
    VecGetArray(f_b, &local_f_b);
919

920
921
    PetscScalar *local_f_primal;
    VecGetArray(f_primal, &local_f_primal);
922

923
924
925
926
927
928
929
    {
      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");
    }
930

931
932
    PetscScalar *local_f;
    VecGetArray(f, &local_f);
933

934
935
936
937
938
    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];
939

940
941
942
    VecRestoreArray(f, &local_f);  
    VecRestoreArray(f_b, &local_f_b);
    VecRestoreArray(f_primal, &local_f_primal);
943

944
945
    
    // === Create solver and solve the overall FETI system. ===
946

947
948
949
950
    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, A, A, SAME_NONZERO_PATTERN);
    KSPSetFromOptions(ksp);
951
952


953
    KSPSolve(ksp, f, b);
954
955


956
957
958
959
960
    // === Reconstruct FETI solution vectors. ===
    
    Vec u_b, u_primal;
    VecDuplicate(f_b, &u_b);
    VecDuplicate(f_primal, &u_primal);
961
962
    

963
964
    PetscScalar *local_b;
    VecGetArray(b, &local_b);
965

966
967
    PetscScalar *local_u_b;
    VecGetArray(u_b, &local_u_b);
968

969
970
    PetscScalar *local_u_primal;
    VecGetArray(u_primal, &local_u_primal);
971

972
973
974
975
976
    index = 0;
    for (int i = 0; i < nRankB * nComponents; i++)
      local_u_b[i] = local_b[index++];
    for (int i = 0; i < nRankPrimals * nComponents; i++)
      local_u_primal[i] = local_b[index++];
977

978
979
980
    VecRestoreArray(f, &local_b);
    VecRestoreArray(u_b, &local_u_b);
    VecRestoreArray(u_primal, &local_u_primal);
981

982
    recoverSolution(u_b, u_primal, vec);
983

984
985
986
987
    VecDestroy(u_b);
    VecDestroy(u_primal);
    VecDestroy(b);
    VecDestroy(f);
988

989
990
    KSPDestroy(ksp);
  }
991
992


993
994
995
  void PetscSolverFeti::solveReducedFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
996

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

999
1000
1001
1002
1003
1004
1005
    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);

    // RHS and solution vector.
    Vec vec_rhs;
1006

1007
1008
1009
1010
1011
1012
1013
1014
    // Some temporary vectors.
    Vec tmp_b0, tmp_b1, tmp_lagrange0, tmp_primal0, tmp_primal1;
    MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange0);
    MatGetVecs(mat_lagrange, PETSC_NULL, &vec_rhs);
    MatGetVecs(mat_b_b, PETSC_NULL, &tmp_b0);
    MatGetVecs(mat_b_b, PETSC_NULL, &tmp_b1);
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal0);
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal1);
1015
1016


1017
    // === Create new rhs ===
1018

1019
1020
1021
    // vec_rhs = L * inv(K_BB) * f_b
    KSPSolve(ksp_b, f_b, tmp_b0);
    MatMult(mat_lagrange, tmp_b0, vec_rhs);
1022

1023
1024
    // tmp_primal0 = M_PiB * inv(K_BB) * f_b
    MatMult(mat_primal_b, tmp_b0, tmp_primal0);
1025

1026
1027
    // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_b
    VecAXPBY(tmp_primal0, -1.0, 1.0, f_primal);
1028

1029
1030
    // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_b)
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1031

1032
1033
1034
1035
    //
    MatMult(mat_b_primal, tmp_primal0, tmp_b0);
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
    MatMult(mat_lagrange, tmp_b0, tmp_lagrange0);
1036

1037
1038
    //
    VecAXPBY(vec_rhs, 1.0, 1.0, tmp_lagrange0);
1039
1040


1041
    // === Solve with FETI-DP operator. ===
1042

1043
    KSPSolve(ksp_feti, vec_rhs, vec_rhs);
1044
1045

   
1046
    // === Solve for u_primals. ===
1047

1048
    VecCopy(f_primal, tmp_primal0);
1049

1050
1051
    KSPSolve(ksp_b, f_b, tmp_b0);
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
1052

1053
    VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1);
1054

1055
1056
1057
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b0);
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
1058

1059
1060
    VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1);
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1061
1062

    
1063
    // === Solve for u_b. ===
1064

1065
1066
1067
    VecCopy(f_b, tmp_b0);
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
1068

1069
1070
    MatMult(mat_b_primal, tmp_primal0, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
1071

1072
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
1073
1074


1075
    // === And recover AMDiS solution vectors. ===
1076
    
1077
    recoverSolution(tmp_b0, tmp_primal0, vec);
1078
1079


1080
    // === Destroy all data structures. ===
1081
    
1082
1083
1084
1085
1086
1087
    VecDestroy(vec_rhs);
    VecDestroy(tmp_b0);
    VecDestroy(tmp_b1);
    VecDestroy(tmp_lagrange0);
    VecDestroy(tmp_primal0);
    VecDestroy(tmp_primal1);
1088
1089
	    

1090
    KSPDestroy(ksp_b);
1091

1092
1093
1094
1095
1096
    MatDestroy(mat_b_b);
    MatDestroy(mat_primal_primal);
    MatDestroy(mat_b_primal);
    MatDestroy(mat_primal_b);
    MatDestroy(mat_lagrange);
1097

1098
1099
    VecDestroy(f_b);
    VecDestroy(f_primal);
1100

1101
    destroySchurPrimalKsp();
1102

1103
1104
    destroyFetiKsp();   
  }
1105

1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121

  void PetscSolverFeti::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverFeti::solvePetscMatrix()");

    int debug = 0;
    Parameters::get("parallel->debug feti", debug);

    if (debug) {
      WARNING("FETI matrix is solved globally, thus without reducing to the lagrange multipliers!\n");

      solveFetiMatrix(vec);
    } else {
      solveReducedFetiMatrix(vec);
    }
      
1122
1123
1124
1125
  }
#endif

}