Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist ü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. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

PetscSolverFeti.cc 41.3 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 233 234 235 236 237 238 239 240 241 242 243
    if (subDomainSolver == NULL) {
      if (meshLevel == 0) {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);
      } else {
	MeshLevelData& levelData = meshDistributor->getMeshLevelData();

	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, 
			      levelData.getMpiComm(meshLevel - 1),
			      levelData.getMpiComm(meshLevel));
      }
    }
244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265

    primalDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(), 
		      true, true);
    dualDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		    false, false);
    lagrangeMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		     true, true);
    localDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		     false, false);
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.init(mpiComm, feSpaces,  meshDistributor->getFeSpaces(),
			  false, false);
  }


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

267 268 269 270 271 272
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
    interiorDofMap.clear();

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

276 277
      createPrimals(feSpace);      
      createDuals(feSpace);
278
      createLagrange(feSpace);      
279 280
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
281

282 283 284 285 286
    primalDofMap.setDofComm(meshDistributor->getSendDofs(),
			    meshDistributor->getRecvDofs());
    lagrangeMap.setDofComm(meshDistributor->getSendDofs(), 
			   meshDistributor->getRecvDofs());

Thomas Witkowski's avatar
Thomas Witkowski committed
287 288 289 290
    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();
291
    if (fetiPreconditioner != FETI_NONE)
292
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
293 294 295

    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
296
      MSG("FETI-DP data for %d-ith FE space:\n", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
297

298
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
299 300
	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
      
301
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
302
	  dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
303

304
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
305
	  lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
306 307 308 309

      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
310
    }
311 312 313


    // If multi level test, inform sub domain solver about coarse space.
314
    subDomainSolver->setDofMapping(&primalDofMap, &localDofMap);
315 316 317
  }


318
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
319
  {
320
    FUNCNAME("PetscSolverFeti::createPrimals()");  
321

322 323 324
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

325 326 327
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
328
    DofContainerSet& vertices = 
329
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
330 331 332
    for (DofContainerSet::iterator it = vertices.begin(); 
	 it != vertices.end(); ++it)
      primals.insert(**it);
333

334 335 336 337

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

338
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
339 340
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
341 342
      else
  	primalDofMap[feSpace].insert(*it);
343 344 345
  }


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

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

    DofContainer allBoundaryDofs;
353
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
354 355 356 357 358 359 360 361 362 363 364 365

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

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

369
    boundaryDofRanks[feSpace].clear();
370

371
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), meshLevel, feSpace); 
372 373
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
374
	if (!isPrimal(feSpace, it.getDofIndex())) {
375 376
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
377 378
	}
      }
379
	
380 381 382 383

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

384
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
385

386
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), meshLevel, feSpace);
387 388
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
389
	if (!isPrimal(feSpace, it.getDofIndex()))
390
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
391 392 393

    stdMpi.updateSendDataSize();

394
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), meshLevel, feSpace); 
395
	 !it.end(); it.nextRank()) {
396
      bool recvFromRank = false;
397
      for (; !it.endDofIter(); it.nextDof()) {
398
	if (!isPrimal(feSpace, it.getDofIndex())) {
399 400 401
	  recvFromRank = true;
	  break;
	}
402
      }
403 404

      if (recvFromRank)
405
	stdMpi.recv(it.getRank());
406
    }
407

408 409
    stdMpi.startCommunication();

410
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), meshLevel, feSpace); 
411
	 !it.end(); it.nextRank()) {
412
      int i = 0;
413
      for (; !it.endDofIter(); it.nextDof())
414
	if (!isPrimal(feSpace, it.getDofIndex()))
415
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
416
	    stdMpi.getRecvData(it.getRank())[i++];
417 418
    }

419 420 421
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

422
    int nRankLagrange = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
423 424
    map<DegreeOfFreedom, MultiIndex>& dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
425
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
426
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
427
	int degree = boundaryDofRanks[feSpace][it->first].size();
428
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
429 430
      } else {
	lagrangeMap[feSpace].insert(it->first);
431 432
      }
    }
433
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
434 435 436
  }


437
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
438
  {
439
    FUNCNAME("PetscSolverFeti::createIndexB()");
440

441
    DOFAdmin* admin = feSpace->getAdmin();
442 443 444 445

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

447
    int nLocalInterior = 0;    
448
    for (int i = 0; i < admin->getUsedSize(); i++) {
449
      if (admin->isDofFree(i) == false && 
450
	  isPrimal(feSpace, i) == false &&
451
	  dualDofMap[feSpace].isSet(i) == false) {
452 453
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);

454
	if (fetiPreconditioner != FETI_NONE)
455 456 457
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);

	nLocalInterior++;
458
      }
459
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
460
    
461 462
    // === And finally, add the global indicies of all dual nodes. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
463
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualDofMap[feSpace].getMap().begin();
464
	 it != dualDofMap[feSpace].getMap().end(); ++it)
465
      localDofMap[feSpace].insertRankDof(it->first);
466 467 468
  }


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

473 474
    // === Create distributed matrix for Lagrange constraints. ===

475
    MatCreateMPIAIJ(mpiComm,
476 477
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
478 479 480
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

481 482 483 484 485 486 487
    // === 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
488 489 490
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
      map<DegreeOfFreedom, MultiIndex> &dualMap = 
	dualDofMap[feSpaces[k]].getMap();
491

Thomas Witkowski's avatar
Thomas Witkowski committed
492 493
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	   it != dualMap.end(); ++it) {
494
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
495 496 497 498 499 500
	  ("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.
501 502
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
503 504 505 506 507 508
	// 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) {
509
	      int colIndex = localDofMap.getMatIndex(k, it->first);
510
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
511
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
512
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
513
	    index++;	      
514 515 516 517 518 519 520 521 522 523
	  }
	}
      }
    }

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


524
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
525 526 527
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

531
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
532
      
533
      VecCreateMPI(mpiComm, 
534
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
535
		   &(schurPrimalData.tmp_vec_b));
536
      VecCreateMPI(mpiComm, 
537
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
538 539
		   &(schurPrimalData.tmp_vec_primal));

540
      MatCreateShell(mpiComm,
541 542
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
543 544 545 546 547
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
548
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
549
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
550 551
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
552 553
      KSPSetFromOptions(ksp_schur_primal);
    } else {
554 555
      MSG("Create direct schur primal solver!\n");

556 557
      double wtime = MPI::Wtime();

558 559
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
560 561
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
562

Thomas Witkowski's avatar
Thomas Witkowski committed
563
      Mat matBPi;
564
      MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
565 566 567 568
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
569

570 571 572
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
573

574 575
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

576 577
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
578
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
579 580 581 582 583

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

584
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
585 586
      }

587 588
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
		primalDofMap.getLocalDofs())
589
	("Should not happen!\n");
590

591 592 593
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
594
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
595 596 597 598 599 600 601 602

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

603
	subDomainSolver->solve(tmpVec, tmpVec);
604

Thomas Witkowski's avatar
Thomas Witkowski committed
605
	PetscScalar *tmpValues;
606
	VecGetArray(tmpVec, &tmpValues);
607
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
608
	  MatSetValue(matBPi, 
609
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
610 611 612
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
613 614 615 616 617
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
618 619
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
620 621

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
622
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
623
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
624 625

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
626
      MatDestroy(&matBPi);
627 628

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

632
      KSPCreate(mpiComm, &ksp_schur_primal);
633
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
634
		      SAME_NONZERO_PATTERN);
635 636 637 638 639 640
      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);
641
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
642

643 644
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
645
    }
646 647 648 649 650 651 652
  }


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

653
    if (schurPrimalSolver == 0) {
654
      schurPrimalData.subSolver = NULL;
655

656 657 658
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
659 660 661
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
662 663 664
  }


665
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
666 667 668
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

669 670
    // === Create FETI-DP solver object. ===

671
    fetiData.mat_lagrange = &mat_lagrange;
672
    fetiData.subSolver = subDomainSolver;
673
    fetiData.ksp_schur_primal = &ksp_schur_primal;
674

675
    VecCreateMPI(mpiComm, 
676
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
677
		 &(fetiData.tmp_vec_b));
678
    VecCreateMPI(mpiComm,
679
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
680
		 &(fetiData.tmp_vec_lagrange));
681
    VecCreateMPI(mpiComm, 
682
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
683
		 &(fetiData.tmp_vec_primal));
684

685
    MatCreateShell(mpiComm,
686 687
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
688
		   &fetiData, &mat_feti);
689 690 691
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


692
    KSPCreate(mpiComm, &ksp_feti);
693
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
694 695 696
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
697
    KSPSetFromOptions(ksp_feti);
698 699


700
    // === Create FETI-DP preconditioner object. ===
701

702 703 704 705
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
706

707 708 709
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
710 711
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
712 713 714 715 716 717
      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);
718 719 720 721 722 723 724 725 726
      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;
      
727
      VecCreateMPI(mpiComm, 
728
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
729
		   &(fetiDirichletPreconData.tmp_vec_b));      
730 731 732 733 734 735
      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));
736 737 738 739 740 741 742 743 744 745 746 747 748

      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex> &dualMap = 
	  dualDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	     it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

749 750 751 752 753 754 755 756 757 758 759
      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
760 761 762 763 764 765 766 767 768 769 770 771
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex> &dualMap = 
	  dualDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	     it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

772
      VecCreateMPI(mpiComm, 
773 774
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
775
		   &(fetiLumpedPreconData.tmp_vec_b));
776 777 778 779
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
780 781 782 783 784 785

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
786 787
      break;
    default:
788 789
      break;
    }
790 791 792 793 794 795 796
  }
  

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

797 798
    // === Destroy FETI-DP solver object. ===

799
    fetiData.mat_lagrange = PETSC_NULL;
800
    fetiData.subSolver = NULL;
801
    fetiData.ksp_schur_primal = PETSC_NULL;
802

Thomas Witkowski's avatar
Thomas Witkowski committed
803 804
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
805
    VecDestroy(&fetiData.tmp_vec_primal);
806 807
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
808 809


810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837
    // === 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;