PetscSolverFeti.cc 54.9 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/PetscSolverGlobalMatrix.h"
19
#include "io/VtkWriter.h"
20 21 22 23 24

namespace AMDiS {

  using namespace std;

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


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

40 41 42 43
    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);
44 45 46 47 48 49 50 51 52
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


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

55 56
    //    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)
57

58 59
    double wtime = MPI::Wtime();

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

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

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

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

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

73
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
74 75

    wtime01 = MPI::Wtime();
76
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
77 78
    FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);

79
    MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
80 81

    wtime01 = MPI::Wtime();
82
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
83 84
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

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

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

89 90
    FetiTimings::fetiSolve += (MPI::Wtime() - wtime);

91 92 93 94
    return 0;
  }


95
  // y = PC * x
96
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
97
  {
98 99
    double wtime = MPI::Wtime();

100
    // Get data for the preconditioner
101 102
    void *ctx;
    PCShellGetContext(pc, &ctx);
103
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
104

105
    // Multiply with scaled Lagrange constraint matrix.
106 107 108
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


109
    // === Restriction of the B nodes to the boundary nodes. ===
110

111 112 113 114
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
115

116 117 118 119
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

120 121 122
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
123 124

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


128
    // === K_DD - K_DI inv(K_II) K_ID ===
129

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

132
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
133
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
134 135 136 137 138 139 140 141 142 143
    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);

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

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

155 156
    FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);

157 158 159 160 161 162 163 164 165 166 167 168 169 170
    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);
171 172


173
    // === Restriction of the B nodes to the boundary nodes. ===
174

175 176 177 178
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
179

180
    PetscScalar *local_b, *local_duals;
181
    VecGetArray(data->tmp_vec_b, &local_b);
182
    VecGetArray(data->tmp_vec_duals0, &local_duals);
183

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

    VecRestoreArray(data->tmp_vec_b, &local_b);
189 190 191 192 193 194 195
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

196

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

199 200 201
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

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

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

209 210

    // Multiply with scaled Lagrange constraint matrix.
211 212 213 214 215 216
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


217 218
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
219
      schurPrimalSolver(0),
220
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
221
      subdomain(NULL),
222 223
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
224 225
      nGlobalOverallInterior(0),
      printTimings(false)
226 227 228 229 230 231
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

    Parameters::get("parallel->multi level test", multiLevelTest);
251 252
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
253 254

    Parameters::get("parallel->print timings", printTimings);
255 256 257
  }


258
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
259
  {
260 261
    FUNCNAME("PetscSolverFeti::initialize()");

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

265
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
266
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
267

Thomas Witkowski's avatar
Thomas Witkowski committed
268 269
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
270

271
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
272
	subdomain->setMeshDistributor(meshDistributor, 
273
				      mpiCommGlobal, mpiCommLocal);
274
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
275
	subdomain->setMeshDistributor(meshDistributor, 
276 277
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
278
	subdomain->setLevel(meshLevel);
279 280
      }
    }
281

282 283 284 285
    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
286

287 288
    if (fullInterface != NULL)
      interfaceDofMap.init(levelData, fullInterface);
Thomas Witkowski's avatar
Thomas Witkowski committed
289

290
    if (fetiPreconditioner != FETI_NONE) {
291 292 293
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

294
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
295
    }
296 297 298 299 300 301 302
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
303 304
    double timeCounter = MPI::Wtime();

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

309 310
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

311 312 313 314
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
315
    if (fetiPreconditioner != FETI_NONE)
316 317
      interiorDofMap.clear();

318 319
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
320

321 322
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
323
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
324
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
325 326
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
327

328 329 330
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
331 332
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

333 334 335 336 337 338
    if (fullInterface != NULL) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

339 340
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
341

342 343
      createPrimals(feSpace);  

344 345 346
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
347

348 349
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
350 351 352 353

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
354
    if (fetiPreconditioner != FETI_NONE)
355
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
356

357 358 359
    if (fullInterface != NULL)
      interfaceDofMap.update();

360 361 362 363
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
364

365 366 367
    lagrangeMap.update();


368 369 370 371 372 373 374 375 376 377 378 379
    // === ===

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

380
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
381 382 383 384 385 386 387 388 389
			   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);
    }

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

395
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
396 397
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
398
      
399
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
400 401
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
402

403
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
404 405
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
406

407
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
408 409
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
410
    }
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412 413 414
    subdomain->setDofMapping(&localDofMap, &primalDofMap);

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
415
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
416
      timeCounter = MPI::Wtime() - timeCounter;
417
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
418 419
	  timeCounter);
    }
420 421 422
  }


423
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
424
  {
425
    FUNCNAME("PetscSolverFeti::createPrimals()");  
426

427 428 429
    if (feSpace == fullInterface)
      return;

430 431 432
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    // Set of DOF indices that are considered to be primal variables.
434
    DofContainerSet& vertices = 
435
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
436 437

    DofIndexSet primals;
438
    for (DofContainerSet::iterator it = vertices.begin(); 
439 440
	 it != vertices.end(); ++it) {
      double e = 1e-8;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
441
      if (meshLevel == 0) {
442
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
443
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
444 445 446
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
447 448 449 450 451 452 453
	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);
	}
      }
454
    }
455

456 457 458 459

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

460
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
461 462
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
463
      else
464
  	primalDofMap[feSpace].insertNonRankDof(*it);
465 466 467
  }


468
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
469 470
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
471

472 473 474
    if (feSpace == fullInterface)
      return;

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

    DofContainer allBoundaryDofs;
478
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
479 480 481 482

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
483 484
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
485 486 487 488
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
489 490 491
  }

  
492 493 494 495 496 497 498 499 500 501 502 503
  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
504 505
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
	interfaceDofMap[feSpace].insertRankDof(**it);
506
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
507
	interfaceDofMap[feSpace].insertNonRankDof(**it);
508 509 510
  }


511 512 513 514
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

515 516 517
    if (feSpace == fullInterface)
      return;

518 519
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
520 521 522
    // 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
523
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
524

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
525
    if (meshLevel > 0) {
526
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
527 528 529 530 531 532 533 534 535

      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
536
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
537
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
538
	  else
539 540 541 542 543 544 545 546 547 548 549 550 551 552 553
	    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
554 555 556 557 558 559
	   !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());
    }
560

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
561 562 563 564 565 566 567 568 569 570 571 572 573 574
    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)).         ===

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

575
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
576
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
577
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
578 579 580 581
	}
      }
    }

582 583 584 585

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

586
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
587

588
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
589 590
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
591
	if (!isPrimal(feSpace, it.getDofIndex()))
592 593 594
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
595 596 597

    stdMpi.updateSendDataSize();

598
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
599
	 !it.end(); it.nextRank()) {
600
      bool recvFromRank = false;
601
      for (; !it.endDofIter(); it.nextDof()) {
602
	if (!isPrimal(feSpace, it.getDofIndex())) {
603 604 605
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
606 607 608
	    recvFromRank = true;
	    break;
	  }
609
	}
610
      }
611 612

      if (recvFromRank)
613
	stdMpi.recv(it.getRank());
614
    }
615

616 617
    stdMpi.startCommunication();

618
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
619
	 !it.end(); it.nextRank()) {
620
      int i = 0;
621
      for (; !it.endDofIter(); it.nextDof())
622
	if (!isPrimal(feSpace, it.getDofIndex()))
623 624 625
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
626 627
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
628 629
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
630 631
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
632

633 634 635
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

636
    int nRankLagrange = 0;
637
    DofMap& dualMap = dualDofMap[feSpace].getMap();
638
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
639 640
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
641
	int degree = boundaryDofRanks[feSpace][it->first].size();
642
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
643
      } else {
644
	lagrangeMap[feSpace].insertNonRankDof(it->first);
645 646
      }
    }
647
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
648 649 650
  }


651
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
652
  {
653
    FUNCNAME("PetscSolverFeti::createIndexB()");
654

655
    DOFAdmin* admin = feSpace->getAdmin();
656 657 658 659

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

661
    int nLocalInterior = 0;    
662
    for (int i = 0; i < admin->getUsedSize(); i++) {
663
      if (admin->isDofFree(i) == false && 
664
	  isPrimal(feSpace, i) == false &&
665 666
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
667 668
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
669

670 671 672 673 674 675 676 677 678
	  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);
679 680 681

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
682
	}
683
      }
684
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
685
    
686 687
    // === And finally, add the global indicies of all dual nodes. ===

688 689 690 691 692 693 694 695 696 697
    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);
      }
698 699 700
  }


701
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
702 703 704
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
705 706
    double wtime = MPI::Wtime();

707 708
    // === Create distributed matrix for Lagrange constraints. ===

709 710 711 712 713
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
714
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
715

716 717 718 719 720 721 722
    // === 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
723
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
724
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
725

726
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
727
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
728 729 730
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
731
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
732

Thomas Witkowski's avatar
Thomas Witkowski committed
733
	// Copy set of all ranks that contain this dual node.
734 735
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
736 737
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
738 739

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
740 741 742 743
	
	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
744 745
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
746

747
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
748
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
749
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
750
	    index++;	      
751 752 753 754 755 756 757
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
759 760
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
761 762
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
763
    }
764 765 766
  }


767
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
768
  {
769
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
770

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

Thomas Witkowski's avatar
Thomas Witkowski committed
774
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
775
      
776
      VecCreateMPI(mpiCommGlobal, 
777
		   localDofMap.getRankDofs(), 
778
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
779
		   &(schurPrimalData.tmp_vec_b));
780
      VecCreateMPI(mpiCommGlobal, 
781 782
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
783 784
		   &(schurPrimalData.tmp_vec_primal));

785
      MatCreateShell(mpiCommGlobal,
786 787 788 789
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
790 791 792 793 794
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
795
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
796
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
797 798
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
799 800
      KSPSetFromOptions(ksp_schur_primal);
    } else {
801 802
      MSG("Create direct schur primal solver!\n");

803 804
      double wtime = MPI::Wtime();

805 806 807
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

808 809 810
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
811