Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

PetscSolverFeti.cc 54.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


#include "parallel/PetscSolverFeti.h"
14
#include "parallel/PetscSolverFetiStructs.h"
15 16
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
17
#include "io/VtkWriter.h"
18 19 20 21 22

namespace AMDiS {

  using namespace std;

23 24 25 26 27 28 29
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
30
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
31 32

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
33
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
34 35 36 37 38 39 40 41 42 43 44
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(*(data->mat_primal_primal), x, y);
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
45 46
    //    F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
47 48 49

    void *ctx;
    MatShellGetContext(mat, &ctx);
50
    FetiData* data = static_cast<FetiData*>(ctx);
51

Thomas Witkowski's avatar
Thomas Witkowski committed
52 53 54
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
55

Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
57
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
58 59 60
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
61

Thomas Witkowski's avatar
Thomas Witkowski committed
62
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
63 64 65 66 67

    return 0;
  }


68
  // y = PC * x
69
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
70
  {
71
    // Get data for the preconditioner
72 73
    void *ctx;
    PCShellGetContext(pc, &ctx);
74
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
75

76
    // Multiply with scaled Lagrange constraint matrix.
77 78 79
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


80
    // === Restriction of the B nodes to the boundary nodes. ===
81

82 83 84 85
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
86

87 88 89 90 91 92
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];
93 94

    VecRestoreArray(data->tmp_vec_b, &local_b);
95
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
96 97


98
    // === K_DD - K_DI inv(K_II) K_ID ===
99

100
    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
101

102
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
103
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

    VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);


    // === Prolongation from local dual nodes to B nodes.

    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];

    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
138 139


140
    // === Restriction of the B nodes to the boundary nodes. ===
141

142 143 144 145
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
146

147
    PetscScalar *local_b, *local_duals;
148
    VecGetArray(data->tmp_vec_b, &local_b);
149
    VecGetArray(data->tmp_vec_duals0, &local_duals);
150

151 152
    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];
153 154

    VecRestoreArray(data->tmp_vec_b, &local_b);
155 156 157 158 159 160 161
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);

162

163
    // === Prolongation from local dual nodes to B nodes.
164

165 166 167 168 169
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];
170

171 172
    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
173

174 175

    // Multiply with scaled Lagrange constraint matrix.
176 177 178 179 180 181
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


182 183
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
Thomas Witkowski's avatar
Thomas Witkowski committed
184 185
      nComponents(-1),
      schurPrimalSolver(0)
186 187 188 189 190 191
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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


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

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
218
   
219 220
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
221

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

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


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

235 236 237
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

238 239 240
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
241
    DofContainerSet& vertices = 
242
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
243 244 245 246
    TEST_EXIT_DBG(vertices.size())("No primal vertices on this rank!\n");
    for (DofContainerSet::iterator it = vertices.begin(); 
	 it != vertices.end(); ++it)
      primals.insert(**it);
247

248 249 250 251

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

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

256
    primalDofMap[feSpace].update();
257

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

262 263 264 265

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

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

278 279
    stdMpi.updateSendDataSize();

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

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

295 296
    stdMpi.startCommunication();

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

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


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

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

    boundaryDofRanks.clear();

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

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

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

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

    stdMpi.updateSendDataSize();

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

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

362 363
    stdMpi.startCommunication();

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

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

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

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

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

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

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

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

397 398

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


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

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

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

428 429
    stdMpi.updateSendDataSize();

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

    stdMpi.startCommunication();

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


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

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

    // === To ensure that all interior node on each rank are listen first in ===
    // === the global index of all B nodes, insert all interior nodes first, ===
    // === without defining a correct index.                                 ===
468 469 470

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

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

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

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

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

    localDofMap[feSpace].update(false);

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


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

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

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

516 517 518 519 520 521 522
    // === Create for all duals the corresponding Lagrange constraints. On ===
    // === each rank we traverse all pairs (n, m) of ranks, with n < m,    ===
    // === that contain this node. If the current rank number is r, and    ===
    // === n == r, the rank sets 1.0 for the corresponding constraint, if  ===
    // === m == r, than the rank sets -1.0 for the corresponding           ===
    // === constraint.                                                     ===

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

      DofMapping &dualMap = dualDofMap[feSpace].getMap();
      for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
	TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	  ("Should not happen!\n");

	// Global index of the first Lagrange constriant for this node.
	int constraintMatIndex = lagrangeMap.mapGlobal(it->first, feIdx);
	// Copy set of all ranks that contain this dual node.
	vector<int> W(boundaryDofRanks[it->first].begin(), 
		      boundaryDofRanks[it->first].end());
	// Number of ranks that contain this dual node.
	int degree = W.size();
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
#if 1
	      // Get global index of the dual variable
	      int colIndex = localDofMap.mapGlobal(it->first, feIdx);
#else
546
	      int colIndex = it->second * nComponents + k;
547
#endif
548
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
549 550
	      MatSetValue(mat_lagrange, constraintMatIndex, colIndex, value, 
			  INSERT_VALUES);	      
551
	    }
552 553
	    
	    constraintMatIndex++;
554 555 556 557 558 559 560 561 562 563
	  }
	}
      }
    }

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


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
573 574 575 576 577 578
      schurPrimalData.mat_primal_primal = &mat_primal_primal;
      schurPrimalData.mat_primal_b = &mat_primal_b;
      schurPrimalData.mat_b_primal = &mat_b_primal;
      schurPrimalData.fetiSolver = this;
      
      VecCreateMPI(PETSC_COMM_WORLD, 
579
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
580 581
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
582
		   primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
583 584 585
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
586 587
		     primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nRankDofs * nComponents,
		     primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
588 589 590 591 592 593 594 595 596 597
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
    } else {
598 599
      MSG("Create direct schur primal solver!\n");

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

      double wtime = MPI::Wtime();

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

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

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

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

622 623
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
624 625 626 627 628 629
	MatGetRow(mat_b_primal, row, &nCols, &cols, &values);

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

Thomas Witkowski's avatar
Thomas Witkowski committed
630
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
631 632 633 634 635 636
      }

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

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

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);

       	KSPSolve(ksp_b, tmpVec, tmpVec);

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

	VecDestroy(&tmpVec);
      }

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

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

      MatInfo minfo;
      MatGetInfo(mat_primal_primal, MAT_GLOBAL_SUM, &minfo);
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_primal_primal, 
		      mat_primal_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
682

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


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

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

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

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


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

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

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

718 719 720
    fetiData.mat_primal_b = &mat_primal_b;
    fetiData.mat_b_primal = &mat_b_primal;
    fetiData.mat_lagrange = &mat_lagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
721
    fetiData.fetiSolver = this;
722
    fetiData.ksp_schur_primal = &ksp_schur_primal;
723

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

    MatCreateShell(PETSC_COMM_WORLD,
737 738 739 740
		   lagrangeMap[feSpace].nRankDofs * nComponents, 
		   lagrangeMap[feSpace].nRankDofs * nComponents,
		   lagrangeMap[feSpace].nOverallDofs * nComponents, 
		   lagrangeMap[feSpace].nOverallDofs * nComponents,
741
		   &fetiData, &mat_feti);
742 743 744 745 746 747 748
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


    KSPCreate(PETSC_COMM_WORLD, &ksp_feti);
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_feti, "solver_feti_");
    KSPSetFromOptions(ksp_feti);
749 750


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

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

758 759 760
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
761 762
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
763 764 765 766 767 768 769 770 771 772
      KSPSetOptionsPrefix(ksp_interior, "solver_interior_");
      KSPSetFromOptions(ksp_interior);
            
      fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
      fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
      fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
      fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
      fetiDirichletPreconData.ksp_interior = &ksp_interior;
      
773
      VecCreateMPI(PETSC_COMM_WORLD, 
774
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
775
		   &(fetiDirichletPreconData.tmp_vec_b));      
776 777 778 779 780 781
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_interior));
782 783 784 785 786 787 788 789 790 791 792 793
      
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
      
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;

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

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

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

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

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

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


834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861
    // === Destroy FETI-DP preconditioner object. ===

    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPDestroy(&ksp_interior);

      fetiDirichletPreconData.mat_lagrange_scaled = NULL;
      fetiDirichletPreconData.mat_interior_interior = NULL;
      fetiDirichletPreconData.mat_duals_duals = NULL;
      fetiDirichletPreconData.mat_interior_duals = NULL;
      fetiDirichletPreconData.mat_duals_interior = NULL;
      fetiDirichletPreconData.ksp_interior = NULL;
      
      VecDestroy(&fetiDirichletPreconData.tmp_vec_b);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals0);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals1);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_interior);
      MatDestroy(&mat_lagrange_scaled);
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = NULL;
      fetiLumpedPreconData.mat_duals_duals = NULL;

      VecDestroy(&fetiLumpedPreconData.tmp_vec_b);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1);
      break;
862 863
    default:
      break;
864
    }
865 866 867 868 869 870 871 872 873
  }


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

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