PetscSolverFeti.cc 48.5 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/SubDomainSolver.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 236 237 238 239 240 241 242 243
    if (subDomainSolver == NULL) {
      if (meshLevel == 0) {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);
      } else {
	subDomainSolver = 
	  new SubDomainSolver(meshDistributor, 
			      levelData.getMpiComm(meshLevel - 1),
			      levelData.getMpiComm(meshLevel));
244
	subDomainSolver->setLevel(meshLevel);
245 246
      }
    }
247

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

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

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

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

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


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

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

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

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

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

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

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

320 321
      createPrimals(feSpace);  

322
      createDuals(feSpace);      
323

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

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

Thomas Witkowski's avatar
Blub  
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);
    }
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
338
    
339 340 341
    lagrangeMap.update();


342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      MSG("RANK %d FROM %d\n",
	  levelData.getMpiComm(1).Get_rank(),
	  levelData.getMpiComm(1).Get_size());

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

      mpi::getDofNumbering(mpiComm, groupRowsInterior,
			   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);

      MSG("COMM TEST FETI-DP: %d %d %d %d %d\n",
	  levelData.getMpiComm(1).Get_size(),
	  localDofMap.getRankDofs(),
	  localDofMap.getOverallDofs(),
	  nGlobalOverallInterior, rStartInterior);
    }

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

379
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
380 381
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
382
      
383
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
384 385
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
386

387
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
388 389
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
390

391
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
392 393
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
394
    }
395 396 397


    // If multi level test, inform sub domain solver about coarse space.
398
    subDomainSolver->setDofMapping(&primalDofMap, &localDofMap);
399 400 401
  }


402
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
403
  {
404
    FUNCNAME("PetscSolverFeti::createPrimals()");  
405

406 407 408
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

409 410
    /// Set of DOF indices that are considered to be primal variables.
    
411
    DofContainerSet& vertices = 
412
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
413 414

    DofIndexSet primals;
415
    for (DofContainerSet::iterator it = vertices.begin(); 
416 417 418
	 it != vertices.end(); ++it) {
      WorldVector<double> c;
      feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
419

420
      double e = 1e-8;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
421
      if (meshLevel == 0) {
422
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
423 424 425 426 427 428 429 430
      } 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);
	}
      }
431
    }
432

433 434 435 436

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

437
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
438 439
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
440
      else
441
  	primalDofMap[feSpace].insertNonRankDof(*it);
442 443 444
  }


445
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
446 447
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
448

449 450 451
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
452
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
453 454 455 456

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
457 458 459 460 461 462
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
	    dualDofMap[feSpace].insertRankDof(**it);
	}	  
463 464 465 466 467 468 469
  }

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

470 471
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
472 473 474 475
    // 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;
476

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
477
    if (meshLevel > 0) {
478 479 480 481 482 483 484 485 486 487
      StdMpi<vector<int> > stdMpi(mpiComm);

      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
488
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
489
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
490
	  else
491 492 493 494 495 496 497 498 499 500 501 502 503 504 505
	    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
506 507 508 509 510 511
	   !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());
    }
512 513


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530
    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);

	  if (meshLevel == 0 ||
	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
531 532 533 534
	}
      }
    }

535 536 537 538

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

539
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
540

541
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
542 543
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
544
	if (!isPrimal(feSpace, it.getDofIndex()))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
545 546
	  if (meshLevel == 0 ||
	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
547
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
548 549 550

    stdMpi.updateSendDataSize();

551
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
552
	 !it.end(); it.nextRank()) {
553
      bool recvFromRank = false;
554
      for (; !it.endDofIter(); it.nextDof()) {
555
	if (!isPrimal(feSpace, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
556 557 558 559 560 561
	  if (meshLevel == 0 ||
	      (meshLevel > 0 && 
	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
	    recvFromRank = true;
	    break;
	  }
562
	}
563
      }
564 565

      if (recvFromRank)
566
	stdMpi.recv(it.getRank());
567
    }
568

569 570
    stdMpi.startCommunication();

571
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
572
	 !it.end(); it.nextRank()) {
573
      int i = 0;
574
      for (; !it.endDofIter(); it.nextDof())
575
	if (!isPrimal(feSpace, it.getDofIndex()))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
576 577 578 579 580
	  if (meshLevel == 0 ||
	      (meshLevel > 0 && 
	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
581 582
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
583

584 585 586
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

587
    int nRankLagrange = 0;
588
    DofMap& dualMap = dualDofMap[feSpace].getMap();
589
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
590 591
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
592
	int degree = boundaryDofRanks[feSpace][it->first].size();
593
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
594
      } else {
595
	lagrangeMap[feSpace].insertNonRankDof(it->first);
596 597
      }
    }
598
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
599 600 601
  }


602
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
603
  {
604
    FUNCNAME("PetscSolverFeti::createIndexB()");
605

606
    DOFAdmin* admin = feSpace->getAdmin();
607 608 609 610

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

612
    int nLocalInterior = 0;    
613
    for (int i = 0; i < admin->getUsedSize(); i++) {
614
      if (admin->isDofFree(i) == false && 
615
	  isPrimal(feSpace, i) == false &&
616
	  isDual(feSpace, i) == false) {
617 618
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
619

620 621 622 623 624 625 626 627 628
	  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);
629 630 631

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
632
	}
633
      }
634
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
635
    
636 637
    // === And finally, add the global indicies of all dual nodes. ===

638 639 640 641 642 643 644 645 646 647
    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);
      }
648 649 650
  }


651
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
652 653 654
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

655 656
    // === Create distributed matrix for Lagrange constraints. ===

657
    MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
658 659
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
660 661 662
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

663 664 665 666 667 668 669
    // === 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
670
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
671
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
672

673
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
674
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
675 676 677
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
678
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
679

Thomas Witkowski's avatar
Thomas Witkowski committed
680
	// Copy set of all ranks that contain this dual node.
681 682
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
683 684
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
685 686

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
687 688 689 690
	
	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
691 692
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
693

694
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
695
	      MSG("SET VALUE: %f\n", value);
696
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
697
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
698
	    index++;	      
699 700 701 702 703 704 705 706 707 708
	  }
	}
      }
    }

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


709
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
710
  {
711
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
712

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

716
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
717
      
718
      VecCreateMPI(mpiComm, 
719
		   localDofMap.getRankDofs(), 
720
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
721
		   &(schurPrimalData.tmp_vec_b));
722
      VecCreateMPI(mpiComm, 
723 724
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
725 726
		   &(schurPrimalData.tmp_vec_primal));

727
      MatCreateShell(mpiComm,
728 729 730 731
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
732 733 734 735 736
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
737
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
738
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
739 740
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
741 742
      KSPSetFromOptions(ksp_schur_primal);
    } else {
743 744
      MSG("Create direct schur primal solver!\n");

745 746
      double wtime = MPI::Wtime();

747 748 749
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

750 751 752
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
753

Thomas Witkowski's avatar
Thomas Witkowski committed
754
      Mat matBPi;
755
      MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
756
		      nRowsRankB, nRowsRankPrimal, 
757
		      nGlobalOverallInterior, nRowsOverallPrimal,
Thomas Witkowski's avatar
Thomas Witkowski committed
758 759
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
760

761 762 763
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
764

765 766
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

767
      for (int i = 0; i < nRowsRankB; i++) {
768
	PetscInt row = localDofMap.getStartDofs() + i;
769
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
770 771 772 773 774

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

775
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
776 777
      }

778
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
779
		primalDofMap.getLocalDofs())
780
	("Should not happen!\n");
781

782 783 784
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
785
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
786 787 788 789 790 791 792 793

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

794
	subDomainSolver->solve(tmpVec, tmpVec);
795

Thomas Witkowski's avatar
Thomas Witkowski committed
796
	PetscScalar *tmpValues;
797
	VecGetArray(tmpVec, &tmpValues);
798
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
799
	  MatSetValue(matBPi, 
800
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
801 802 803
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
804 805 806 807 808
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
809 810
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
811 812

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
813
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
814
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
815 816

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
817
      MatDestroy(&matBPi);
818 819

      MatInfo minfo;
820
      MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
821 822
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

823
      KSPCreate(mpiComm, &ksp_schur_primal);
824
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
825
		      SAME_NONZERO_PATTERN);
826 827 828 829 830 831
      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);
832
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
833

834 835
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
836
    }
837 838 839 840 841 842 843
  }


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

844
    if (schurPrimalSolver == 0) {
845
      schurPrimalData.subSolver = NULL;
846

847 848 849
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
850 851 852
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
853 854 855
  }


856
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
857 858 859
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

860 861
    // === Create FETI-DP solver object. ===

862
    fetiData.mat_lagrange = &mat_lagrange;
863
    fetiData.subSolver = subDomainSolver;
864
    fetiData.ksp_schur_primal = &ksp_schur_primal;