PetscSolverFeti.cc 54.4 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 266
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

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

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

281 282
    primalDofMap.init(levelData,
		      feSpaces, meshDistributor->getFeSpaces(), 
283
		      true, true);
284 285
    primalDofMap.setComputeMatIndex(true);
    
286 287
    dualDofMap.init(levelData, 
		    feSpaces, meshDistributor->getFeSpaces(),
288
		    false, false);
289 290
    dualDofMap.setComputeMatIndex(true);

291 292 293 294 295 296 297 298 299 300 301
    if (meshLevel == 0) {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       false, false);
      localDofMap.setComputeMatIndex(true);
    } else {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       true, true);
      localDofMap.setComputeMatIndex(true);
    }
302

Thomas Witkowski's avatar
Thomas Witkowski committed
303 304 305 306 307 308
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
		     true, true);
    lagrangeMap.setComputeMatIndex(true);


309
    if (fetiPreconditioner != FETI_NONE) {
310 311 312
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

313 314
      interiorDofMap.init(levelData, 
			  feSpaces,  meshDistributor->getFeSpaces(),
315
			  false, false);
316 317
      interiorDofMap.setComputeMatIndex(true);
    }
318 319 320 321 322 323 324
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
325 326
    double timeCounter = MPI::Wtime();

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

331 332
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

333 334 335 336
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
337
    if (fetiPreconditioner != FETI_NONE)
338 339
      interiorDofMap.clear();

340 341
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
342

343 344
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
345
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
346
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
347 348
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
349

350 351 352
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
353 354
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

355 356
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
357

358 359
      createPrimals(feSpace);  

360
      createDuals(feSpace);      
361

362 363
      createIndexB(feSpace);
    }
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 374
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
375

376 377 378
    lagrangeMap.update();


379 380 381 382 383 384 385 386 387 388 389 390
    // === ===

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

391
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
392 393 394 395 396 397 398 399 400
			   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);
    }

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

406
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
407 408
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
409
      
410
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
411 412
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
413

414
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
415 416
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
417

418
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
419 420
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
421
    }
422

Thomas Witkowski's avatar
Thomas Witkowski committed
423 424 425
    subdomain->setDofMapping(&localDofMap, &primalDofMap);

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
426
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
427
      timeCounter = MPI::Wtime() - timeCounter;
428
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
429 430
	  timeCounter);
    }
431 432 433
  }


434
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
435
  {
436
    FUNCNAME("PetscSolverFeti::createPrimals()");  
437

438 439 440
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
441
    // Set of DOF indices that are considered to be primal variables.
442
    DofContainerSet& vertices = 
443
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
444 445

    DofIndexSet primals;
446
    for (DofContainerSet::iterator it = vertices.begin(); 
447 448
	 it != vertices.end(); ++it) {
      double e = 1e-8;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
449
      if (meshLevel == 0) {
450
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
451
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
452 453 454
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
455 456 457 458 459 460 461
	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);
	}
      }
462
    }
463

464 465 466 467

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

468
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
469 470
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
471
      else
472
  	primalDofMap[feSpace].insertNonRankDof(*it);
473 474 475
  }


476
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
477 478
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
479

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

    DofContainer allBoundaryDofs;
483
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
484 485 486 487

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
488 489
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
490 491 492 493
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
494 495 496 497 498 499 500
  }

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

501 502
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
503 504 505
    // 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
506
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
507

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
508
    if (meshLevel > 0) {
509
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
510 511 512 513 514 515 516 517 518

      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
519
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
520
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
521
	  else
522 523 524 525 526 527 528 529 530 531 532 533 534 535 536
	    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
537 538 539 540 541 542
	   !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());
    }
543

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
544 545 546 547 548 549 550 551 552 553 554 555 556 557
    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);

558
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
559
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
560
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
561 562 563 564
	}
      }
    }

565 566 567 568

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

569
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
570

571
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
572 573
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
574
	if (!isPrimal(feSpace, it.getDofIndex()))
575 576 577
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
578 579 580

    stdMpi.updateSendDataSize();

581
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
582
	 !it.end(); it.nextRank()) {
583
      bool recvFromRank = false;
584
      for (; !it.endDofIter(); it.nextDof()) {
585
	if (!isPrimal(feSpace, it.getDofIndex())) {
586 587 588
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
589 590 591
	    recvFromRank = true;
	    break;
	  }
592
	}
593
      }
594 595

      if (recvFromRank)
596
	stdMpi.recv(it.getRank());
597
    }
598

599 600
    stdMpi.startCommunication();

601
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
602
	 !it.end(); it.nextRank()) {
603
      int i = 0;
604
      for (; !it.endDofIter(); it.nextDof())
605
	if (!isPrimal(feSpace, it.getDofIndex()))
606 607 608
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
609 610
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
611 612
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
613 614
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
615

616 617 618
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

619
    int nRankLagrange = 0;
620
    DofMap& dualMap = dualDofMap[feSpace].getMap();
621
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
622 623
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
624
	int degree = boundaryDofRanks[feSpace][it->first].size();
625
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
626
      } else {
627
	lagrangeMap[feSpace].insertNonRankDof(it->first);
628 629
      }
    }
630
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
631 632 633
  }


634
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
635
  {
636
    FUNCNAME("PetscSolverFeti::createIndexB()");
637

638
    DOFAdmin* admin = feSpace->getAdmin();
639 640 641 642

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

644
    int nLocalInterior = 0;    
645
    for (int i = 0; i < admin->getUsedSize(); i++) {
646
      if (admin->isDofFree(i) == false && 
647
	  isPrimal(feSpace, i) == false &&
648
	  isDual(feSpace, i) == false) {
649 650
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
651

652 653 654 655 656 657 658 659 660
	  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);
661 662 663

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
664
	}
665
      }
666
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
667
    
668 669
    // === And finally, add the global indicies of all dual nodes. ===

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


683
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
684 685 686
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
687 688
    double wtime = MPI::Wtime();

689 690
    // === Create distributed matrix for Lagrange constraints. ===

691 692 693 694 695
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
696
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
697

698 699 700 701 702 703 704
    // === 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
705
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
706
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
707

708
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
709
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
710 711 712
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
713
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
714

Thomas Witkowski's avatar
Thomas Witkowski committed
715
	// Copy set of all ranks that contain this dual node.
716 717
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
718 719
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
720 721

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
722 723 724 725
	
	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
726 727
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
728

729
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
730
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
731
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
732
	    index++;	      
733 734 735 736 737 738 739
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
741 742
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
743 744
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
745
    }
746 747 748
  }


749
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
750
  {
751
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
752

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

Thomas Witkowski's avatar
Thomas Witkowski committed
756
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
757
      
758
      VecCreateMPI(mpiCommGlobal, 
759
		   localDofMap.getRankDofs(), 
760
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
761
		   &(schurPrimalData.tmp_vec_b));
762
      VecCreateMPI(mpiCommGlobal, 
763 764
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
765 766
		   &(schurPrimalData.tmp_vec_primal));

767
      MatCreateShell(mpiCommGlobal,
768 769 770 771
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
772 773 774 775 776
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
777
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
778
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
779 780
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
781 782
      KSPSetFromOptions(ksp_schur_primal);
    } else {
783 784
      MSG("Create direct schur primal solver!\n");

785 786
      double wtime = MPI::Wtime();

787 788 789
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

790 791 792
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
793

Thomas Witkowski's avatar
Thomas Witkowski committed
794
      Mat matBPi;
795 796 797
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
798 799 800
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

801

802 803 804
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
805

806 807
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

808
      for (int i = 0; i < nRowsRankB; i++) {