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


13
#include "AMDiS.h"
14
#include "parallel/PetscSolverFeti.h"
15
#include "parallel/PetscSolverFetiStructs.h"
16 17
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
18
#include "parallel/PetscSolverGlobalMatrix.h"
19
#include "io/VtkWriter.h"
20 21 22 23 24

namespace AMDiS {

  using namespace std;

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


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

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

40 41 42 43
    MatMult(data->subSolver->getMatIntCoarse(), x, data->tmp_vec_b);
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarseCoarse(), x, y);
44 45 46 47 48 49 50 51 52
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

91 92 93 94
    return 0;
  }


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

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

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


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

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

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

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

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


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

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

132
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
133
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
134 135 136 137 138 139 140 141 142 143
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

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


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

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

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

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


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

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

157 158 159 160 161 162 163 164 165 166 167 168 169 170
    return 0;
  }


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

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


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

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

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

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

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


    // === K_DD ===

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

196

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

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

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

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

209 210

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

    return 0;
  }


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

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

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

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

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


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

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

265 266
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

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

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

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

291 292
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
293
		     true, true);
294
    lagrangeMap.setComputeMatIndex(true);
295 296 297 298 299 300 301 302 303 304 305 306

    if (meshLevel == 0) {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       false, false);
      localDofMap.setComputeMatIndex(true);
    } else {
      localDofMap.init(levelData, 
		       feSpaces, meshDistributor->getFeSpaces(),
		       true, true);
      localDofMap.setComputeMatIndex(true);
    }
307

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

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


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

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

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

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

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

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

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

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

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

357 358
      createPrimals(feSpace);  

359
      createDuals(feSpace);      
360

361 362
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
363 364 365 366

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

370 371 372 373
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
374

375 376 377
    lagrangeMap.update();


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

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

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

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

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

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

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

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

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


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

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

439 440
    /// Set of DOF indices that are considered to be primal variables.
    
441
    DofContainerSet& vertices = 
442
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
443 444

    DofIndexSet primals;
445
    for (DofContainerSet::iterator it = vertices.begin(); 
446 447 448
	 it != vertices.end(); ++it) {
      WorldVector<double> c;
      feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
449

450
      double e = 1e-8;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
451
      if (meshLevel == 0) {
452
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
453 454 455 456 457 458 459 460
      } else {
	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);
	}
      }
461
    }
462

463 464 465 466

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

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


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

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

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

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

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

500 501
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
502 503 504 505
    // Stores for all rank own communication DOFs, if the counterpart is
    // a rank owmed DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFS that fulfill this requirenment.
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
506

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

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

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

557 558 559 560 561
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex()))) {
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
 	  }

562 563 564 565
	}
      }
    }

566 567 568 569

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

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

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

    stdMpi.updateSendDataSize();

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

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

600 601
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
616

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

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


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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
741 742 743 744

    if (printTimings)
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
745 746 747
  }


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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
800
      Mat matPrimal;
801

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

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

808
      for (int i = 0; i < nRowsRankB; i++) {
809
	PetscInt row = localDofMap.getStartDofs() + i;
Thomas Witkowski's avatar
Thomas Witkowski committed
810
	MatGetRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
811 812 813 814 815

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

Thomas Witkowski's avatar
Thomas Witkowski committed
816
	MatRestoreRow(subdomain->getMatIntCoarse(), row, &nCols, &cols, &values);
817 818
      }

819
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
820
		primalDofMap.getLocalDofs())
821
	("Should not happen!\n");
822

823 824 825
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
826
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
827 828 829 830 831 832 833 834

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
835
	subdomain->solve(tmpVec, tmpVec);
836

Thomas Witkowski's avatar
Thomas Witkowski committed
837
	PetscScalar *tmpValues;
838
	VecGetArray(tmpVec, &tmpValues);
839
	for (int i  = 0; i < nRowsRankB; i++)