PetscSolverFeti.cc 56.2 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"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscSolverFeti.h"
16
#include "parallel/PetscSolverFetiStructs.h"
17 18
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
19
#include "parallel/PetscSolverGlobalMatrix.h"
20
#include "io/VtkWriter.h"
21 22 23 24 25

namespace AMDiS {

  using namespace std;

26 27 28 29 30 31
  double FetiTimings::fetiSolve = 0.0;
  double FetiTimings::fetiSolve01 = 0.0;
  double FetiTimings::fetiSolve02 = 0.0;
  double FetiTimings::fetiPreconditioner = 0.0;


32 33 34 35 36 37 38
  // 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);
39
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
40

41
    MatMult(data->subSolver->getMatInteriorCoarse(), x, data->tmp_vec_b);
42
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
43 44 45
    MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b, 
	    data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarse(), x, y);
46 47 48 49 50 51 52 53 54
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
55 56
    FUNCNAME("petscMultMatFeti()");

57 58
    //    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)
59

60 61
    double wtime = MPI::Wtime();

62 63
    void *ctx;
    MatShellGetContext(mat, &ctx);
64
    FetiData* data = static_cast<FetiData*>(ctx);
65

Thomas Witkowski's avatar
Thomas Witkowski committed
66
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
67 68

    double wtime01 = MPI::Wtime();
69
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
70 71 72

    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
73
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
74

75 76
    MatMult(data->subSolver->getMatCoarseInterior(), 
	    data->tmp_vec_b, data->tmp_vec_primal);
77 78

    wtime01 = MPI::Wtime();
79
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
80 81
    FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);

82 83
    MatMult(data->subSolver->getMatInteriorCoarse(), 
	    data->tmp_vec_primal, data->tmp_vec_b);
84 85

    wtime01 = MPI::Wtime();
86
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
87 88
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
89
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
90

Thomas Witkowski's avatar
Thomas Witkowski committed
91
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
92

93 94
    FetiTimings::fetiSolve += (MPI::Wtime() - wtime);

95 96 97 98
    return 0;
  }


99
  // y = PC * x
100
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
101
  {
102 103
    double wtime = MPI::Wtime();

104
    // Get data for the preconditioner
105 106
    void *ctx;
    PCShellGetContext(pc, &ctx);
107
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
108

109
    // Multiply with scaled Lagrange constraint matrix.
110 111 112
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


113
    // === Restriction of the B nodes to the boundary nodes. ===
114

115 116 117 118
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
119

120 121 122 123
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

124 125 126
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
127 128

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


132
    // === K_DD - K_DI inv(K_II) K_ID ===
133

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

136
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
137
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
138 139 140 141 142 143 144 145 146 147
    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);

148 149 150
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
151 152 153 154 155 156 157 158

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

159 160
    FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);

161 162 163 164 165 166 167 168 169 170 171 172 173 174
    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);
175 176


177
    // === Restriction of the B nodes to the boundary nodes. ===
178

179 180 181 182
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
183

184
    PetscScalar *local_b, *local_duals;
185
    VecGetArray(data->tmp_vec_b, &local_b);
186
    VecGetArray(data->tmp_vec_duals0, &local_duals);
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188 189 190
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
191 192

    VecRestoreArray(data->tmp_vec_b, &local_b);
193 194 195 196 197 198 199
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

200

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

203 204 205
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

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

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

213 214

    // Multiply with scaled Lagrange constraint matrix.
215 216 217 218 219 220
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


221 222
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
223
      schurPrimalSolver(0),
224
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
225
      subdomain(NULL),
226 227
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
228
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
229 230
      printTimings(false),
      enableStokesMode(false)
231 232 233 234 235 236
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
237
      MSG("Create FETI-DP solver with no preconditioner!\n");
238 239
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
240
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
241 242
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
243
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
244 245
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
246 247
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
248
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
249 250 251 252 253

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

    Parameters::get("parallel->multi level test", multiLevelTest);
256 257
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
258 259

    Parameters::get("parallel->print timings", printTimings);
260 261 262
  }


263
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
264
  {
265 266
    FUNCNAME("PetscSolverFeti::initialize()");

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

270
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
271
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
272

Thomas Witkowski's avatar
Thomas Witkowski committed
273 274 275 276 277 278 279 280 281

    enableStokesMode = false;
    Parameters::get("parallel->stokes mode", enableStokesMode);
    if (enableStokesMode) {
      MSG("========= !!!! SUPER WARNING !!! ========\n");
      MSG("Use make use of stokes mode which is work in progress! We guarantee for nothing and even less than this!\n");
      fullInterface = feSpaces[2];
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
282 283
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
284

285
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
286
	subdomain->setMeshDistributor(meshDistributor, 
287
				      mpiCommGlobal, mpiCommLocal);
288
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
289
	subdomain->setMeshDistributor(meshDistributor, 
290 291
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
292
	subdomain->setLevel(meshLevel);
293 294
      }
    }
295

296 297 298 299
    primalDofMap.init(levelData, feSpaces, uniqueFe);   
    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
    lagrangeMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
300

301
    if (fullInterface != NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
302
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
303

304
    if (fetiPreconditioner != FETI_NONE) {
305 306 307
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

308
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
309
    }
310 311 312 313 314 315 316
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
317 318
    double timeCounter = MPI::Wtime();

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

323 324
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

325 326 327 328
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
329
    if (fetiPreconditioner != FETI_NONE)
330 331
      interiorDofMap.clear();

332 333
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
334

335 336
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
337
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
338
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
339 340
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
341

342 343 344
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
345 346
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

347 348 349 350 351 352
    if (fullInterface != NULL) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

353 354
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
355
    
356 357
      createPrimals(feSpace);  

358 359 360
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
361

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
362
      createIndexB(feSpace);     
363
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
364 365 366 367

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
368
    if (fetiPreconditioner != FETI_NONE)
369
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
370

371 372 373
    if (fullInterface != NULL)
      interfaceDofMap.update();

374 375 376 377
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
378

379 380 381
    lagrangeMap.update();


382 383 384 385 386 387 388 389 390 391 392 393
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

394
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
395 396 397 398 399 400 401 402 403
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	tmp = rStartInterior;

      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429
      if (feSpace == fullInterface) {
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
	MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
	    primalDofMap[feSpace].nRankDofs, 
	    primalDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
	    dualDofMap[feSpace].nRankDofs, 
	    dualDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
	    lagrangeMap[feSpace].nRankDofs, 
	    lagrangeMap[feSpace].nOverallDofs);

	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
430
    }
431

Thomas Witkowski's avatar
Thomas Witkowski committed
432
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
433 434 435 436 437 438 439 440 441

    if (enableStokesMode) {
      MSG("WARNING: MAKE THIS MORE GENERAL!!!!!\n");
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 0);
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 1);
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, 2);
    } else {	  
      subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
442 443

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
444
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
445
      timeCounter = MPI::Wtime() - timeCounter;
446
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
447 448
	  timeCounter);
    }
449 450 451
  }


452
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
453
  {
454
    FUNCNAME("PetscSolverFeti::createPrimals()");  
455

456 457 458
    if (feSpace == fullInterface)
      return;

459 460 461
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
462
    // Set of DOF indices that are considered to be primal variables.
463
    DofContainerSet& vertices = 
464
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
465 466

    DofIndexSet primals;
467
    for (DofContainerSet::iterator it = vertices.begin(); 
468 469
	 it != vertices.end(); ++it) {
      double e = 1e-8;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
470
      if (meshLevel == 0) {
471
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
472
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
473 474 475
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
476 477 478 479 480 481 482
	if (fabs(c[0]) < e || fabs(c[1]) < e ||
	    fabs(c[0] - 25.0) < e || fabs(c[1] - 25.0) < e ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 12.5) < e)) {
	  MSG("PRIMAL COORD %f %f\n", c[0], c[1]);
	  primals.insert(**it);
	}
      }
483
    }
484

485 486 487 488

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

489
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
490 491
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
492
      else
493
  	primalDofMap[feSpace].insertNonRankDof(*it);
494 495 496
  }


497
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
498 499
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
500

501 502 503
    if (feSpace == fullInterface)
      return;

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

    DofContainer allBoundaryDofs;
507
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
508 509 510 511

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
512 513
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
514 515 516 517
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
518 519 520
  }

  
521 522 523 524 525 526 527 528 529 530 531 532
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

    if (feSpace != fullInterface)
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
533 534
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it)) {
	MSG("INSERT INTERFACE RANK DOF %d\n", **it);
Thomas Witkowski's avatar
Thomas Witkowski committed
535
	interfaceDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Thomas Witkowski committed
536 537
      } else {
	MSG("INSERT INTERFACE NON-RANK DOF %d\n", **it);
Thomas Witkowski's avatar
Thomas Witkowski committed
538
	interfaceDofMap[feSpace].insertNonRankDof(**it);
Thomas Witkowski's avatar
Thomas Witkowski committed
539
      }
540 541 542
  }


543 544 545 546
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

547 548 549
    if (feSpace == fullInterface)
      return;

550 551
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
552 553 554
    // Stores for all rank owned communication DOFs, if the counterpart is
    // a rank owned DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFs that fulfill this requirenment.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
555
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
556

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
557
    if (meshLevel > 0) {
558
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
559 560 561 562 563 564 565 566 567

      for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank()) {

	vector<int> subdomainRankDofs;
	subdomainRankDofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
568
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
569
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
570
	  else
571 572 573 574 575 576 577 578 579 580 581 582 583 584 585
	    subdomainRankDofs.push_back(0);
	}

	stdMpi.send(it.getRank(), subdomainRankDofs);
      }	     

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace); 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
586 587 588 589 590 591
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex()))
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
592

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
593 594 595 596 597 598 599
    if (dualDofMap[feSpace].nLocalDofs == 0)
      return;


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

600
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
601 602 603 604 605 606 607
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(feSpace, it.getDofIndex())) {
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

608
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
609
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
610
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
611 612 613 614
	}
      }
    }

615 616 617 618

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

619
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
620

621
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
622 623
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
624
	if (!isPrimal(feSpace, it.getDofIndex()))
625 626 627
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
628 629 630

    stdMpi.updateSendDataSize();

631
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
632
	 !it.end(); it.nextRank()) {
633
      bool recvFromRank = false;
634
      for (; !it.endDofIter(); it.nextDof()) {
635
	if (!isPrimal(feSpace, it.getDofIndex())) {
636 637 638
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
639 640 641
	    recvFromRank = true;
	    break;
	  }
642
	}
643
      }
644 645

      if (recvFromRank)
646
	stdMpi.recv(it.getRank());
647
    }
648

649 650
    stdMpi.startCommunication();

651
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
652
	 !it.end(); it.nextRank()) {
653
      int i = 0;
654
      for (; !it.endDofIter(); it.nextDof())
655
	if (!isPrimal(feSpace, it.getDofIndex()))
656 657 658
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
659 660
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
661 662
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
663 664
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
665

666 667 668
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

669
    int nRankLagrange = 0;
670
    DofMap& dualMap = dualDofMap[feSpace].getMap();
671
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
672 673
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
674
	int degree = boundaryDofRanks[feSpace][it->first].size();
675
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
676
      } else {
677
	lagrangeMap[feSpace].insertNonRankDof(it->first);
678 679
      }
    }
680
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
681 682 683
  }


684
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
685
  {
686
    FUNCNAME("PetscSolverFeti::createIndexB()");
687

688
    DOFAdmin* admin = feSpace->getAdmin();
689 690 691 692

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

694
    int nLocalInterior = 0;    
695
    for (int i = 0; i < admin->getUsedSize(); i++) {
696
      if (admin->isDofFree(i) == false && 
697
	  isPrimal(feSpace, i) == false &&
698 699
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
700 701
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
702

703 704 705 706 707 708 709 710 711
	  if (fetiPreconditioner != FETI_NONE)
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
712 713 714

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
715
	}
716
      }
717
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
718
    
719 720
    // === And finally, add the global indicies of all dual nodes. ===

721 722 723 724 725 726 727 728 729 730
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      if (meshLevel == 0)
	localDofMap[feSpace].insertRankDof(it->first);
      else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
731 732 733
  }


734
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
735 736 737
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
738
    double wtime = MPI::Wtime();
739
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
740

741 742
    // === Create distributed matrix for Lagrange constraints. ===

743 744 745 746 747
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
748
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
749

750 751 752 753 754 755 756
    // === 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
757
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
758
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
759

760
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
761
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
762 763 764
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
765
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
766

Thomas Witkowski's avatar
Thomas Witkowski committed
767
	// Copy set of all ranks that contain this dual node.
768 769
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
770 771
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
772 773

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
774 775 776 777
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
778 779
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
780

781
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
782
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
783
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
784
	    index++;	      
785 786 787 788 789 790 791
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
792

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
793 794
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
795 796
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
797
    }
798 799 800
  }


801
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
802
  {
803
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
804

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

Thomas Witkowski's avatar
Thomas Witkowski committed
808
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
809
      
810
      VecCreateMPI(mpiCommGlobal, 
811
		   localDofMap.getRankDofs(), 
812
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
813
		   &(schurPrimalData.tmp_vec_b));
814
      VecCreateMPI(mpiCommGlobal, 
815 816
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
817 818
		   &(schurPrimalData.tmp_vec_primal));

819
      MatCreateShell(mpiCommGlobal,
820 821 822 823
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
824 825 826 827 828
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
829
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
830
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
831 832
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
833 834
      KSPSetFromOptions(ksp_schur_primal);
    } else {
835 836
      MSG("Create direct schur primal solver!\n");

837 838
      double wtime = MPI::Wtime();

839 840 841
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

842 843 844
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
845

Thomas Witkowski's avatar
Thomas Witkowski committed
846
      Mat matBPi;
847 848 849
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
850 851 852
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

853

854 855 856
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed