PetscSolverFeti.cc 43.1 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
      schurPrimalSolver(0)
185
186
187
188
189
190
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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


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

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
217
   
218
219
220
221
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createPrimals(feSpace);      
      createDuals(feSpace);
222
      createLagrange(feSpace);      
223
224
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
225

226
227
228
229
230
    primalDofMap.setDofComm(meshDistributor->getSendDofs(),
			    meshDistributor->getRecvDofs());
    lagrangeMap.setDofComm(meshDistributor->getSendDofs(), 
			   meshDistributor->getRecvDofs());

Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
233
234
    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();
235
236
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
237
238
239
240
241
242
243
244

    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",
245
	  dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
246
247

      MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
248
	  lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
249
250
251
252

      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
253
    }
254
255
256
  }


257
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
258
  {
259
    FUNCNAME("PetscSolverFeti::createPrimals()");  
260

261
262
263
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

264
265
266
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
267
    DofContainerSet& vertices = 
268
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
269
270
271
272
    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);
273

274
275
276
277

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

278
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
279
280
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
281
282
      else
  	primalDofMap[feSpace].insert(*it);
283
284
285
  }


286
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
287
288
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
289

290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);

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

  
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

306
307
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
308

309
    boundaryDofRanks[feSpace].clear();
310

311
312
313
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
314
	if (!isPrimal(feSpace, it.getDofIndex())) {
315
316
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
317
318
	}
      }
319
	
320
321
322
323

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

324
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
325
326
327
328

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
329
	if (!isPrimal(feSpace, it.getDofIndex()))
330
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
331
332
333

    stdMpi.updateSendDataSize();

334
335
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
336
      bool recvFromRank = false;
337
      for (; !it.endDofIter(); it.nextDof()) {
338
	if (!isPrimal(feSpace, it.getDofIndex())) {
339
340
341
	  recvFromRank = true;
	  break;
	}
342
      }
343
344

      if (recvFromRank)
345
	stdMpi.recv(it.getRank());
346
    }
347

348
349
    stdMpi.startCommunication();

350
351
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
352
      int i = 0;
353
      for (; !it.endDofIter(); it.nextDof())
354
	if (!isPrimal(feSpace, it.getDofIndex()))
355
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
356
	    stdMpi.getRecvData(it.getRank())[i++];
357
358
359
    }


360
361
362
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


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

382
    DOFAdmin* admin = feSpace->getAdmin();
383
384
385
386

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

388
    int nLocalInterior = 0;    
389
    for (int i = 0; i < admin->getUsedSize(); i++) {
390
      if (admin->isDofFree(i) == false && 
391
	  isPrimal(feSpace, i) == false &&
392
	  dualDofMap[feSpace].isSet(i) == false) {
393
394
395
396
397
398
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);

	if (fetiPreconditioner == FETI_DIRICHLET)
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);

	nLocalInterior++;
399
      }
400
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
401
    
402
403
    // === And finally, add the global indicies of all dual nodes. ===

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


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

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

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

422
423
424
425
426
427
428
    // === 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
429
430
431
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
      map<DegreeOfFreedom, MultiIndex> &dualMap = 
	dualDofMap[feSpaces[k]].getMap();
432

Thomas Witkowski's avatar
Thomas Witkowski committed
433
434
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	   it != dualMap.end(); ++it) {
435
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
439
440
441
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
	int index = lagrangeMap.getMatIndex(k, it->first);
	
	// Copy set of all ranks that contain this dual node.
442
443
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
446
447
448
449
	// 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) {
450
	      int colIndex = localDofMap.getMatIndex(k, it->first);
451
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
452
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
453
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
454
	    index++;	      
455
456
457
458
459
460
461
462
463
464
	  }
	}
      }
    }

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


465
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
466
467
468
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
472
473
474
475
476
477
      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, 
478
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
479
480
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
481
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
482
483
484
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
485
486
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
487
488
489
490
491
492
493
494
495
496
		     &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 {
497
498
      MSG("Create direct schur primal solver!\n");

499
500
501
502
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

503
504
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
505
506
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
507

Thomas Witkowski's avatar
Thomas Witkowski committed
508
509
510
511
512
513
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
514

515
516
517
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
518

519
520
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

521
522
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
523
524
525
526
527
528
	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
529
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
530
531
532
533
534
      }

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

535
536
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
		primalDofMap.getLocalDofs())
537
	("Should not happen!\n");
538

539
540
541
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
542
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
543
544
545
546
547
548
549
550
551
552

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
566
567
568
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
569
570
571
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
572
      MatDestroy(&matBPi);
573
574
575
576
577
578
579
580
581
582

      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
583

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


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

594
595
596
597
598
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
599

600
601
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
602

603
604
605
606
607
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
608
609
610
  }


611
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
612
613
614
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

615
616
    // === Create FETI-DP solver object. ===

617
618
619
    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
620
    fetiData.fetiSolver = this;
621
    fetiData.ksp_schur_primal = &ksp_schur_primal;
622

623
    VecCreateMPI(PETSC_COMM_WORLD, 
624
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
625
626
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
627
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
628
		 &(fetiData.tmp_vec_lagrange));
629
    VecCreateMPI(PETSC_COMM_WORLD, 
630
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
631
		 &(fetiData.tmp_vec_primal));
632
633

    MatCreateShell(PETSC_COMM_WORLD,
634
635
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
636
		   &fetiData, &mat_feti);
637
638
639
640
641
642
643
    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);
644
645


646
    // === Create FETI-DP preconditioner object. ===
647

648
649
650
651
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
652

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

689
      VecCreateMPI(PETSC_COMM_WORLD, 
690
691
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
692
		   &(fetiLumpedPreconData.tmp_vec_b));
693
694
695
696
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
697
698
699
700
701
702

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
703
704
      break;
    default:
705
706
      break;
    }
707
708
709
710
711
712
713
  }
  

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

714
715
    // === Destroy FETI-DP solver object. ===

716
717
718
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
719
    fetiData.fetiSolver = NULL;
720
    fetiData.ksp_schur_primal = PETSC_NULL;
721

Thomas Witkowski's avatar
Thomas Witkowski committed
722
723
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
724
    VecDestroy(&fetiData.tmp_vec_primal);
725
726
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
727
728


729
730
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
    // === 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;
757
758
    default:
      break;
759
    }
760
761
762
763
764
765
766
767
768
  }


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

769
    // === Get local part of the solution for B variables. ===
770
771
772
773
774

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


775
776
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
777
    
Thomas Witkowski's avatar
Thomas Witkowski committed
778
779
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);

780
    vector<PetscInt> globalIsIndex, localIsIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
781
782
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());
783
784

    {
Thomas Witkowski's avatar
Thomas Witkowski committed
785
786
787
788
789
790
791
792
      int cnt = 0;
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex>& dofMap = 
	  primalDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin();
	     it != dofMap.end(); ++it) {
	  globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
	  localIsIndex.push_back(cnt++);
793
794
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
795
796

      TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs())("Should not happen!\n");
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
    }
    
    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

Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
    int cnt = 0;
    for (int i = 0; i < vec.getSize(); i++) {
833
834
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));

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

Thomas Witkowski's avatar
Thomas Witkowski committed
841
842
843
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpaces[i]].getMap().begin();
	   it != primalDofMap[feSpaces[i]].getMap().end(); ++it)
	dofVec[it->first] = localSolPrimal[cnt++];
844
845
846
847
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
848
    VecDestroy(&local_sol_primal);
849
850
851
  }


852
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
853
854
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
855
    
856
    // === Create all sets and indices. ===
857
    
Thomas Witkowski's avatar
Thomas Witkowski committed
858
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
859
    
860
861
862
863
    primalDofMap.init(mpiComm, feSpaces, true, true);
    dualDofMap.init(mpiComm, feSpaces, false, false);
    lagrangeMap.init(mpiComm, feSpaces, true, true);
    localDofMap.init(mpiComm, feSpaces, false, false);
864
865
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.init(mpiComm, feSpaces, false, false);
866

867
    updateDofData();
868
    
869
    // === Create matrices for the FETI-DP method. ===
870
    
871
872
    int nRowsRankB = localDofMap.getRankDofs();
    int nRowsOverallB = localDofMap.getOverallDofs();
873
874
    int nRowsRankPrimal = primalDofMap.getRankDofs();
    int nRowsOverallPrimal = primalDofMap.getOverallDofs();
875
    
Thomas Witkowski's avatar
Thomas Witkowski committed
876
877
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
878
    
879
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
880
881
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
882
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
883
    
884
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
885
886
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
887
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
888
    
889
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
890
891
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
892
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
893
894
    
    
895
    // === Create matrices for FETI-DP preconditioner. ===
896
897
898
    
    if (fetiPreconditioner != FETI_NONE) {
      int nRowsDual = dualDofMap.getRankDofs();
899
900

      MatCreateSeqAIJ(PETSC_COMM_SELF,
901
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
902
903
		      &mat_duals_duals);

904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
      if (fetiPreconditioner == FETI_DIRICHLET) {
	int nRowsInterior = interiorDofMap.getRankDofs();

	MatCreateSeqAIJ(PETSC_COMM_SELF,
			nRowsInterior, nRowsInterior, 30, PETSC_NULL,
			&mat_interior_interior);
	
	MatCreateSeqAIJ(PETSC_COMM_SELF,
			nRowsInterior, nRowsDual, 30, PETSC_NULL,
			&mat_interior_duals);
	
	MatCreateSeqAIJ(PETSC_COMM_SELF,
			nRowsDual, nRowsInterior, 30, PETSC_NULL,
			&mat_duals_interior);
      }
919
    }
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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
949
    int nComponents = mat->getSize();
950
951
952
953
954
955
956
957
958
959
960
    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) {
961

Thomas Witkowski's avatar
Thomas Witkowski committed
962
	  bool rowPrimal = isPrimal(feSpaces[i], *cursor);
963
  
964
965
	  cols.clear();
	  colsOther.clear();
966
	  values.clear();	  
967
	  valuesOther.clear();
968
969
970
971
972
973

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

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