Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

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


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

namespace AMDiS {

  using namespace std;

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


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

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

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

    return 0;
  }


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

95 96 97 98
    return 0;
  }


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

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

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


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

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

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

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

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


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

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

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

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


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

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

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

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


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

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

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


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

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


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

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

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

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

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


    // === K_DD ===

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

200

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

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

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

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

213 214

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

    return 0;
  }


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

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

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

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

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


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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

359 360 361
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
362

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

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

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

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

380 381 382
    lagrangeMap.update();


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

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

395
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
396 397 398 399 400 401 402 403 404
			   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);
    }

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

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

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

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

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

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


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

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

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

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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
477 478 479 480 481 482 483
	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);
	}
      }
484
    }
485

486 487 488 489

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

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


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

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

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

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

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

  
522 523 524 525 526 527 528 529 530 531 532 533
  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
534 535
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
	interfaceDofMap[feSpace].insertRankDof(**it);
536
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
537
	interfaceDofMap[feSpace].insertNonRankDof(**it);
538 539 540
  }


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

545 546 547
    if (feSpace == fullInterface)
      return;

548 549
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
550 551 552
    // 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
553
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
554

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

      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
566
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
567
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
568
	  else
569 570 571 572 573 574 575 576 577 578 579 580 581 582 583
	    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
584 585 586 587 588 589
	   !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());
    }
590

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
591 592 593 594 595 596 597
    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)).         ===

598
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
599 600 601 602 603 604 605
    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);

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

613 614 615 616

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

617
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
618

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
644
	stdMpi.recv(it.getRank());
645
    }
646

647 648
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
663

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

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


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

686
    DOFAdmin* admin = feSpace->getAdmin();
687 688 689 690

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

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

701 702 703 704 705 706 707 708 709
	  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);
710 711 712

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

719 720 721 722 723 724 725 726 727 728
    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);
      }
729 730 731
  }


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

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

739 740
    // === Create distributed matrix for Lagrange constraints. ===

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

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

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
772 773 774 775
	
	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
776 777
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
778

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

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

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


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

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

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

817
      MatCreateShell(mpiCommGlobal,
818 819 820 821
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
822 823 824 825 826
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
827
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);