PetscSolverFeti.cc 42.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// 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.


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

namespace AMDiS {

  using namespace std;

25 26 27 28 29 30 31
  // 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);
32
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
33

34 35 36 37
    MatMult(data->subSolver->getMatIntCoarse(), x, data->tmp_vec_b);
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarseCoarse(), x, y);
38 39 40 41 42 43 44 45 46
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
47 48
    //    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)
49 50 51

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

Thomas Witkowski's avatar
Thomas Witkowski committed
54
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
55
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
57

58
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
59
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
60
   MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
61
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
63

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

    return 0;
  }


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

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


82
    // === Restriction of the B nodes to the boundary nodes. ===
83

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

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

93 94 95
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
96 97

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


101
    // === K_DD - K_DI inv(K_II) K_ID ===
102

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

105
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
106
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
107 108 109 110 111 112 113 114 115 116
    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);

117 118 119
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141

    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);
142 143


144
    // === Restriction of the B nodes to the boundary nodes. ===
145

146 147 148 149
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
150

151
    PetscScalar *local_b, *local_duals;
152
    VecGetArray(data->tmp_vec_b, &local_b);
153
    VecGetArray(data->tmp_vec_duals0, &local_duals);
154

Thomas Witkowski's avatar
Thomas Witkowski committed
155 156 157
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
158 159

    VecRestoreArray(data->tmp_vec_b, &local_b);
160 161 162 163 164 165 166
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

167

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

170 171 172
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

Thomas Witkowski's avatar
Thomas Witkowski committed
173 174 175
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
176

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

180 181

    // Multiply with scaled Lagrange constraint matrix.
182 183 184 185 186 187
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


188 189
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
190
      schurPrimalSolver(0),
191
      multiLevelTest(false),
192 193
      subDomainSolver(NULL),
      meshLevel(0)
194 195 196 197 198 199
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
200
      MSG("Create FETI-DP solver with no preconditioner!\n");
201 202
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
203
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
204 205
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
206
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
207 208
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
209 210
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
211
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
212 213 214 215 216

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

    Parameters::get("parallel->multi level test", multiLevelTest);
219 220
    if (multiLevelTest)
      meshLevel = 1;
221 222 223
  }


224
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
225
  {
226 227
    FUNCNAME("PetscSolverFeti::initialize()");

228 229 230
    TEST_EXIT_DBG(meshLevel + 1 == meshDistributor->getMeshLevelData().getLevelNumber())
      ("Mesh hierarchy does not contain %d levels!\n", meshLevel + 1);

231 232
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

233 234 235 236 237 238 239 240 241
    if (subDomainSolver == NULL) {
      if (meshLevel == 0) {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);
      } else {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, 
			      levelData.getMpiComm(meshLevel - 1),
			      levelData.getMpiComm(meshLevel));
242
	subDomainSolver->setLevel(meshLevel);
243 244
      }
    }
245

246 247
    primalDofMap.init(levelData,
		      feSpaces, meshDistributor->getFeSpaces(), 
248
		      true, true);
249 250
    primalDofMap.setComputeMatIndex(true);
    
251 252
    dualDofMap.init(levelData, 
		    feSpaces, meshDistributor->getFeSpaces(),
253
		    false, false);
254 255
    dualDofMap.setComputeMatIndex(true);

256 257
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
258
		     true, true);
259 260
    lagrangeMap.setComputeMatIndex(true);
 
261 262
    localDofMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
263
		     false, false);
264 265 266
    localDofMap.setComputeMatIndex(true);

    if (fetiPreconditioner == FETI_DIRICHLET) {
267 268
      interiorDofMap.init(levelData, 
			  feSpaces,  meshDistributor->getFeSpaces(),
269
			  false, false);
270 271
      interiorDofMap.setComputeMatIndex(true);
    }
272 273 274 275 276 277 278 279 280 281
  }


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

    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
282

283 284 285 286
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
287 288 289 290 291 292 293
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.clear();

    primalDofMap.setDofComm(meshDistributor->getSendDofs(),
			    meshDistributor->getRecvDofs());
    lagrangeMap.setDofComm(meshDistributor->getSendDofs(), 
			   meshDistributor->getRecvDofs());
294

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

298 299
      createPrimals(feSpace);  

300
      createDuals(feSpace);      
301

302 303
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
304 305 306 307

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
308
    if (fetiPreconditioner != FETI_NONE)
309
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
310

311 312 313 314 315 316 317 318 319
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);

      MSG("TEST DUALS %d %d %d\n",
	  dualDofMap[feSpace].nLocalDofs[meshLevel], 
	  dualDofMap[feSpace].nRankDofs[meshLevel], 
	  dualDofMap[feSpace].nOverallDofs[meshLevel]);
    }	

320 321 322 323 324 325 326
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
    lagrangeMap.update();


327
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
Thomas Witkowski's avatar
Thomas Witkowski committed
328 329
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
330
      MSG("FETI-DP data for %d-ith FE space:\n", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
331

332
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
333 334
	  primalDofMap[feSpace].nRankDofs[meshLevel], 
	  primalDofMap[feSpace].nOverallDofs[meshLevel]);
Thomas Witkowski's avatar
Thomas Witkowski committed
335
      
336
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
337 338
	  dualDofMap[feSpace].nRankDofs[meshLevel], 
	  dualDofMap[feSpace].nOverallDofs[meshLevel]);
Thomas Witkowski's avatar
Thomas Witkowski committed
339

340
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
341 342
	  lagrangeMap[feSpace].nRankDofs[meshLevel], 
	  lagrangeMap[feSpace].nOverallDofs[meshLevel]);
343

344
      TEST_EXIT_DBG(localDofMap[feSpace].size(meshLevel) + primalDofMap[feSpace].size(meshLevel) == 
345 346
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
347
    }
348 349 350


    // If multi level test, inform sub domain solver about coarse space.
351
    subDomainSolver->setDofMapping(&primalDofMap, &localDofMap);
352 353 354
  }


355
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
356
  {
357
    FUNCNAME("PetscSolverFeti::createPrimals()");  
358

359 360 361
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

362 363 364
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
365
    DofContainerSet& vertices = 
366
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
367 368 369
    for (DofContainerSet::iterator it = vertices.begin(); 
	 it != vertices.end(); ++it)
      primals.insert(**it);
370

371 372 373 374

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

375
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
376 377
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it, meshLevel))
	primalDofMap[feSpace].insertRankDof(meshLevel, *it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
378
      else
379
  	primalDofMap[feSpace].insertNonRankDof(meshLevel, *it);
380 381 382
  }


383
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
384 385
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
386

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

    DofContainer allBoundaryDofs;
390
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
391 392 393 394

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
395
	dualDofMap[feSpace].insertRankDof(meshLevel, **it);
396 397 398 399 400 401 402
  }

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

403 404 405 406 407
    boundaryDofRanks[feSpace].clear();

    if (dualDofMap[feSpace].nLocalDofs[meshLevel] == 0)
      return;

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

411
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), meshLevel, feSpace); 
412 413
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
414
	if (!isPrimal(feSpace, it.getDofIndex())) {
415 416
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
417 418
	}
      }
419
	
420 421 422 423

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

424
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
425

426
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), meshLevel, feSpace);
427 428
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
429
	if (!isPrimal(feSpace, it.getDofIndex()))
430
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
431 432 433

    stdMpi.updateSendDataSize();

434
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), meshLevel, feSpace); 
435
	 !it.end(); it.nextRank()) {
436
      bool recvFromRank = false;
437
      for (; !it.endDofIter(); it.nextDof()) {
438
	if (!isPrimal(feSpace, it.getDofIndex())) {
439 440 441
	  recvFromRank = true;
	  break;
	}
442
      }
443 444

      if (recvFromRank)
445
	stdMpi.recv(it.getRank());
446
    }
447

448 449
    stdMpi.startCommunication();

450
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), meshLevel, feSpace); 
451
	 !it.end(); it.nextRank()) {
452
      int i = 0;
453
      for (; !it.endDofIter(); it.nextDof())
454
	if (!isPrimal(feSpace, it.getDofIndex()))
455
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
456
	    stdMpi.getRecvData(it.getRank())[i++];
457 458
    }

459 460 461
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

462
    int nRankLagrange = 0;
463
    DofMap& dualMap = dualDofMap[feSpace].getMap(meshLevel);
464
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
465

466 467
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first, meshLevel)) {
	lagrangeMap[feSpace].insertRankDof(meshLevel, it->first, nRankLagrange);
468
	int degree = boundaryDofRanks[feSpace][it->first].size();
469
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
470
      } else {
471
	lagrangeMap[feSpace].insertNonRankDof(meshLevel, it->first);
472 473
      }
    }
474
    lagrangeMap[feSpace].nRankDofs[meshLevel] = nRankLagrange;
475 476 477
  }


478
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
479
  {
480
    FUNCNAME("PetscSolverFeti::createIndexB()");
481

482
    DOFAdmin* admin = feSpace->getAdmin();
483 484 485 486

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

488
    int nLocalInterior = 0;    
489
    for (int i = 0; i < admin->getUsedSize(); i++) {
490
      if (admin->isDofFree(i) == false && 
491
	  isPrimal(feSpace, i) == false &&
492 493
	  isDual(feSpace, i) == false) {
	localDofMap[feSpace].insertRankDof(meshLevel, i, nLocalInterior);
494

495
	if (fetiPreconditioner != FETI_NONE)
496
	  interiorDofMap[feSpace].insertRankDof(meshLevel, i, nLocalInterior);
497 498

	nLocalInterior++;
499
      }
500
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
501
    
502 503
    // === And finally, add the global indicies of all dual nodes. ===

504 505 506
    for (DofMap::iterator it = dualDofMap[feSpace].getMap(meshLevel).begin();
	 it != dualDofMap[feSpace].getMap(meshLevel).end(); ++it)
      localDofMap[feSpace].insertRankDof(meshLevel, it->first);
507 508 509
  }


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

514 515
    // === Create distributed matrix for Lagrange constraints. ===

516
    MatCreateMPIAIJ(mpiComm,
517 518 519 520
		    lagrangeMap.getRankDofs(meshLevel), 
		    localDofMap.getRankDofs(meshLevel),
		    lagrangeMap.getOverallDofs(meshLevel), 
		    localDofMap.getOverallDofs(meshLevel),
521 522 523
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

524 525 526 527 528 529 530
    // === 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
531
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
532
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap(meshLevel);
533

534
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
535
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
536 537 538
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
539
	int index = lagrangeMap.getMatIndex(meshLevel, k, it->first);
Thomas Witkowski's avatar
Thomas Witkowski committed
540 541
	
	// Copy set of all ranks that contain this dual node.
542 543
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
544 545 546 547 548 549
	// 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) {
550
	      int colIndex = localDofMap.getMatIndex(meshLevel, k, it->first);
551
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
552
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
553
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
554
	    index++;	      
555 556 557 558 559 560 561 562 563 564
	  }
	}
      }
    }

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


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

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

572
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
573
      
574
      VecCreateMPI(mpiComm, 
575 576
		   localDofMap.getRankDofs(meshLevel), 
		   localDofMap.getOverallDofs(meshLevel),
Thomas Witkowski's avatar
Thomas Witkowski committed
577
		   &(schurPrimalData.tmp_vec_b));
578
      VecCreateMPI(mpiComm, 
579 580
		   primalDofMap.getRankDofs(meshLevel), 
		   primalDofMap.getOverallDofs(meshLevel),
Thomas Witkowski's avatar
Thomas Witkowski committed
581 582
		   &(schurPrimalData.tmp_vec_primal));

583
      MatCreateShell(mpiComm,
584 585 586 587
		     primalDofMap.getRankDofs(meshLevel), 
		     primalDofMap.getRankDofs(meshLevel), 
		     primalDofMap.getOverallDofs(meshLevel), 
		     primalDofMap.getOverallDofs(meshLevel),
Thomas Witkowski's avatar
Thomas Witkowski committed
588 589 590 591 592
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
593
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
594
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
595 596
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
597 598
      KSPSetFromOptions(ksp_schur_primal);
    } else {
599 600
      MSG("Create direct schur primal solver!\n");

601 602
      double wtime = MPI::Wtime();

603 604 605 606
      int nRowsRankPrimal = primalDofMap.getRankDofs(meshLevel);
      int nRowsOverallPrimal = primalDofMap.getOverallDofs(meshLevel);
      int nRowsRankB = localDofMap.getRankDofs(meshLevel);
      int nRowsOverallB = localDofMap.getOverallDofs(meshLevel);
607

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

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

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

621
      for (int i = 0; i < nRowsRankB; i++) {
622
	PetscInt row = localDofMap.getStartDofs(meshLevel) + i;
623
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
624 625 626 627 628

	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]));

629
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
630 631
      }

632
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
633
		primalDofMap.getLocalDofs(meshLevel))
634
	("Should not happen!\n");
635

636 637 638
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
639
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
640 641 642 643 644 645 646 647

 	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);

648
	subDomainSolver->solve(tmpVec, tmpVec);
649

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
663 664
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
665 666

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
667
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
668
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
669 670

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
671
      MatDestroy(&matBPi);
672 673

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

677
      KSPCreate(mpiComm, &ksp_schur_primal);
678
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
679
		      SAME_NONZERO_PATTERN);
680 681 682 683 684 685
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPPREONLY);
      PC pc_schur_primal;      
      KSPGetPC(ksp_schur_primal, &pc_schur_primal);
      PCSetType(pc_schur_primal, PCLU);
      PCFactorSetMatSolverPackage(pc_schur_primal, MATSOLVERMUMPS);
686
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
687

688 689
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
690
    }
691 692 693 694 695 696 697
  }


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

698
    if (schurPrimalSolver == 0) {
699
      schurPrimalData.subSolver = NULL;
700

701 702 703
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
704 705 706
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
707 708 709
  }


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

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

716
    fetiData.mat_lagrange = &mat_lagrange;
717
    fetiData.subSolver = subDomainSolver;
718
    fetiData.ksp_schur_primal = &ksp_schur_primal;
719

720
    VecCreateMPI(mpiComm, 
721 722
		 localDofMap.getRankDofs(meshLevel), 
		 localDofMap.getOverallDofs(meshLevel),
Thomas Witkowski's avatar
Thomas Witkowski committed
723
		 &(fetiData.tmp_vec_b));
724
    VecCreateMPI(mpiComm,
725 726
		 lagrangeMap.getRankDofs(meshLevel), 
		 lagrangeMap.getOverallDofs(meshLevel),
Thomas Witkowski's avatar
Thomas Witkowski committed
727
		 &(fetiData.tmp_vec_lagrange));
728
    VecCreateMPI(mpiComm, 
729 730
		 primalDofMap.getRankDofs(meshLevel), 
		 primalDofMap.getOverallDofs(meshLevel),
731
		 &(fetiData.tmp_vec_primal));
732

733
    MatCreateShell(mpiComm,
734 735 736 737
		   lagrangeMap.getRankDofs(meshLevel), 
		   lagrangeMap.getRankDofs(meshLevel),
		   lagrangeMap.getOverallDofs(meshLevel), 
		   lagrangeMap.getOverallDofs(meshLevel),
738
		   &fetiData, &mat_feti);
739 740 741
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


742
    KSPCreate(mpiComm, &ksp_feti);
743
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
744 745 746
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
747
    KSPSetFromOptions(ksp_feti);
748 749


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

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

757 758 759
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
760 761
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
762 763 764 765 766 767
      KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
      KSPSetType(ksp_interior, KSPPREONLY);
      PC pc_interior;
      KSPGetPC(ksp_interior, &pc_interior);
      PCSetType(pc_interior, PCLU);
      PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK);
768 769 770 771 772 773 774 775 776
      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;
      
777
      VecCreateMPI(mpiComm, 
778 779
		   localDofMap.getRankDofs(meshLevel),
		   localDofMap.getOverallDofs(meshLevel),
780
		   &(fetiDirichletPreconData.tmp_vec_b));      
781 782 783 784 785 786
      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));
787 788

      for (unsigned int i = 0; i < feSpaces.size(); i++) {
789
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(meshLevel);
790
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
791
	  DegreeOfFreedom d = it->first;
792 793
	  int matIndexLocal = localDofMap.getLocalMatIndex(meshLevel, i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(meshLevel, i, d);
794 795 796 797
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

798 799 800 801 802 803 804 805 806 807 808
      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;

Thomas Witkowski's avatar
Thomas Witkowski committed
809
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
810
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(meshLevel);
811
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
812
	  DegreeOfFreedom d = it->first;
813 814
	  int matIndexLocal = localDofMap.