PetscSolverFeti.cc 54.6 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
224
225
      createPrimals(feSpace);      
      createDuals(feSpace);
      createIndexB(feSpace);
    }
226
227

    createLagrange(meshDistributor->getFeSpaces());
228
229
230
  }


231
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
232
  {
233
    FUNCNAME("PetscSolverFeti::createPrimals()");  
234

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
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
253
254
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
255

256
    primalDofMap[feSpace].update();
257

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

262
263
264
265

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

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

278
279
    stdMpi.updateSendDataSize();

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

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

295
296
    stdMpi.startCommunication();

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

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


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

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

    boundaryDofRanks.clear();

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

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

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

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
359
	stdMpi.recv(it.getRank());
360
    }
361

362
363
    stdMpi.startCommunication();

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

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

    DofContainer allBoundaryDofs;
376
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
377
378

    for (DofContainer::iterator it = allBoundaryDofs.begin();
379
380
381
	 it != allBoundaryDofs.end(); ++it)
      if (primalDofMap[feSpace].isSet(**it) == false)
	dualDofMap[feSpace].insertRankDof(**it);
382

383
    dualDofMap[feSpace].update(false);
384

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

  
390
  void PetscSolverFeti::createLagrange(vector<const FiniteElemSpace*> &feSpaces)
391
392
393
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

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

397
398

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


415
    // === Communicate lagrangeMap to all other ranks. ===
416
417

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

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

428
429
    stdMpi.updateSendDataSize();

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

    stdMpi.startCommunication();

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


459
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
460
  {
461
    FUNCNAME("PetscSolverFeti::createIndexB()");
462

463
    DOFAdmin* admin = feSpace->getAdmin();
464
465
466
467

    // === 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.                                 ===
468
469
470

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
471
472
      if (admin->isDofFree(i) == false && 
	  primalDofMap[feSpace].isSet(i) == false && 
473
474
475
476
	  dualDofMap[feSpace].isSet(i) == false) {
	localDofMap[feSpace].insertRankDof(i);
	nLocalInterior++;
      }
477
478
    }

479
480
    localDofMap[feSpace].update(false);

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

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

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

    localDofMap[feSpace].update(false);

495
#if 0
496
497
    dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs + 
				  nLocalInterior);
498
#endif
499
500
501
  }


502
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
503
504
505
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

506
507
    // === Create distributed matrix for Lagrange constraints. ===

508
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
509
		    lagrangeMap.getRankDofs(feSpaces),
510
		    localDofMap.getRankDofs(),
511
		    lagrangeMap.getOverallDofs(feSpaces),
512
		    localDofMap.getOverallDofs(),	
513
514
515
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

516
517
518
519
520
521
522
    // === 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.                                                     ===

523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
    for (unsigned int feIdx = 0; feIdx < feSpaces.size(); feIdx++) {
      const FiniteElemSpace *feSpace = feSpaces[feIdx];

      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 constraintMatIndex = lagrangeMap.mapGlobal(it->first, feIdx);
	// 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) {
#if 1
	      // Get global index of the dual variable
	      int colIndex = localDofMap.mapGlobal(it->first, feIdx);
#else
546
	      int colIndex = it->second * nComponents + k;
547
#endif
548
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
549
550
	      MatSetValue(mat_lagrange, constraintMatIndex, colIndex, value, 
			  INSERT_VALUES);	      
551
	    }
552
553
	    
	    constraintMatIndex++;
554
555
556
557
558
559
560
561
562
563
	  }
	}
      }
    }

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


564
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
565
566
567
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

568
569
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
573
574
575
576
577
578
      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, 
579
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
580
581
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
582
		   primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
583
584
585
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
586
587
		     primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nRankDofs * nComponents,
		     primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
588
589
590
591
592
593
594
595
596
597
		     &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 {
598
599
      MSG("Create direct schur primal solver!\n");

600
601
602
603
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

604
605
      int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents;
      int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents;
606
607
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
608

Thomas Witkowski's avatar
Thomas Witkowski committed
609
610
611
612
613
614
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
615

616
617
618
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
619

620
621
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

622
623
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
624
625
626
627
628
629
	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
630
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
631
632
633
634
635
636
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
637
		primalDofMap[feSpace].size() * nComponents)("Should not happen!\n");
638
639
640
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
641
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
642
643
644
645
646
647
648
649
650
651

 	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
652
	PetscScalar *tmpValues;
653
	VecGetArray(tmpVec, &tmpValues);
654
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
655
	  MatSetValue(matBPi, 
656
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
659
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
660
661
662
663
664
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
667
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
668
669
670
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
671
      MatDestroy(&matBPi);
672
673
674
675
676
677
678
679
680
681

      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
682

683
684
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
685
    }
686
687
688
689
690
691
692
  }


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

693
694
695
696
697
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
698

699
700
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
701

702
703
704
705
706
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
707
708
709
  }


710
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
711
712
713
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

714
715
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

716
717
    // === Create FETI-DP solver object. ===

718
719
720
    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
721
    fetiData.fetiSolver = this;
722
    fetiData.ksp_schur_primal = &ksp_schur_primal;
723

724
    VecCreateMPI(PETSC_COMM_WORLD, 
725
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
726
727
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
728
729
		 lagrangeMap[feSpace].nRankDofs * nComponents, 
		 lagrangeMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
730
		 &(fetiData.tmp_vec_lagrange));
731
    VecCreateMPI(PETSC_COMM_WORLD, 
732
733
		 primalDofMap[feSpace].nRankDofs * nComponents, 
		 primalDofMap[feSpace].nOverallDofs * nComponents,
734
		 &(fetiData.tmp_vec_primal));
735
736

    MatCreateShell(PETSC_COMM_WORLD,
737
738
739
740
		   lagrangeMap[feSpace].nRankDofs * nComponents, 
		   lagrangeMap[feSpace].nRankDofs * nComponents,
		   lagrangeMap[feSpace].nOverallDofs * nComponents, 
		   lagrangeMap[feSpace].nOverallDofs * nComponents,
741
		   &fetiData, &mat_feti);
742
743
744
745
746
747
748
    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);
749
750


751
    // === Create FETI-DP preconditioner object. ===
752

753
754
755
756
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
757

758
759
760
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
761
762
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
763
764
765
766
767
768
769
770
771
772
      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;
      
773
      VecCreateMPI(PETSC_COMM_WORLD, 
774
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
775
		   &(fetiDirichletPreconData.tmp_vec_b));      
776
777
778
779
780
781
      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));
782
783
784
785
786
787
788
789
790
791
792
793
      
      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;

794
      VecCreateMPI(PETSC_COMM_WORLD, 
795
796
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
797
		   &(fetiLumpedPreconData.tmp_vec_b));
798
799
800
801
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
802
803
804
805
806
807

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
808
809
      break;
    default:
810
811
      break;
    }
812
813
814
815
816
817
818
  }
  

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

819
820
    // === Destroy FETI-DP solver object. ===

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

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


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
    // === 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;
862
863
    default:
      break;
864
    }
865
866
867
868
869
870
871
872
873
  }


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

874
875
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

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


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

    {
      int counter = 0;
891
892
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
893
894
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
	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);

923
924
925
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
926
927
928
929

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

930
    // === And copy from PETSc local vectors to the DOF vectors. ===
931
932
933
934

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

935
936
      for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin();
	   it != localDofMap[feSpace].getMap().end(); ++it) {
937
938
939
940

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

944
945
946
947
	dofVec[it->first] = localSolB[petscIndex];
      }

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

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
957
    VecDestroy(&local_sol_primal);
958
959
960
  }


961
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
962
963
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
964