PetscSolverFeti.cc 55.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscSolverFeti.h"
16
#include "parallel/PetscSolverFetiStructs.h"
17 18
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
19
#include "parallel/PetscSolverGlobalMatrix.h"
20
#include "io/VtkWriter.h"
21 22 23 24 25

namespace AMDiS {

  using namespace std;

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


32 33 34 35 36 37 38
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
39
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
40

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

    return 0;
  }


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

57 58
    //    F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
59

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

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

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

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

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

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

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

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

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

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

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

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

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

95 96 97 98
    return 0;
  }


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

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

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


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

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

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

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

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


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

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

136
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
137
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
138 139 140 141 142 143 144 145 146 147
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

    VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);


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

    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

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

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


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

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

161 162 163 164 165 166 167 168 169 170 171 172 173 174
    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
175 176


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

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

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

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

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


    // === K_DD ===

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

200

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

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

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

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

213 214

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

    return 0;
  }


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

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

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

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

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


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
272 273
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
274

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

286 287 288 289
    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
290

291 292
    if (fullInterface != NULL)
      interfaceDofMap.init(levelData, fullInterface);
Thomas Witkowski's avatar
Thomas Witkowski committed
293

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

298
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
299
    }
300 301 302 303 304 305 306
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
307 308
    double timeCounter = MPI::Wtime();

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

313 314
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

315 316 317 318
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
319
    if (fetiPreconditioner != FETI_NONE)
320 321
      interiorDofMap.clear();

322 323
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
324

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

332 333 334
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
335 336
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

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

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

348 349 350
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
351

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
352
      createIndexB(feSpace);     
353
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
354 355 356 357

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
358
    if (fetiPreconditioner != FETI_NONE)
359
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
360

361 362 363
    if (fullInterface != NULL)
      interfaceDofMap.update();

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

369 370 371
    lagrangeMap.update();


372 373 374 375 376 377 378 379 380 381 382 383
    // === ===

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

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

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

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

407
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
408 409
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
410

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

Thomas Witkowski's avatar
Thomas Witkowski committed
416 417
    subdomain->setDofMapping(&localDofMap);
    subdomain->setCoarseSpaceDofMapping(&primalDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
418 419

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


428
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
429
  {
430
    FUNCNAME("PetscSolverFeti::createPrimals()");  
431

432 433 434
    if (feSpace == fullInterface)
      return;

435 436 437
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

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

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

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

461 462 463 464

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

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


473
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
474 475
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
476

477 478 479
    if (feSpace == fullInterface)
      return;

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


516 517 518 519
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

520 521 522
    if (feSpace == fullInterface)
      return;

523 524
    boundaryDofRanks[feSpace].clear();

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

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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
566 567 568 569 570 571 572
    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)).         ===

573
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
574 575 576 577 578 579 580
    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);

581
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
582
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
583
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
584 585 586 587
	}
      }
    }

588 589 590 591

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

592
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
593

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
619
	stdMpi.recv(it.getRank());
620
    }
621

622 623
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
638

639 640 641
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


657
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
658
  {
659
    FUNCNAME("PetscSolverFeti::createIndexB()");
660

661
    DOFAdmin* admin = feSpace->getAdmin();
662 663 664 665

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

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

676 677 678 679 680 681 682 683 684
	  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);
685 686 687

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

694 695 696 697 698 699 700 701 702 703
    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);
      }
704 705 706
  }


707
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
708 709 710
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
711
    double wtime = MPI::Wtime();
712
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
713

714 715
    // === Create distributed matrix for Lagrange constraints. ===

716 717 718 719 720
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
721
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
722

723 724 725 726 727 728 729
    // === 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
730
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
731
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
732

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

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

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

754
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
755
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
756
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
757
	    index++;	      
758 759 760 761 762 763 764
	  }
	}
      }
    }

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

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


774
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
775
  {
776
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
777

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

Thomas Witkowski's avatar
Thomas Witkowski committed
781
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
782
      
783
      VecCreateMPI(mpiCommGlobal, 
784
		   localDofMap.getRankDofs(), 
785
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
786
		   &(schurPrimalData.tmp_vec_b));
787
      VecCreateMPI(mpiCommGlobal, 
788 789
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
790 791
		   &(schurPrimalData.tmp_vec_primal));

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

810 811
      double wtime = MPI::Wtime();

812 813 814
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

815 816 817
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
818

Thomas Witkowski's avatar
Thomas Witkowski committed
819
      Mat matBPi;
820 821 822
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
823 824 825
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

826

827 828 829
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
830

831 832
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

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