PetscSolverFeti.cc 46.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
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
221

222
223
      createPrimals(feSpace);      
      createDuals(feSpace);
224
      createLagrange(feSpace);      
225
226
      createIndexB(feSpace);
    }
227
228
229
  }


230
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
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
    DofContainerSet& vertices = 
241
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
242
243
244
245
    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
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
252
253
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
254

255
    primalDofMap[feSpace].update();
256

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

261
262
263
264

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

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

277
278
    stdMpi.updateSendDataSize();

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

      if (recvFromRank) 
291
	stdMpi.recv(it.getRank());
292
    }
293

294
295
    stdMpi.startCommunication();

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

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


314
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
315
316
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
317

318
319
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
320
321
322

    boundaryDofRanks.clear();

323
324
325
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
326
	if (!isPrimal(feSpace, it.getDofIndex())) {
327
328
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
329
330
	}
      }
331
	
332
333
334
335

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

336
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
337
338
339
340

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
341
	if (!isPrimal(feSpace, it.getDofIndex()))
342
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
343
344
345

    stdMpi.updateSendDataSize();

346
347
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
348
      bool recvFromRank = false;
349
      for (; !it.endDofIter(); it.nextDof()) {
350
	if (!isPrimal(feSpace, it.getDofIndex())) {
351
352
353
	  recvFromRank = true;
	  break;
	}
354
      }
355
356

      if (recvFromRank)
357
	stdMpi.recv(it.getRank());
358
    }
359

360
361
    stdMpi.startCommunication();

362
363
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
364
      int i = 0;
365
      for (; !it.endDofIter(); it.nextDof())
366
	if (!isPrimal(feSpace, it.getDofIndex()))
367
368
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
369
370
371
372
373
    }

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

    DofContainer allBoundaryDofs;
374
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
375
376

    for (DofContainer::iterator it = allBoundaryDofs.begin();
377
	 it != allBoundaryDofs.end(); ++it)
378
      if (!isPrimal(feSpace, **it))
379
	dualDofMap[feSpace].insertRankDof(**it);
380

381
    dualDofMap[feSpace].update(false);
382

383
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
384
385
	dualDofMap[feSpace].nRankDofs, 
	dualDofMap[feSpace].nOverallDofs);
386
387
388
  }

  
389
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
390
391
392
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

393
394
395
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

396
397
398
399
    int nRankLagrange = 0;
    DofMapping& dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
400
	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
401
	int degree = boundaryDofRanks[it->first].size();
402
403
404
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
405
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
406
    lagrangeMap[feSpace].update();
407
    
408
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
409
410
	lagrangeMap[feSpace].nRankDofs, 
	lagrangeMap[feSpace].nOverallDofs);
411
412


413
    // === Communicate lagrangeMap to all other ranks. ===
414
415

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

417
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
418
	 !it.end(); it.nextRank()) {
419
      for (; !it.endDofIter(); it.nextDof())
420
	if (!isPrimal(feSpace, it.getDofIndex())) {
421
	  DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
422
	  stdMpi.getSendData(it.getRank()).push_back(d);
423
	}
424
    }
425

426
427
    stdMpi.updateSendDataSize();

428
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
429
	 !it.end(); it.nextRank()) {
430
      bool recvData = false;
431
      for (; !it.endDofIter(); it.nextDof())
432
	if (!isPrimal(feSpace, it.getDofIndex())) {
433
434
435
436
437
	  recvData = true;
	  break;
	}
	  
      if (recvData)
438
	stdMpi.recv(it.getRank());
439
440
441
442
    }

    stdMpi.startCommunication();

443
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
444
	 !it.end(); it.nextRank()) {
445
      int counter = 0;
446
      for (; !it.endDofIter(); it.nextDof()) {
447
	if (!isPrimal(feSpace, it.getDofIndex())) {
448
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
449
	  lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
450
451
	}
      }
452
    }
453
454
455
  }


456
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
457
  {
458
    FUNCNAME("PetscSolverFeti::createIndexB()");
459

460
    DOFAdmin* admin = feSpace->getAdmin();
461
462
463
464

    // === 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.                                 ===
465
466
467

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
468
      if (admin->isDofFree(i) == false && 
469
	  isPrimal(feSpace, i) == false &&
470
471
472
473
	  dualDofMap[feSpace].isSet(i) == false) {
	localDofMap[feSpace].insertRankDof(i);
	nLocalInterior++;
      }
474
475
    }

476
477
    localDofMap[feSpace].update(false);

478
479
    TEST_EXIT_DBG(nLocalInterior + 
		  primalDofMap[feSpace].size() + 
480
		  dualDofMap[feSpace].size() == 
481
482
483
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

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

486
487
    for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
488
489
490
491
492
493
      localDofMap[feSpace].insertRankDof(it->first);

    localDofMap[feSpace].update(false);

    dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs + 
				  nLocalInterior);
494
495
496
  }


497
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
498
499
500
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

501
502
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

503
504
    // === Create distributed matrix for Lagrange constraints. ===

505
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
506
507
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
508
509
510
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

511
512
513
514
515
516
517
    // === 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.                                                     ===

518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
    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");

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

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
537
	      int colIndex = it->second * nComponents + k;
538
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
539
540
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
541
542
	    }
	  }
543
544

	  index++;
545
546
547
548
549
550
551
552
553
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
  }


554
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
555
556
557
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

558
559
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
563
564
565
566
567
568
      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, 
569
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
570
571
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
572
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
573
574
575
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
576
577
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
578
579
580
581
582
583
584
585
586
587
		     &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 {
588
589
      MSG("Create direct schur primal solver!\n");

590
591
592
593
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

594
595
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
596
597
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
598

Thomas Witkowski's avatar
Thomas Witkowski committed
599
600
601
602
603
604
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
605

606
607
608
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
609

610
611
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

612
613
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
614
615
616
617
618
619
	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
620
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
621
622
623
624
625
626
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
627
628
		primalDofMap[feSpace].size() * nComponents)
	("Should not happen!\n");
629
630
631
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
632
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
633
634
635
636
637
638
639
640
641
642

 	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
643
	PetscScalar *tmpValues;
644
	VecGetArray(tmpVec, &tmpValues);
645
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
646
	  MatSetValue(matBPi, 
647
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
648
649
650
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
651
652
653
654
655
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
658
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
659
660
661
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
662
      MatDestroy(&matBPi);
663
664
665
666
667
668
669
670
671
672

      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
673

674
675
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
676
    }
677
678
679
680
681
682
683
  }


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

684
685
686
687
688
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
689

690
691
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
692

693
694
695
696
697
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
698
699
700
  }


701
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
702
703
704
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

705
706
    // === Create FETI-DP solver object. ===

707
708
709
    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
710
    fetiData.fetiSolver = this;
711
    fetiData.ksp_schur_primal = &ksp_schur_primal;
712

713
    VecCreateMPI(PETSC_COMM_WORLD, 
714
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
715
716
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
717
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
718
		 &(fetiData.tmp_vec_lagrange));
719
    VecCreateMPI(PETSC_COMM_WORLD, 
720
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
721
		 &(fetiData.tmp_vec_primal));
722
723

    MatCreateShell(PETSC_COMM_WORLD,
724
725
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
726
		   &fetiData, &mat_feti);
727
728
729
730
731
732
733
    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);
734
735


736
    // === Create FETI-DP preconditioner object. ===
737

738
739
740
741
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
742

743
744
745
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
746
747
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
748
749
750
751
752
753
754
755
756
757
      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;
      
758
      VecCreateMPI(PETSC_COMM_WORLD, 
759
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
760
		   &(fetiDirichletPreconData.tmp_vec_b));      
761
762
763
764
765
766
      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));
767
768
769
770
771
772
773
774
775
776
777
778
      
      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;

779
      VecCreateMPI(PETSC_COMM_WORLD, 
780
781
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
782
		   &(fetiLumpedPreconData.tmp_vec_b));
783
784
785
786
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
787
788
789
790
791
792

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
793
794
      break;
    default:
795
796
      break;
    }
797
798
799
800
801
802
803
  }
  

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

804
805
    // === Destroy FETI-DP solver object. ===

806
807
808
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
809
    fetiData.fetiSolver = NULL;
810
    fetiData.ksp_schur_primal = PETSC_NULL;
811

Thomas Witkowski's avatar
Thomas Witkowski committed
812
813
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
814
    VecDestroy(&fetiData.tmp_vec_primal);
815
816
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
817
818


819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
    // === 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;
847
848
    default:
      break;
849
    }
850
851
852
853
854
855
856
857
858
  }


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

859
860
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

861
    // === Get local part of the solution for B variables. ===
862
863
864
865
866

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


867
868
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
869
870
    
    vector<PetscInt> globalIsIndex, localIsIndex;
871
872
    globalIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
    localIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
873
874
875

    {
      int counter = 0;
876
877
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
	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);

908
909
910
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
911
912
913
914

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

915
    // === And copy from PETSc local vectors to the DOF vectors. ===
916
917
918
919

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

920
921
      for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin();
	   it != localDofMap[feSpace].getMap().end(); ++it) {
922
923
924
925

#if 1
	int petscIndex = localDofMap.mapLocal(it->first, i);
#else	  
Thomas Witkowski's avatar
Thomas Witkowski committed
926
	int petscIndex = it->second * nComponents + i;
927
928
#endif

929
930
931
932
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
933
934
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
935
936
937
938
939
940
941
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
942
    VecDestroy(&local_sol_primal);
943
944
945
  }


946
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
947
948
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
949

950
951
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);

952
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
953
    nComponents = mat->getSize();
954
955
956

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

Thomas Witkowski's avatar