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

    return 0;
  }


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

92 93 94 95
    return 0;
  }


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

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

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


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

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

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

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

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


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

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

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

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

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

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

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


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

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

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

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

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


    // === K_DD ===

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

197

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

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

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

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

210 211

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

    return 0;
  }


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

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

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

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

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


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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

340 341
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
342
    
343 344
      createPrimals(feSpace);  

345 346 347
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
348

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

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

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

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

366 367 368
    lagrangeMap.update();


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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
413 414
    subdomain->setDofMapping(&localDofMap);
    subdomain->setCoarseSpaceDofMapping(&primalDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
415 416

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


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

429 430 431
    if (feSpace == fullInterface)
      return;

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

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

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

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

458 459 460 461

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

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


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

474 475 476
    if (feSpace == fullInterface)
      return;

477 478 479
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
480
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
481 482 483 484

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

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


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

517 518 519
    if (feSpace == fullInterface)
      return;

520 521
    boundaryDofRanks[feSpace].clear();

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

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

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

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

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

584 585 586 587

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

588
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
589

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
615
	stdMpi.recv(it.getRank());
616
    }
617

618 619
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
634

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

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


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

657
    DOFAdmin* admin = feSpace->getAdmin();
658 659 660 661

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

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

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

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

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


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

Thomas Witkowski's avatar
Thomas Witkowski committed
707 708
    double wtime = MPI::Wtime();

709 710
    // === Create distributed matrix for Lagrange constraints. ===

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

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

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

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

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

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

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

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


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

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

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

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

805 806
      double wtime = MPI::Wtime();

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