PetscSolverFeti.cc 46.8 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 31
  // 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);
32
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
33

34 35 36 37
    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);
38 39 40 41 42 43 44 45 46
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
47 48
    //    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)
49 50 51

    void *ctx;
    MatShellGetContext(mat, &ctx);
52
    FetiData* data = static_cast<FetiData*>(ctx);
53

Thomas Witkowski's avatar
Thomas Witkowski committed
54
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
55
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
57

58
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
59
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
60
    MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
61
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
63

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
65 66 67 68 69

    return 0;
  }


70
  // y = PC * x
71
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
72
  {
73
    // Get data for the preconditioner
74 75
    void *ctx;
    PCShellGetContext(pc, &ctx);
76
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
77

78
    // Multiply with scaled Lagrange constraint matrix.
79 80 81
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


82
    // === Restriction of the B nodes to the boundary nodes. ===
83

84 85 86 87
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
88

89 90 91 92
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

93 94 95
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
96 97

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


101
    // === K_DD - K_DI inv(K_II) K_ID ===
102

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

105
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
106
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
107 108 109 110 111 112 113 114 115 116
    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);

117 118 119
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141

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

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


144
    // === Restriction of the B nodes to the boundary nodes. ===
145

146 147 148 149
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
150

151
    PetscScalar *local_b, *local_duals;
152
    VecGetArray(data->tmp_vec_b, &local_b);
153
    VecGetArray(data->tmp_vec_duals0, &local_duals);
154

Thomas Witkowski's avatar
Thomas Witkowski committed
155 156 157
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
158 159

    VecRestoreArray(data->tmp_vec_b, &local_b);
160 161 162 163 164 165 166
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

167

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

170 171 172
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

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

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

180 181

    // Multiply with scaled Lagrange constraint matrix.
182 183 184 185 186 187
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


188 189
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
190
      schurPrimalSolver(0),
191
      multiLevelTest(false),
192
      subDomainSolver(NULL),
193 194 195
      meshLevel(0),
      rStartInterior(0),
      nGlobalOverallInterior(0)
196 197 198 199 200 201
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
202
      MSG("Create FETI-DP solver with no preconditioner!\n");
203 204
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
205
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
206 207
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
208
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
209 210
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
211 212
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
213
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
214 215 216 217 218

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

    Parameters::get("parallel->multi level test", multiLevelTest);
221 222
    if (multiLevelTest)
      meshLevel = 1;
223 224 225
  }


226
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
227
  {
228 229
    FUNCNAME("PetscSolverFeti::initialize()");

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

233 234
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

235
    if (subDomainSolver == NULL) {
236
      subDomainSolver = new PetscSolverGlobalMatrix();
237

238
      if (meshLevel == 0) {
239 240
	subDomainSolver->setMeshDistributor(meshDistributor, 
					    mpiCommGlobal, mpiCommLocal);
241
      } else {
242 243 244
	subDomainSolver->setMeshDistributor(meshDistributor, 
					    levelData.getMpiComm(meshLevel - 1),
					    levelData.getMpiComm(meshLevel));
245
	subDomainSolver->setLevel(meshLevel);
246 247
      }
    }
248

249 250
    primalDofMap.init(levelData,
		      feSpaces, meshDistributor->getFeSpaces(), 
251
		      true, true);
252 253
    primalDofMap.setComputeMatIndex(true);
    
254 255
    dualDofMap.init(levelData, 
		    feSpaces, meshDistributor->getFeSpaces(),
256
		    false, false);
257 258
    dualDofMap.setComputeMatIndex(true);

259 260
    lagrangeMap.init(levelData, 
		     feSpaces, meshDistributor->getFeSpaces(),
261
		     true, true);
262
    lagrangeMap.setComputeMatIndex(true);
263 264 265 266 267 268 269 270 271 272 273 274

    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);
    }
275 276

    if (fetiPreconditioner == FETI_DIRICHLET) {
277 278 279
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

280 281
      interiorDofMap.init(levelData, 
			  feSpaces,  meshDistributor->getFeSpaces(),
282
			  false, false);
283 284
      interiorDofMap.setComputeMatIndex(true);
    }
285 286 287 288 289 290 291 292 293 294
  }


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

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

296 297
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

298 299 300 301
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
302 303 304
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.clear();

305 306
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
307

308 309
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
310
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
311
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
312 313
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
314 315 316 317

    if (meshLevel > 0)
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

318 319
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
320

321 322
      createPrimals(feSpace);  

323
      createDuals(feSpace);      
324

325 326
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
327 328 329 330

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
331
    if (fetiPreconditioner != FETI_NONE)
332
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
333

334 335 336 337
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
338

339 340 341
    lagrangeMap.update();


342 343 344 345 346 347 348 349 350 351 352 353
    // === ===

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

354
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
355 356 357 358 359 360 361 362 363
			   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);
    }

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

369
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
370 371
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
372
      
373
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
374 375
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
376

377
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
378 379
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
380

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

    // If multi level test, inform sub domain solver about coarse space.
387
    subDomainSolver->setDofMapping(&localDofMap, &primalDofMap);
388 389 390
  }


391
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
392
  {
393
    FUNCNAME("PetscSolverFeti::createPrimals()");  
394

395 396 397
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

398 399
    /// Set of DOF indices that are considered to be primal variables.
    
400
    DofContainerSet& vertices = 
401
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
402 403

    DofIndexSet primals;
404
    for (DofContainerSet::iterator it = vertices.begin(); 
405 406 407
	 it != vertices.end(); ++it) {
      WorldVector<double> c;
      feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
408

409
      double e = 1e-8;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
410
      if (meshLevel == 0) {
411
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
412 413 414 415 416 417 418 419
      } 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);
	}
      }
420
    }
421

422 423 424 425

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

426
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
427 428
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
429
      else
430
  	primalDofMap[feSpace].insertNonRankDof(*it);
431 432 433
  }


434
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
435 436
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
437

438 439 440
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
441
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
442 443 444 445

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
446 447
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
448 449 450 451
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
452 453 454 455 456 457 458
  }

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

459 460
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
461 462 463 464
    // 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;
465

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
466
    if (meshLevel > 0) {
467
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
468 469 470 471 472 473 474 475 476

      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
477
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
478
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
479
	  else
480 481 482 483 484 485 486 487 488 489 490 491 492 493 494
	    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
495 496 497 498 499 500
	   !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());
    }
501

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
502 503 504 505 506 507 508 509 510 511 512 513 514 515
    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);

516 517 518 519 520
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex()))) {
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
 	  }

521 522 523 524
	}
      }
    }

525 526 527 528

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

529
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
530

531
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
532 533
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
534
	if (!isPrimal(feSpace, it.getDofIndex()))
535 536 537
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
538 539 540

    stdMpi.updateSendDataSize();

541
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
542
	 !it.end(); it.nextRank()) {
543
      bool recvFromRank = false;
544
      for (; !it.endDofIter(); it.nextDof()) {
545
	if (!isPrimal(feSpace, it.getDofIndex())) {
546 547 548
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
549 550 551
	    recvFromRank = true;
	    break;
	  }
552
	}
553
      }
554 555

      if (recvFromRank)
556
	stdMpi.recv(it.getRank());
557
    }
558

559 560
    stdMpi.startCommunication();

561
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
562
	 !it.end(); it.nextRank()) {
563
      int i = 0;
564
      for (; !it.endDofIter(); it.nextDof())
565
	if (!isPrimal(feSpace, it.getDofIndex()))
566 567 568
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
569 570
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
571 572
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
573 574
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
575

576 577 578
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

579
    int nRankLagrange = 0;
580
    DofMap& dualMap = dualDofMap[feSpace].getMap();
581
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
582 583
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
584
	int degree = boundaryDofRanks[feSpace][it->first].size();
585
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
586
      } else {
587
	lagrangeMap[feSpace].insertNonRankDof(it->first);
588 589
      }
    }
590
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
591 592 593
  }


594
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
595
  {
596
    FUNCNAME("PetscSolverFeti::createIndexB()");
597

598
    DOFAdmin* admin = feSpace->getAdmin();
599 600 601 602

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

604
    int nLocalInterior = 0;    
605
    for (int i = 0; i < admin->getUsedSize(); i++) {
606
      if (admin->isDofFree(i) == false && 
607
	  isPrimal(feSpace, i) == false &&
608
	  isDual(feSpace, i) == false) {
609 610
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
611

612 613 614 615 616 617 618 619 620
	  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);
621 622 623

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
624
	}
625
      }
626
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
627
    
628 629
    // === And finally, add the global indicies of all dual nodes. ===

630 631 632 633 634 635 636 637 638 639
    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);
      }
640 641 642
  }


643
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
644 645 646
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

647 648
    // === Create distributed matrix for Lagrange constraints. ===

649
    MatCreateMPIAIJ(mpiCommGlobal,
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
650 651
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
652 653 654
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

655 656 657 658 659 660 661
    // === 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
662
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
663
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
664

665
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
666
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
667 668 669
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
670
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
671

Thomas Witkowski's avatar
Thomas Witkowski committed
672
	// Copy set of all ranks that contain this dual node.
673 674
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
675 676
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
677 678

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
679 680 681 682
	
	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
683 684
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
685

686
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
687
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
688
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
689
	    index++;	      
690 691 692 693 694 695 696 697 698 699
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
  }


700
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
701
  {
702
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
703

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

707
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
708
      
709
      VecCreateMPI(mpiCommGlobal, 
710
		   localDofMap.getRankDofs(), 
711
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
712
		   &(schurPrimalData.tmp_vec_b));
713
      VecCreateMPI(mpiCommGlobal, 
714 715
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
716 717
		   &(schurPrimalData.tmp_vec_primal));

718
      MatCreateShell(mpiCommGlobal,
719 720 721 722
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
723 724 725 726 727
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
728
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
729
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
730 731
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
732 733
      KSPSetFromOptions(ksp_schur_primal);
    } else {
734 735
      MSG("Create direct schur primal solver!\n");

736 737
      double wtime = MPI::Wtime();

738 739 740
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

741 742 743
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
744

Thomas Witkowski's avatar
Thomas Witkowski committed
745
      Mat matBPi;
746
      MatCreateMPIAIJ(mpiCommGlobal,
Thomas Witkowski's avatar
Thomas Witkowski committed
747
		      nRowsRankB, nRowsRankPrimal, 
748
		      nGlobalOverallInterior, nRowsOverallPrimal,
Thomas Witkowski's avatar
Thomas Witkowski committed
749 750
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
751

752 753 754
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
755

756 757
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

758
      for (int i = 0; i < nRowsRankB; i++) {
759
	PetscInt row = localDofMap.getStartDofs() + i;
760
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
761 762 763 764 765

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

766
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
767 768
      }

769
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
770
		primalDofMap.getLocalDofs())
771
	("Should not happen!\n");
772

773 774 775
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
776
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
777 778 779 780 781 782 783 784

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

785
	subDomainSolver->solve(tmpVec, tmpVec);
786

Thomas Witkowski's avatar
Thomas Witkowski committed
787
	PetscScalar *tmpValues;
788
	VecGetArray(tmpVec, &tmpValues);
789
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
790
	  MatSetValue(matBPi, 
791
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
792 793 794
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
795 796 797 798 799
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
800 801
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
802 803

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
804
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
805
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
806 807

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
808
      MatDestroy(&matBPi);
809 810

      MatInfo minfo;
811
      MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
812 813
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

814
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
815
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
816
		      SAME_NONZERO_PATTERN);
817 818 819 820 821 822
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPPREONLY);
      PC pc_schur_primal;      
      KSPGetPC(ksp_schur_primal, &pc_schur_primal);
      PCSetType(pc_schur_primal, PCLU);
      PCFactorSetMatSolverPackage(pc_schur_primal, MATSOLVERMUMPS);
823
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
824

825 826
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
827
    }
828 829 830 831 832 833 834
  }


  void PetscSolverFeti::destroySchurPrimalKsp()
  {
    FUNCNAME("PetscSolverFeti::destroySchurPrimal()");

835
    if (schurPrimalSolver == 0) {
836
      schurPrimalData.subSolver = NULL;
837

838 839 840
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }