PetscSolverFeti.cc 44.7 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);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
254
255
      else
  	primalDofMap[feSpace].insert(*it);
256

Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
257
258
259
    primalDofMap[feSpace].setOverlap(true);
    primalDofMap[feSpace].setDofComm(meshDistributor->getSendDofs(),
  				     meshDistributor->getRecvDofs());
260
    primalDofMap[feSpace].update();
261

262
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
263
	primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
264

265
    TEST_EXIT_DBG(primals.size() == primalDofMap[feSpace].size())
266
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
267
       primals.size(), primalDofMap[feSpace].size());
268
269
270
  }


271
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
272
273
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
274

275
276
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
277
278
279

    boundaryDofRanks.clear();

280
281
282
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
283
	if (!isPrimal(feSpace, it.getDofIndex())) {
284
285
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
286
287
	}
      }
288
	
289
290
291
292

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

293
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
294
295
296
297

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
298
	if (!isPrimal(feSpace, it.getDofIndex()))
299
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
300
301
302

    stdMpi.updateSendDataSize();

303
304
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
305
      bool recvFromRank = false;
306
      for (; !it.endDofIter(); it.nextDof()) {
307
	if (!isPrimal(feSpace, it.getDofIndex())) {
308
309
310
	  recvFromRank = true;
	  break;
	}
311
      }
312
313

      if (recvFromRank)
314
	stdMpi.recv(it.getRank());
315
    }
316

317
318
    stdMpi.startCommunication();

319
320
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
321
      int i = 0;
322
      for (; !it.endDofIter(); it.nextDof())
323
	if (!isPrimal(feSpace, it.getDofIndex()))
324
325
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
326
327
328
329
330
    }

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

    DofContainer allBoundaryDofs;
331
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
332
333

    for (DofContainer::iterator it = allBoundaryDofs.begin();
334
	 it != allBoundaryDofs.end(); ++it)
335
      if (!isPrimal(feSpace, **it))
336
	dualDofMap[feSpace].insertRankDof(**it);
337

338
    dualDofMap[feSpace].update(false);
339

340
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
341
342
	dualDofMap[feSpace].nRankDofs, 
	dualDofMap[feSpace].nOverallDofs);
343
344
345
  }

  
346
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
347
348
349
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

350
351
352
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

353
354
355
356
    int nRankLagrange = 0;
    DofMapping& dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
357
	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
358
	int degree = boundaryDofRanks[it->first].size();
359
360
361
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
362
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
363
    lagrangeMap[feSpace].update();
364
    
365
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
366
367
	lagrangeMap[feSpace].nRankDofs, 
	lagrangeMap[feSpace].nOverallDofs);
368
369


370
    // === Communicate lagrangeMap to all other ranks. ===
371
372

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

374
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
375
	 !it.end(); it.nextRank()) {
376
      for (; !it.endDofIter(); it.nextDof())
377
	if (!isPrimal(feSpace, it.getDofIndex())) {
378
	  DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
379
	  stdMpi.getSendData(it.getRank()).push_back(d);
380
	}
381
    }
382

383
384
    stdMpi.updateSendDataSize();

385
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
386
	 !it.end(); it.nextRank()) {
387
      bool recvData = false;
388
      for (; !it.endDofIter(); it.nextDof())
389
	if (!isPrimal(feSpace, it.getDofIndex())) {
390
391
392
393
394
	  recvData = true;
	  break;
	}
	  
      if (recvData)
395
	stdMpi.recv(it.getRank());
396
397
398
399
    }

    stdMpi.startCommunication();

400
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
401
	 !it.end(); it.nextRank()) {
402
      int counter = 0;
403
      for (; !it.endDofIter(); it.nextDof()) {
404
	if (!isPrimal(feSpace, it.getDofIndex())) {
405
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
406
	  lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
407
408
	}
      }
409
    }
410
411
412
  }


413
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
414
  {
415
    FUNCNAME("PetscSolverFeti::createIndexB()");
416

417
    DOFAdmin* admin = feSpace->getAdmin();
418
419
420
421

    // === 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.                                 ===
422
423
424

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
425
      if (admin->isDofFree(i) == false && 
426
	  isPrimal(feSpace, i) == false &&
427
428
429
430
	  dualDofMap[feSpace].isSet(i) == false) {
	localDofMap[feSpace].insertRankDof(i);
	nLocalInterior++;
      }
431
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
432
    
433
434
    localDofMap[feSpace].update(false);

Thomas Witkowski's avatar
Thomas Witkowski committed
435
    TEST_EXIT_DBG(nLocalInterior + primalDofMap[feSpace].size() + 
436
		  dualDofMap[feSpace].size() == 
437
438
439
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

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

442
443
    for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
444
445
446
447
448
449
      localDofMap[feSpace].insertRankDof(it->first);

    localDofMap[feSpace].update(false);

    dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs + 
				  nLocalInterior);
450
451
452
  }


453
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
454
455
456
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

457
458
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

459
460
    // === Create distributed matrix for Lagrange constraints. ===

461
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
462
463
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
464
465
466
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

467
468
469
470
471
472
473
    // === 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.                                                     ===

474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
    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;
493
	      int colIndex = it->second * nComponents + k;
494
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
495
496
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
497
498
	    }
	  }
499
500

	  index++;
501
502
503
504
505
506
507
508
509
	}
      }
    }

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


510
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
511
512
513
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

514
515
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
519
520
521
522
523
524
      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, 
525
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
526
527
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
528
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
529
530
531
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
532
533
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
536
537
538
539
540
541
542
543
		     &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 {
544
545
      MSG("Create direct schur primal solver!\n");

546
547
548
549
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

550
551
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
552
553
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
554

Thomas Witkowski's avatar
Thomas Witkowski committed
555
556
557
558
559
560
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
561

562
563
564
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
565

566
567
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

568
569
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
570
571
572
573
574
575
	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
576
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
577
578
579
580
581
582
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
583
584
		primalDofMap[feSpace].size() * nComponents)
	("Should not happen!\n");
585
586
587
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
588
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
589
590
591
592
593
594
595
596
597
598

 	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
599
	PetscScalar *tmpValues;
600
	VecGetArray(tmpVec, &tmpValues);
601
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
602
	  MatSetValue(matBPi, 
603
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
604
605
606
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
607
608
609
610
611
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
612
613
614
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
615
616
617
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
618
      MatDestroy(&matBPi);
619
620
621
622
623
624
625
626
627
628

      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
629

630
631
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
632
    }
633
634
635
636
637
638
639
  }


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

640
641
642
643
644
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
645

646
647
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
648

649
650
651
652
653
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
654
655
656
  }


657
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
658
659
660
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

661
662
    // === Create FETI-DP solver object. ===

663
664
665
    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
666
    fetiData.fetiSolver = this;
667
    fetiData.ksp_schur_primal = &ksp_schur_primal;
668

669
    VecCreateMPI(PETSC_COMM_WORLD, 
670
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
671
672
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
673
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
674
		 &(fetiData.tmp_vec_lagrange));
675
    VecCreateMPI(PETSC_COMM_WORLD, 
676
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
677
		 &(fetiData.tmp_vec_primal));
678
679

    MatCreateShell(PETSC_COMM_WORLD,
680
681
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
682
		   &fetiData, &mat_feti);
683
684
685
686
687
688
689
    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);
690
691


692
    // === Create FETI-DP preconditioner object. ===
693

694
695
696
697
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
698

699
700
701
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
702
703
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
704
705
706
707
708
709
710
711
712
713
      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;
      
714
      VecCreateMPI(PETSC_COMM_WORLD, 
715
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
716
		   &(fetiDirichletPreconData.tmp_vec_b));      
717
718
719
720
721
722
      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));
723
724
725
726
727
728
729
730
731
732
733
734
      
      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;

735
      VecCreateMPI(PETSC_COMM_WORLD, 
736
737
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
738
		   &(fetiLumpedPreconData.tmp_vec_b));
739
740
741
742
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
743
744
745
746
747
748

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
749
750
      break;
    default:
751
752
      break;
    }
753
754
755
756
757
758
759
  }
  

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

760
761
    // === Destroy FETI-DP solver object. ===

762
763
764
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
765
    fetiData.fetiSolver = NULL;
766
    fetiData.ksp_schur_primal = PETSC_NULL;
767

Thomas Witkowski's avatar
Thomas Witkowski committed
768
769
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
770
    VecDestroy(&fetiData.tmp_vec_primal);
771
772
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
773
774


775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
    // === 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;
803
804
    default:
      break;
805
    }
806
807
808
809
810
811
812
813
814
  }


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

815
816
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

817
    // === Get local part of the solution for B variables. ===
818
819
820
821
822

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


823
824
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
825
826
    
    vector<PetscInt> globalIsIndex, localIsIndex;
827
828
    globalIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
    localIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
829
830
831

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

864
865
866
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
867
868
869
870

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

871
    // === And copy from PETSc local vectors to the DOF vectors. ===
872
873
874
875

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

876
877
      for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin();
	   it != localDofMap[feSpace].getMap().end(); ++it) {
878
	int petscIndex = localDofMap.mapLocal(it->first, i);
879
880
881
882
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
883
884
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
885
886
887
888
889
890
891
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
892
    VecDestroy(&local_sol_primal);
893
894
895
  }


896
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
897
898
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
899

900
901
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);

902
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
903
    nComponents = mat->getSize();
904
905
906

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

907
908
    primalDofMap.setMpiComm(mpiComm);
    dualDofMap.setMpiComm(mpiComm);
909
910
    lagrangeMap.setMpiComm(mpiComm);
    localDofMap.setMpiComm(mpiComm);
911
912
913
914
915
916

    primalDofMap.setFeSpaces(feSpaces);
    dualDofMap.setFeSpaces(feSpaces);
    lagrangeMap.setFeSpaces(feSpaces);
    localDofMap.setFeSpaces(feSpaces);

917
918
    updateDofData();

919
920
921
922
923
    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();

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

926
927
    int nRowsRankB = localDofMap.getRankDofs();
    int nRowsOverallB = localDofMap.getOverallDofs();
928
929
930
    int nRowsRankPrimal = primalDofMap.getRankDofs();
    int nRowsOverallPrimal = primalDofMap.getOverallDofs();
    int nRowsDual = dualDofMap.getRankDofs();
931
932
    int nRowsInterior = nRowsRankB - nRowsDual;

Thomas Witkowski's avatar
Thomas Witkowski committed
933
934
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
935
936

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
937
938
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
939
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
940
941

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
942
943
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
944
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
945
946

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
947
948
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
949
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
950

951

952
953
954
955
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,