PetscSolverFeti.cc 43.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
221
222
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createPrimals(feSpace);      
      createDuals(feSpace);
223
      createLagrange(feSpace);      
224
225
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245

    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();

    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);

      MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
      
      MSG("nRankDuals = %d   nOverallDuals = %d\n",
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);

      MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
    }
246
247
248
  }


249
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
250
  {
251
    FUNCNAME("PetscSolverFeti::createPrimals()");  
252

253
254
255
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

256
257
258
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
259
    DofContainerSet& vertices = 
260
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
261
262
263
264
    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);
265

266
267
268
269

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

270
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
271
272
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
273
274
      else
  	primalDofMap[feSpace].insert(*it);
275

Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
276
277
278
    primalDofMap[feSpace].setOverlap(true);
    primalDofMap[feSpace].setDofComm(meshDistributor->getSendDofs(),
  				     meshDistributor->getRecvDofs());
279
280
281
  }


282
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
283
284
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
285

286
287
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
288
289
290

    boundaryDofRanks.clear();

291
292
293
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
294
	if (!isPrimal(feSpace, it.getDofIndex())) {
295
296
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
297
298
	}
      }
299
	
300
301
302
303

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

304
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
305
306
307
308

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
309
	if (!isPrimal(feSpace, it.getDofIndex()))
310
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
311
312
313

    stdMpi.updateSendDataSize();

314
315
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
316
      bool recvFromRank = false;
317
      for (; !it.endDofIter(); it.nextDof()) {
318
	if (!isPrimal(feSpace, it.getDofIndex())) {
319
320
321
	  recvFromRank = true;
	  break;
	}
322
      }
323
324

      if (recvFromRank)
325
	stdMpi.recv(it.getRank());
326
    }
327

328
329
    stdMpi.startCommunication();

330
331
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
332
      int i = 0;
333
      for (; !it.endDofIter(); it.nextDof())
334
	if (!isPrimal(feSpace, it.getDofIndex()))
335
336
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
337
338
339
340
341
    }

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

    DofContainer allBoundaryDofs;
342
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
343
344

    for (DofContainer::iterator it = allBoundaryDofs.begin();
345
	 it != allBoundaryDofs.end(); ++it)
346
      if (!isPrimal(feSpace, **it))
347
	dualDofMap[feSpace].insertRankDof(**it);
348
349
350
  }

  
351
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
352
353
354
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

355
356
357
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

358
    int nRankLagrange = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
    map<DegreeOfFreedom, MultiIndex>& dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
361
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
362
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
363
	int degree = boundaryDofRanks[it->first].size();
364
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
      } else {
	lagrangeMap[feSpace].insert(it->first);
367
368
      }
    }
369
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
372
    lagrangeMap[feSpace].setOverlap(true);
    lagrangeMap[feSpace].setDofComm(meshDistributor->getSendDofs(),
				    meshDistributor->getRecvDofs());
373
374
375
  }


376
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
377
  {
378
    FUNCNAME("PetscSolverFeti::createIndexB()");
379

380
    DOFAdmin* admin = feSpace->getAdmin();
381
382
383
384

    // === 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.                                 ===
385
386
387

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
388
      if (admin->isDofFree(i) == false && 
389
	  isPrimal(feSpace, i) == false &&
390
	  dualDofMap[feSpace].isSet(i) == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
391
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
392
393
	nLocalInterior++;
      }
394
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
395
    
Thomas Witkowski's avatar
Thomas Witkowski committed
396
    TEST_EXIT_DBG(nLocalInterior + primalDofMap[feSpace].size() + 
397
		  dualDofMap[feSpace].size() == 
398
399
400
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
403
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualDofMap[feSpace].getMap().begin();
404
	 it != dualDofMap[feSpace].getMap().end(); ++it)
405
      localDofMap[feSpace].insertRankDof(it->first);
406
407
408
  }


409
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
410
411
412
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

413
414
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

415
416
    // === Create distributed matrix for Lagrange constraints. ===

417
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
418
419
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
420
421
422
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

423
424
425
426
427
428
429
    // === 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.                                                     ===

Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
    map<DegreeOfFreedom, MultiIndex> &dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
432
433
434
435
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");

      // Global index of the first Lagrange constriant for this node.
Thomas Witkowski's avatar
Thomas Witkowski committed
436
      int index = lagrangeMap[feSpace][it->first].local;
437
438
439
440
441
442
443
444
445
446
447
448
      // 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;
Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
451
452
	      //	      int colIndex = it->second * nComponents + k;
	      int colIndex = 
		(localDofMap[feSpace].rStartDofs +
		 localDofMap[feSpace][it->first].local) * nComponents + k;
453
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
454
455
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
456
457
	    }
	  }
458
459

	  index++;
460
461
462
463
464
465
466
467
468
	}
      }
    }

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


469
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
470
471
472
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

473
474
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
478
479
480
481
482
483
      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, 
484
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
487
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
488
489
490
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
491
492
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
493
494
495
496
497
498
499
500
501
502
		     &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 {
503
504
      MSG("Create direct schur primal solver!\n");

505
506
507
508
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

509
510
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
511
512
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
513

Thomas Witkowski's avatar
Thomas Witkowski committed
514
515
516
517
518
519
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
520

521
522
523
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
524

525
526
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

527
528
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
529
530
531
532
533
534
	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
535
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
536
537
538
539
540
      }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
541
      TEST_EXIT(mat_b_primal_cols.size() == primalDofMap.getLocalDofs())
542
	("Should not happen!\n");
543
544
545
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
546
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
547
548
549
550
551
552
553
554
555
556

 	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
557
	PetscScalar *tmpValues;
558
	VecGetArray(tmpVec, &tmpValues);
559
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
560
	  MatSetValue(matBPi, 
561
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
562
563
564
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
565
566
567
568
569
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
570
571
572
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
573
574
575
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
576
      MatDestroy(&matBPi);
577
578
579
580
581
582
583
584
585
586

      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
587

588
589
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
590
    }
591
592
593
594
595
596
597
  }


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

598
599
600
601
602
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
603

604
605
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
606

607
608
609
610
611
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
612
613
614
  }


615
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
616
617
618
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

619
620
    // === Create FETI-DP solver object. ===

621
622
623
    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
624
    fetiData.fetiSolver = this;
625
    fetiData.ksp_schur_primal = &ksp_schur_primal;
626

627
    VecCreateMPI(PETSC_COMM_WORLD, 
628
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
631
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
632
		 &(fetiData.tmp_vec_lagrange));
633
    VecCreateMPI(PETSC_COMM_WORLD, 
634
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
635
		 &(fetiData.tmp_vec_primal));
636
637

    MatCreateShell(PETSC_COMM_WORLD,
638
639
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
640
		   &fetiData, &mat_feti);
641
642
643
644
645
646
647
    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);
648
649


650
    // === Create FETI-DP preconditioner object. ===
651

652
653
654
655
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
656

657
658
659
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
660
661
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
662
663
664
665
666
667
668
669
670
671
      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;
      
672
      VecCreateMPI(PETSC_COMM_WORLD, 
673
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
674
		   &(fetiDirichletPreconData.tmp_vec_b));      
675
676
677
678
679
680
      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));
681
682
683
684
685
686
687
688
689
690
691
692
      
      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;

693
      VecCreateMPI(PETSC_COMM_WORLD, 
694
695
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
696
		   &(fetiLumpedPreconData.tmp_vec_b));
697
698
699
700
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
701
702
703
704
705
706

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
707
708
      break;
    default:
709
710
      break;
    }
711
712
713
714
715
716
717
  }
  

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

718
719
    // === Destroy FETI-DP solver object. ===

720
721
722
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
723
    fetiData.fetiSolver = NULL;
724
    fetiData.ksp_schur_primal = PETSC_NULL;
725

Thomas Witkowski's avatar
Thomas Witkowski committed
726
727
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
728
    VecDestroy(&fetiData.tmp_vec_primal);
729
730
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
731
732


733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
    // === 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;
761
762
    default:
      break;
763
    }
764
765
766
767
768
769
770
771
772
  }


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

773
774
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

775
    // === Get local part of the solution for B variables. ===
776
777
778
779
780

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


781
782
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
783
784
    
    vector<PetscInt> globalIsIndex, localIsIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
785
786
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());
787
788
789

    {
      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
790
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
791
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
792
	for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
793
	  globalIsIndex.push_back(it->second.local * nComponents + i);
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
	  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);

822
823
824
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
825
826
827
828

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

829
    // === And copy from PETSc local vectors to the DOF vectors. ===
830
831
832
833

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

Thomas Witkowski's avatar
Thomas Witkowski committed
834
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpace].getMap().begin();
835
	   it != localDofMap[feSpace].getMap().end(); ++it) {
836
	int petscIndex = localDofMap.mapLocal(it->first, i);
837
838
839
840
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
841
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
842
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
843
844
845
846
847
848
849
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
850
    VecDestroy(&local_sol_primal);
851
852
853
  }


854
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
855
856
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
857

858
859
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);

860
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
861
    nComponents = mat->getSize();
862
863
864

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

Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
867
868
    primalDofMap.init(mpiComm, feSpaces, true);
    dualDofMap.init(mpiComm, feSpaces, false);
    lagrangeMap.init(mpiComm, feSpaces, true);
    localDofMap.init(mpiComm, feSpaces, false);
869

870
871
    updateDofData();

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

874
875
    int nRowsRankB = localDofMap.getRankDofs();
    int nRowsOverallB = localDofMap.getOverallDofs();
876
877
878
    int nRowsRankPrimal = primalDofMap.getRankDofs();
    int nRowsOverallPrimal = primalDofMap.getOverallDofs();
    int nRowsDual = dualDofMap.getRankDofs();
879
880
    int nRowsInterior = nRowsRankB - nRowsDual;

Thomas Witkowski's avatar
Thomas Witkowski committed
881
882
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
883
884

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
885
886
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
887
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
888
889

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
890
891
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
892
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
893
894

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
895
896
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
897
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
898

899

900
901
902
903
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
904
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
905
906
907
908
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
909
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
910
911
912
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
913
		      nRowsInterior, nRowsDual, 30, PETSC_NULL,
914
915
916
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
917
		      nRowsDual, nRowsInterior, 30, PETSC_NULL,
918
919
		      &mat_duals_interior);
    }
920

921
922
    
    // === Prepare traverse of sequentially created matrices. ===
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

    vector<int> cols, colsOther;
    vector<double> values, valuesOther;
    cols.reserve(300);
    colsOther.reserve(300);
    values.reserve(300);
    valuesOther.reserve(300);

938
939
940
941
942
943
944
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

945
946
947
948
949
950
951
952
953
954
955
956
957
958
959

    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	if (!(*mat)[i][j])
	  continue;

	traits::col<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
	traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
	
	// Traverse all rows.
	for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()), 
	       cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
960

961
	  bool rowPrimal = primalDofMap[feSpace].isSet(*cursor);
962
  
963
964
	  cols.clear();
	  colsOther.clear();
965
	  values.clear();	  
966
	  valuesOther.clear();
967
968
969
970
971
972

	  colsLocal.clear();
	  colsLocalOther.clear();
	  valuesLocal.clear();
	  valuesLocalOther.clear();

973
974
975
976
977
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

978
	    bool colPrimal = primalDofMap[feSpace].isSet(col(*icursor));
979
980

	    if (colPrimal) {
981
	      if (rowPrimal) {
982
		cols.push_back(col(*icursor));
983
		values.push_back(value(*icursor));
984
	      } else {
985
		colsOther.push_back(col(*icursor));
Thomas Witkowski's avatar