PetscSolverFeti.cc 52.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


#include "parallel/PetscSolverFeti.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
16
#include "io/VtkWriter.h"
17
18
19
20
21

namespace AMDiS {

  using namespace std;

22
#ifdef HAVE_PETSC_DEV 
23

24
25
26
27
28
29
30
  // 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);
31
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
32
33

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
34
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
35
36
37
38
39
40
41
42
43
44
45
    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

Thomas Witkowski's avatar
Thomas Witkowski committed
53
54
55
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
56

Thomas Witkowski's avatar
Thomas Witkowski committed
57
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
58
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
59
60
61
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
62

Thomas Witkowski's avatar
Thomas Witkowski committed
63
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
64
65
66
67
68

    return 0;
  }


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

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


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

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

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

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


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

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

103
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
104
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
    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);
139
140


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

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

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

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

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


    // === K_DD ===

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

163

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

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

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

175
176

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

    return 0;
  }


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

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

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


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

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
219
   
220
221
222
223
224
225
226
    createPrimals();

    createDuals();

    createLagrange();

    createIndexB();
227
228
229
  }


230
  void PetscSolverFeti::createPrimals()
231
  {
232
    FUNCNAME("PetscSolverFeti::createPrimals()");  
233

234
235
236
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

237
238
239
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
240
241
242
243
244
245
    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);
246

247
248
249
250

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

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

259

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

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

268
269
270

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

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

275
276
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	nRankPrimals, nOverallPrimals);
277

278
279
280
281
282

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

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

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

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

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int i = 0;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false)
	  globalPrimalIndex[**dofIt] = stdMpi.getRecvData(it->first)[i++];
      }
    }

    TEST_EXIT_DBG(primals.size() == globalPrimalIndex.size())
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
       primals.size(), globalPrimalIndex.size());
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
    
331
332
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
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
342
	if (globalPrimalIndex.count(**dofIt) == 0) {
343
344
345
346
347
348
	  boundaryDofRanks[**dofIt].insert(mpiRank);
	  boundaryDofRanks[**dofIt].insert(it->first);
	}
      }
    }

349
350
351
352

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

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

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

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

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

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

429
430
431
432
433
434
435
436
437
    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;
      }
    }

438
439
440
441
442

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

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

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

452
453
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	nRankLagrange, nOverallLagrange);
454
455


456
    // === Communicate dofFirstLagrange to all other ranks. ===
457
458
459
460
461
462
463

    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) {
464
	if (globalPrimalIndex.count(**dofIt) == 0) {
465
466
467
468
469
470
471
472
473
474
475
476
	  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)
477
	if (globalPrimalIndex.count(**dofIt) == 0) {
478
479
480
481
482
483
484
485
486
487
488
489
490
491
	  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++)
492
	if (globalPrimalIndex.count(*(it->second[i])) == 0)
493
	  dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++];
494
    }     
495
496
497
498
499
500
501
  }


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

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

    // === 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.                                 ===
508
509
    
    for (int i = 0; i < admin->getUsedSize(); i++)
510
511
      if (admin->isDofFree(i) == false && globalPrimalIndex.count(i) == 0)
	if (duals.count(i) == 0 && globalPrimalIndex.count(i) == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
512
	  localIndexB[i] = -1;
513

514
515
516

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

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

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

529
530
531

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

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


538
  void PetscSolverFeti::createMatLagrange()
539
540
541
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

542
543
    // === Create distributed matrix for Lagrange constraints. ===

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

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

559
560
561
562
    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");

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


596
  void PetscSolverFeti::createSchurPrimalKsp()
597
598
599
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
      schurPrimalData.mat_primal_primal = &mat_primal_primal;
      schurPrimalData.mat_primal_b = &mat_primal_b;
      schurPrimalData.mat_b_primal = &mat_b_primal;
      schurPrimalData.fetiSolver = this;
      
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankPrimals * nComponents, nOverallPrimals * nComponents,
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
		     nRankPrimals * nComponents, nRankPrimals * nComponents,
		     nOverallPrimals * nComponents, nOverallPrimals * nComponents,
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
    } else {
628
629
      MSG("Create direct schur primal solver!\n");

630
631
632
633
634
635
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

      int nRowsRankPrimal = nRankPrimals * nComponents;
      int nRowsOverallPrimal = nOverallPrimals * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
636
637
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
638

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

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

650
651
652
653
654
655
656
657
658
659
660
      int nLocalPrimals = globalPrimalIndex.size() * nComponents;
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

      for (int i = 0; i < nRankB * nComponents; i++) {
	PetscInt row = rStartB * nComponents + i;
	MatGetRow(mat_b_primal, row, &nCols, &cols, &values);

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

Thomas Witkowski's avatar
Thomas Witkowski committed
661
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
		globalPrimalIndex.size() * nComponents)("Should not happen!\n");
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
	VecCreateSeq(PETSC_COMM_SELF, nRankB * nComponents, &tmpVec);

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);

       	KSPSolve(ksp_b, tmpVec, tmpVec);

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

	VecDestroy(&tmpVec);
      }

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

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

      MatInfo minfo;
      MatGetInfo(mat_primal_primal, MAT_GLOBAL_SUM, &minfo);
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_primal_primal, 
		      mat_primal_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
713

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


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

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

730
731
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
732

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


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

745
746
    // === Create FETI-DP solver object. ===

747
748
749
    fetiData.mat_primal_b = &mat_primal_b;
    fetiData.mat_b_primal = &mat_b_primal;
    fetiData.mat_lagrange = &mat_lagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
750
    fetiData.fetiSolver = this;
751
    fetiData.ksp_schur_primal = &ksp_schur_primal;
752

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

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


776
    // === Create FETI-DP preconditioner object. ===
777

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

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

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

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

838
839
    // === Destroy FETI-DP solver object. ===

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

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


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


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

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

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


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

940
941
942
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
943
944
945
946

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

947
    // === And copy from PETSc local vectors to the DOF vectors. ===
948
949
950
951

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

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


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

976
    nComponents = mat->getSize();
977
978
979

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

980
981
    updateDofData();

982
983
984
985
986
987
988

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

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

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

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
996
997
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
998
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
999
1000

    MatCreateMPIAIJ(PETSC_COMM_WORLD,