PetscSolverFeti.cc 52.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
//
// 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"
14
#include "parallel/PetscSolverFetiStructs.h"
15
16
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
17
#include "io/VtkWriter.h"
18
19
20
21
22

namespace AMDiS {

  using namespace std;

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
54
    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
55

Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
57
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
60
    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);
61

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

    return 0;
  }


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

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


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

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

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

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


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

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

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


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

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

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

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

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


    // === K_DD ===

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

162

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

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

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

174
175

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

    return 0;
  }


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

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

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


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

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

    createDuals();

    createLagrange();
224
    
225
    createIndexB();
226
227
228
  }


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

233
234
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

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

248
249
250
251

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

252
    primalDofMap.addFeSpace(feSpace);
253
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
254
255
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
256

257
    primalDofMap[feSpace].update();
258

259
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
260
261
	primalDofMap[feSpace].nRankDofs, 
	primalDofMap[feSpace].nOverallDofs);
262

263
264
265
266

    // === Communicate primal's global index from ranks that own the     ===
    // === primals to ranks that contain this primals but are not owning ===
    // === them.                                                         ===
267
268
    
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
269

270
    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
271
272
273
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
274
275
276
277
	if (primalDofMap[feSpace].isSet(it.getDofIndex())) {
	  DegreeOfFreedom d = primalDofMap[feSpace][it.getDofIndex()];
	  stdMpi.getSendData(it.getRank()).push_back(d);
	}
278

279
280
    stdMpi.updateSendDataSize();

281
282
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
283
      bool recvFromRank = false;
284
285
286
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
287
288
289
	  recvFromRank = true;
	  break;
	}
290
      }
291
292

      if (recvFromRank) 
293
	stdMpi.recv(it.getRank());
294
    }
295

296
297
    stdMpi.startCommunication();

298
299
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
300
      int i = 0;
301
302
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
303
304
305
306
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[i++];
	  primalDofMap[feSpace].insert(it.getDofIndex(), d);
	}
307
308
309
      }
    }

310
    TEST_EXIT_DBG(primals.size() == primalDofMap[feSpace].size())
311
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
312
       primals.size(), primalDofMap[feSpace].size());
313
314
315
316
317
318
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
319
320

    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
321
    
322
323
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
324
325
326

    boundaryDofRanks.clear();

327
328
329
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
330
	// If DOF is not primal, i.e., its a dual node
331
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
332
333
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
334
335
	}
      }
336
	
337
338
339
340

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

341
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
342
343
344
345

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
346
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
347
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
348
349
350

    stdMpi.updateSendDataSize();

351
352
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
353
      bool recvFromRank = false;
354
      for (; !it.endDofIter(); it.nextDof()) {
355
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
356
357
358
	  recvFromRank = true;
	  break;
	}
359
      }
360
361

      if (recvFromRank)
362
	stdMpi.recv(it.getRank());
363
    }
364

365
366
    stdMpi.startCommunication();

367
368
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
369
      int i = 0;
370
      for (; !it.endDofIter(); it.nextDof())
371
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
372
373
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
374
375
376
377
    }

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

378
    dualDofMap.addFeSpace(feSpace);
379

380
    int nRankAllDofs = feSpace->getAdmin()->getUsedDofs();
381
    nRankB = nRankAllDofs - primalDofMap[feSpace].size();
382
383
384
385
386
    nOverallB = 0;
    rStartB = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankB, rStartB, nOverallB);
    DofContainer allBoundaryDofs;
387
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
388
389
390
    int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size();

    for (DofContainer::iterator it = allBoundaryDofs.begin();
391
392
393
	 it != allBoundaryDofs.end(); ++it)
      if (primalDofMap[feSpace].isSet(**it) == false)
	dualDofMap[feSpace].insertRankDof(**it);
394

395
396
    dualDofMap[feSpace].update(false);
    dualDofMap[feSpace].addOffset(rStartB + nRankInteriorDofs);
397

398
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
399
400
	dualDofMap[feSpace].nRankDofs, 
	dualDofMap[feSpace].nOverallDofs);
401
402

    nLocalDuals = dualDofMap[feSpace].size();
403
404
405
406
407
408
409
  }

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

410
411
412
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

413
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
414
    lagrangeMap.addFeSpace(feSpace);
415
416
417
418
419

    int nRankLagrange = 0;
    DofMapping& dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
420
	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
421
	int degree = boundaryDofRanks[it->first].size();
422
423
424
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
425
426
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
    lagrangeMap[feSpace].update();
427
    
428
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
429
430
	lagrangeMap[feSpace].nRankDofs, 
	lagrangeMap[feSpace].nOverallDofs);
431
432


433
    // === Communicate lagrangeMap to all other ranks. ===
434
435

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
436

437
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
438
	 !it.end(); it.nextRank()) {
439
      for (; !it.endDofIter(); it.nextDof())
440
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
441
	  DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
442
	  stdMpi.getSendData(it.getRank()).push_back(d);
443
	}
444
    }
445

446
447
    stdMpi.updateSendDataSize();

448
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
449
	 !it.end(); it.nextRank()) {
450
      bool recvData = false;
451
      for (; !it.endDofIter(); it.nextDof())
452
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
453
454
455
456
457
	  recvData = true;
	  break;
	}
	  
      if (recvData)
458
	stdMpi.recv(it.getRank());
459
460
461
462
    }

    stdMpi.startCommunication();

463
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
464
	 !it.end(); it.nextRank()) {
465
      int counter = 0;
466
467
468
      for (; !it.endDofIter(); it.nextDof()) {
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
469
	  lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
470
471
	}
      }
472
    }
473
474
475
476
477
  }


  void PetscSolverFeti::createIndexB()
  {
478
    FUNCNAME("PetscSolverFeti::createIndexB()");
479

480
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
481
    localDofMap.addFeSpace(feSpace);
482
    DOFAdmin* admin = feSpace->getAdmin();
483
484
485
486

    // === 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.                                 ===
487
488
489

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
490
491
      if (admin->isDofFree(i) == false && 
	  primalDofMap[feSpace].isSet(i) == false && 
492
493
494
495
	  dualDofMap[feSpace].isSet(i) == false) {
	localDofMap[feSpace].insertRankDof(i);
	nLocalInterior++;
      }
496
497
    }

498
499
    TEST_EXIT_DBG(nLocalInterior + 
		  primalDofMap[feSpace].size() + 
500
		  nLocalDuals == 
501
502
503
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

504
505
506

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

507
508
    for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
509
510
      localDofMap[feSpace].insert(it->first, it->second - rStartB);
      //      localDofMap[feSpace].insertRankDof(it->first);
511
512
513
  }


514
  void PetscSolverFeti::createMatLagrange()
515
516
517
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

518
519
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

520
521
    // === Create distributed matrix for Lagrange constraints. ===

522
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
523
		    lagrangeMap[feSpace].nRankDofs * nComponents, 
524
		    nRankB * nComponents,
525
		    lagrangeMap[feSpace].nOverallDofs * nComponents, 
526
		    nOverallB * nComponents,
527
528
529
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

530
531
532
533
534
535
536
    // === 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.                                                     ===

537
538
539
540
    DofMapping &dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");
541

542
      // Global index of the first Lagrange constriant for this node.
543
      int index = lagrangeMap[feSpace][it->first];
544
      // Copy set of all ranks that contain this dual node.
545
546
      vector<int> W(boundaryDofRanks[it->first].begin(), 
		    boundaryDofRanks[it->first].end());
547
      // Number of ranks that contain this dual node.
548
549
550
551
      int degree = W.size();

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
552
553
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
554
555
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
556
	      int colIndex = it->second * nComponents + k;
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
	      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);
  }


573
  void PetscSolverFeti::createSchurPrimalKsp()
574
575
576
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

577
578
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
582
583
584
585
586
587
588
589
590
      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, 
591
		   primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
592
593
594
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
595
596
		     primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nRankDofs * nComponents,
		     primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
597
598
599
600
601
602
603
604
605
606
		     &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 {
607
608
      MSG("Create direct schur primal solver!\n");

609
610
611
612
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

613
614
      int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents;
      int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
615
616
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
617

Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
620
621
622
623
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
624

625
626
627
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
628

629
630
631
632
633
634
635
636
637
638
      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
639
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
640
641
642
643
644
645
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
646
		primalDofMap[feSpace].size() * nComponents)("Should not happen!\n");
647
648
649
650
651
652
653
654
655
656
657
658
659
660
      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
661
	PetscScalar *tmpValues;
662
663
	VecGetArray(tmpVec, &tmpValues);
	for (int i  = 0; i < nRankB * nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
667
668
	  MatSetValue(matBPi, 
		      rStartB * nComponents + i,
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
669
670
671
672
673
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
674
675
676
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
677
678
679
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
680
      MatDestroy(&matBPi);
681
682
683
684
685
686
687
688
689
690

      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
691

692
693
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
694
    }
695
696
697
698
699
700
701
  }


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

702
703
704
705
706
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
707

708
709
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
710

711
712
713
714
715
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
716
717
718
719
720
721
722
  }


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

723
724
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

725
726
    // === Create FETI-DP solver object. ===

727
728
729
    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
730
    fetiData.fetiSolver = this;
731
    fetiData.ksp_schur_primal = &ksp_schur_primal;
732

733
734
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
737
738
		 lagrangeMap[feSpace].nRankDofs * nComponents, 
		 lagrangeMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
739
		 &(fetiData.tmp_vec_lagrange));
740
    VecCreateMPI(PETSC_COMM_WORLD, 
741
		 primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
742
		 &(fetiData.tmp_vec_primal));
743
744

    MatCreateShell(PETSC_COMM_WORLD,
745
746
747
748
		   lagrangeMap[feSpace].nRankDofs * nComponents, 
		   lagrangeMap[feSpace].nRankDofs * nComponents,
		   lagrangeMap[feSpace].nOverallDofs * nComponents, 
		   lagrangeMap[feSpace].nOverallDofs * nComponents,
749
		   &fetiData, &mat_feti);
750
751
752
753
754
755
756
    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);
757
758


759
    // === Create FETI-DP preconditioner object. ===
760

761
762
763
764
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
765

766
767
768
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
769
770
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
771
772
773
774
775
776
777
778
779
780
      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;
      
781
782
783
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiDirichletPreconData.tmp_vec_b));      
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
      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;

799
800
801
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiLumpedPreconData.tmp_vec_b));
802
803
804
805
806
807
808
809
      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);
      
810
811
      break;
    default:
812
813
      break;
    }
814
815
816
817
818
819
820
  }
  

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

821
822
    // === Destroy FETI-DP solver object. ===

823
824
825
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
826
    fetiData.fetiSolver = NULL;
827
    fetiData.ksp_schur_primal = PETSC_NULL;
828

Thomas Witkowski's avatar
Thomas Witkowski committed
829
830
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
831
    VecDestroy(&fetiData.tmp_vec_primal);
832
833
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
834
835


836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
    // === 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;
864
865
    default:
      break;
866
    }
867
868
869
870
871
872
873
874
875
  }


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

876
877
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

878
    // === Get local part of the solution for B variables. ===
879
880
881
882
883

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


884
885
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
886
887
    
    vector<PetscInt> globalIsIndex, localIsIndex;
888
889
    globalIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
    localIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
890
891
892

    {
      int counter = 0;
893
894
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
	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);

925
926
927
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
928
929
930
931

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

932
    // === And copy from PETSc local vectors to the DOF vectors. ===
933
934
935
936

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

937
938
      for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin();
	   it != localDofMap[feSpace].getMap().end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
939
	int petscIndex = it->second * nComponents + i;
940
941
942
943
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
944
945
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
946
947
948
949
950
951
952
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
953
    VecDestroy(&local_sol_primal);
954
955
956
  }


957
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
958
959
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
960

961
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
962
    nComponents = mat->getSize();
963
964
965

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

966
967
    primalDofMap.setMpiComm(mpiComm);
    dualDofMap.setMpiComm(mpiComm);
968
969
    lagrangeMap.setMpiComm(mpiComm);
    localDofMap.setMpiComm(mpiComm);
970
971
    updateDofData();

972
973
974
975
    // === Create matrices for the FETI-DP method. ===

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
976
977
    int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents;
    int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents;