PetscSolverFeti.cc

/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors:
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 *
 ******************************************************************************/


#include "MatrixVector.h"
#include "parallel/PetscHelper.h"
#include "parallel/PetscSolverFeti.h"
#include "parallel/PetscSolverFetiDebug.h"
#include "parallel/PetscSolverFetiMonitor.h"
#include "parallel/PetscSolverFetiStructs.h"
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
#include "parallel/PetscSolverGlobalMatrix.h"
#include "io/VtkWriter.h"
#include "utility/PetscWrapper.h"

namespace AMDiS { namespace Parallel {

  using namespace std;

  PetscSolverFeti::PetscSolverFeti(string name)
    : PetscSolver(name),
      fetiSolverType(EXACT),
      dofMapSubDomain(FESPACE_WISE, true),
      primalDofMap(COMPONENT_WISE),
      dualDofMap(COMPONENT_WISE),
      interfaceDofMap(COMPONENT_WISE),
      localDofMap(COMPONENT_WISE),
      lagrangeMap(COMPONENT_WISE),
      interiorDofMap(COMPONENT_WISE),
      schurPrimalSolver(0),
      levelMode(1),
      subDomainIsLocal(true),
      subdomain(NULL),
      massMatrixSolver(NULL),
      printTimings(false),
      augmentedLagrange(false),
      nRankEdges(0),
      nOverallEdges(0),
      dirichletMode(0),
      stokesMode(false),
      pressureComponent(-1)
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get(name + "->left precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "no") {
      MSG("Create FETI-DP solver with no preconditioner!\n");
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
      fetiPreconditioner = FETI_LUMPED;
    } else {
      ERROR_EXIT("Preconditioner \"%s\" not available!\n",
		 preconditionerName.c_str());
    }

    preconditionerName = "";
    Parameters::get(name + "->right precon", preconditionerName);
    if (preconditionerName != "" && preconditionerName != "no") {
      ERROR_EXIT("FETI-DP does not support right preconditioning! (parameter \"%s->right precon\" has value \"%s\")\n",
		 name.c_str(), preconditionerName.c_str());
    }

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

    Parameters::get(name + "->feti->stokes mode", stokesMode);
    if (stokesMode) {
      Parameters::get(name + "->feti->pressure component", pressureComponent);
      TEST_EXIT(pressureComponent >= 0)
	("FETI-DP in Stokes mode, no pressure component defined!\n");
    }

    Parameters::get(name + "->feti->augmented lagrange", augmentedLagrange);

    Parameters::get(name + "->feti->symmetric", isSymmetric);

    {
      MSG("WARNING: CHECK THIS HERE BEFORE GOING INTO RUNNING MULTILEVEL FETI-DP!\n");
      Parameters::get("parallel->level mode", levelMode);
    }

    {
      int tmp = 0;
      Parameters::get(name + "->feti->inexact", tmp);
      if (tmp == 1)
	fetiSolverType = INEXACT;
      if (tmp == 2)
	fetiSolverType = INEXACT_REDUCED;
    }

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


  void PetscSolverFeti::init(vector<const FiniteElemSpace*> &fe0,
			     vector<const FiniteElemSpace*> &fe1,
			     bool createGlobalMapping)
  {
    FUNCNAME_DBG("PetscSolverFeti::init()");

    super::init(fe0, fe1, createGlobalMapping);

    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
    int nLevels = levelData.getNumberOfLevels();
    TEST_EXIT_DBG(nLevels >= 1)("nLevels < 1! Should not happen!\n");

    if (createGlobalMapping) {
      if (meshLevel + 1 < nLevels &&
	  levelData.getMpiComm(meshLevel + 1) != MPI::COMM_SELF) {
	dofMapSubDomain.init(componentSpaces, feSpaces);
	dofMapSubDomain.setMpiComm(levelData.getMpiComm(meshLevel + 1));
	dofMapSubDomain.setDofComms(meshDistributor->getDofComms(), meshLevel + 1);
	dofMapSubDomain.clear();
	meshDistributor->registerDofMap(dofMapSubDomain);
      }
    }
  }


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

#ifndef NDEBUG
    MSG("Init FETI-DP on mesh level %d\n", meshLevel);
#endif

    TEST_EXIT_DBG(meshLevel + 2 <=
		  meshDistributor->getMeshLevelData().getNumberOfLevels())
      ("Mesh hierarchy does not contain at least %d levels!\n", meshLevel + 2);

    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

    subDomainIsLocal = (levelData.getMpiComm(meshLevel + 1) == MPI::COMM_SELF);

    if (subdomain == NULL) {
      string subSolverInitStr = name + "->subsolver";
      string solverType = "petsc";
      Parameters::get(subSolverInitStr, solverType);
      solverType = "p_" + solverType;
      LinearSolverCreator *solverCreator =
	dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolverInterface>::getCreator(solverType, name));
      TEST_EXIT(solverCreator)
	("No valid solver type found in parameter \"%s\"\n",
	 name.c_str());
      solverCreator->setName(subSolverInitStr);
      subdomain = dynamic_cast<PetscSolver*>(solverCreator->create());
      subdomain->setSymmetric(isSymmetric);
      subdomain->setHandleDirichletRows(dirichletMode == 0);
      subdomain->setMeshDistributor(meshDistributor, meshLevel + 1);
      subdomain->init(componentSpaces, feSpaces);

      delete solverCreator;
    }

    primalDofMap.init(componentSpaces, feSpaces);
    dualDofMap.init(componentSpaces, feSpaces, false);
    localDofMap.init(componentSpaces, feSpaces, !subDomainIsLocal);
    lagrangeMap.init(componentSpaces, feSpaces);

    if (stokesMode)
      interfaceDofMap.init(componentSpaces, feSpaces);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      TEST_EXIT(levelMode == 1)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

      interiorDofMap.init(componentSpaces, feSpaces, false);
    }
  }


  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

    if (dirichletMode == 1) {
      int nComponents = mat.getNumRows();
      for (int component = 0; component < nComponents; component++) {
	DOFMatrix* dofMat = mat[component][component];
	if (!dofMat)
	  continue;

	dirichletRows[component] = dofMat->getDirichletRows();
      }
    }
  }


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

    double timeCounter = MPI::Wtime();

    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.clear();

    primalDofMap.setDofComms(meshDistributor->getDofComms(), meshLevel);
    lagrangeMap.setDofComms(meshDistributor->getDofComms(), meshLevel);

    primalDofMap.setMpiComm(levelData.getMpiComm(meshLevel));
    dualDofMap.setMpiComm(levelData.getMpiComm(meshLevel));
    lagrangeMap.setMpiComm(levelData.getMpiComm(meshLevel));
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel + 1));
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel + 1));

    localDofMap.setDofComms(meshDistributor->getDofComms(), meshLevel + 1);

    if (stokesMode) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComms(meshDistributor->getDofComms(), meshLevel);
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0));
    }

    int nComponents = componentSpaces.size();
    for (int component = 0; component < nComponents; component++) {
      createPrimals(component);
      createDuals(component);
      createInterfaceNodes(component);
      createIndexB(component);
    }

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();

    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.update();

    if (stokesMode)
      interfaceDofMap.update();

    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
    }

    lagrangeMap.update();


    // === ===

    if (subDomainIsLocal) {
      MSG("WARNING: MAKE GENERAL!\n");

      rStartInterior = 0;
      int localDofs = localDofMap.getOverallDofs();

      mpi::getDofNumbering(domainComm, localDofs,
			   rStartInterior, nGlobalOverallInterior);
    } else {
      MSG("WARNING: MAKE GENERAL!\n");

      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

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

      mpi::getDofNumbering(domainComm, 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);
    }


    for (int i = 0; i < static_cast<int>(componentSpaces.size()); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];

      MSG("FETI-DP data for %d-ith component (FE space %p):\n", i, feSpace);

      if (i == pressureComponent) {
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[i].nRankDofs,
	    interfaceDofMap[i].nOverallDofs);
      } else {
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n",
	    primalDofMap[i].nRankDofs,
	    primalDofMap[i].nLocalDofs,
	    primalDofMap[i].nOverallDofs);

	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
	    dualDofMap[i].nRankDofs,
	    dualDofMap[i].nOverallDofs);

	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
	    lagrangeMap[i].nRankDofs,
	    lagrangeMap[i].nOverallDofs);

	MSG("  nRankLocal = %d  nOverallLocal = %d\n",
	    localDofMap[i].nRankDofs,
	    localDofMap[i].nOverallDofs);
      }
    }

    subdomain->setDofMapping(&localDofMap);
    subdomain->setCoarseSpaceDofMapping(&primalDofMap);
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);

    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      timeCounter = MPI::Wtime() - timeCounter;
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n",
	  timeCounter);
    }


    bool writePrimals = false;
    Parameters::get("parallel->debug->write primals", writePrimals);
    if (writePrimals)
      PetscSolverFetiDebug::writePrimalFiles(*this);
  }


  void PetscSolverFeti::createPrimals(int component)
  {
    FUNCNAME("PetscSolverFeti::createPrimals()");

    if (component == pressureComponent)
      return;

    const FiniteElemSpace *feSpace = componentSpaces[component];

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

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

    DofIndexSet primals;
    for (DofContainerSet::iterator it = vertices.begin();
	 it != vertices.end(); ++it) {

      if (dirichletRows[component].count(**it))
	continue;

      if (meshLevel == 1 && not (*interiorMap)[component].isSet(**it))
	continue;

      if (subDomainIsLocal) {
	primals.insert(**it);
      } else {
	double e = 1e-8;
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

	if ((fabs(c[0]) < e && fabs(c[1] - 12.5) < e) ||
	    (fabs(c[0] - 25.0) < e && fabs(c[1] - 12.5) < e) ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1]) < e) ||
	    (fabs(c[0] - 12.5) < 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);
	} else {
	  MSG("OMMIT SOME PRIMAL!\n");
	}
      }
    }

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

    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) {
      if (dofMap[feSpace].isRankDof(*it)) {
	primalDofMap[component].insertRankDof(*it);
      } else {
  	primalDofMap[component].insertNonRankDof(*it);
      }
    }
  }


  void PetscSolverFeti::createDuals(int component)
  {
    FUNCNAME("PetscSolverFeti::createDuals()");

    if (component == pressureComponent)
      return;

    const FiniteElemSpace *feSpace = componentSpaces[component];

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

    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[component].count(**it))
	continue;

      if (isPrimal(component, **it))
	continue;

      if (meshLevel == 1 && not (*interiorMap)[component].isSet(**it))
	continue;

      if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(**it))
	dualDofMap[component].insertRankDof(**it);
    }
  }


  void PetscSolverFeti::createInterfaceNodes(int component)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

    if (component != pressureComponent)
      return;

    const FiniteElemSpace *feSpace = componentSpaces[component];
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[component].count(**it))
	continue;

      if (dofMap[feSpace].isRankDof(**it))
	interfaceDofMap[component].insertRankDof(**it);
      else
	interfaceDofMap[component].insertNonRankDof(**it);
    }
  }


  void PetscSolverFeti::createLagrange(int component)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

    if (component == pressureComponent)
      return;

    const FiniteElemSpace *feSpace = componentSpaces[component];
    Mesh* mesh = feSpace->getMesh();
    boundaryDofRanks[feSpace].clear();

    // Stores for all rank owned communication DOFs, if the counterpart is
    // a rank owned DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFs that fulfill this requirenment.
    map<int, std::set<DegreeOfFreedom> > subDomainRankDofs;

    if (not subDomainIsLocal) {
      StdMpi<vector<int> > stdMpi(domainComm);

      for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getRecvDofs(), feSpace);
	   !it.end(); it.nextRank()) {

	vector<int> dofs;
	dofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
	  if (dofMapSubDomain[feSpace].isRankDof(it.getDofIndex()))
	    dofs.push_back(1);
	  else
	    dofs.push_back(0);
	}

	stdMpi.send(it.getRank(), dofs);
      }

      for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getSendDofs(), feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getSendDofs(), feSpace);
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(component, it.getDofIndex()) &&
	      stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	    subDomainRankDofs[it.getRank()].insert(it.getDofIndex());
    }

    if (dualDofMap[component].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)).         ===

    int mpiRank = domainComm.Get_rank();
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getSendDofs(), feSpace);
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(component, it.getDofIndex())) {
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

	  // If the subdomain is local, always add the counterpart rank,
	  // otherwise check if the other rank is the owner of the DOF in
	  // its subdomain.
 	  if (subDomainIsLocal ||
	      subDomainRankDofs[it.getRank()].count(it.getDofIndex()))
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
	}
      }
    }


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

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

    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
	if (!isPrimal(component, it.getDofIndex()))
 	  if (subDomainIsLocal || subDomainRankDofs[it.getRank()].count(it.getDofIndex()))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);

    stdMpi.updateSendDataSize();

    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
      bool recvFromRank = false;
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(component, it.getDofIndex())) {
 	  if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(it.getDofIndex())) {
	    recvFromRank = true;
	    break;
	  }
	}
      }

      if (recvFromRank)
	stdMpi.recv(it.getRank());
    }

    stdMpi.startCommunication();

    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
      int i = 0;
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(component, it.getDofIndex())) {
 	  if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(it.getDofIndex())) {
	    boundaryDofRanks[feSpace][it.getDofIndex()] =
	      stdMpi.getRecvData(it.getRank())[i++];
	  } else {
	    lagrangeMap[component].insertNonRankDof(it.getDofIndex());
	  }
	}
      }
    }


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

    int nRankLagrange = 0;
    DofMap& dualMap = dualDofMap[component].getMap();
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (dofMap[feSpace].isRankDof(it->first)) {
	lagrangeMap[component].insertRankDof(it->first, nRankLagrange);
	int degree = boundaryDofRanks[feSpace][it->first].size();
	nRankLagrange += (degree * (degree - 1)) / 2;
      } else {
	lagrangeMap[component].insertNonRankDof(it->first);
      }
    }
    lagrangeMap[component].nRankDofs = nRankLagrange;
  }


  void PetscSolverFeti::createAugmentedLagrange(int component)
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


  void PetscSolverFeti::createIndexB(int component)
  {
    FUNCNAME("PetscSolverFeti::createIndexB()");

    const FiniteElemSpace *feSpace = componentSpaces[component];
    DOFAdmin* admin = feSpace->getAdmin();

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

    int nLocalInterior = 0;
    for (int i = 0; i < admin->getUsedSize(); i++) {
      if (admin->isDofFree(i) ||
	  isPrimal(component, i) ||
	  isDual(component, i) ||
	  isInterface(component, i) ||
	  dirichletRows[component].count(i))
	continue;

      if (meshLevel == 1 &&  not (*interiorMap)[component].isSet(i))
	continue;

      if (subDomainIsLocal) {
	localDofMap[component].insertRankDof(i, nLocalInterior);

	if (fetiPreconditioner == FETI_DIRICHLET)
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);

	nLocalInterior++;
      } else {
	if (dofMapSubDomain[feSpace].isRankDof(i))
	  localDofMap[component].insertRankDof(i);
	else
	  localDofMap[component].insertNonRankDof(i);

	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");
      }
    }

    // === And finally, add the global indicies of all dual nodes. ===

    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
	 it != dualDofMap[component].getMap().end(); ++it) {
      if (subDomainIsLocal) {
	localDofMap[component].insertRankDof(it->first);
      } else {
	if (dofMapSubDomain[feSpace].isRankDof(it->first))
	  localDofMap[component].insertRankDof(it->first);
	else
	  localDofMap[component].insertNonRankDof(it->first);
      }
    }
  }


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

    double wtime = MPI::Wtime();
    int mpiRank = domainComm.Get_rank();

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

    MatCreateAIJ(domainComm,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);


    Vec vec_scale_lagrange;
    createVec(lagrangeMap, vec_scale_lagrange);


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

    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
      DofMap &dualMap = dualDofMap[component].getMap();

      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
	TEST_EXIT_DBG(boundaryDofRanks[componentSpaces[component]].count(it->first))
	  ("Should not happen!\n");

	// Global index of the first Lagrange constriant for this node.
	int index = lagrangeMap.getMatIndex(component, it->first);

	// Copy set of all ranks that contain this dual node.
	vector<int> W(boundaryDofRanks[componentSpaces[component]][it->first].begin(),
		      boundaryDofRanks[componentSpaces[component]][it->first].end());
	// Number of ranks that contain this dual node.
	int degree = W.size();

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");

	int counter = 0;
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
	      MatSetValue(mat_lagrange,
			  index + counter,
			  localDofMap.getMatIndex(component, it->first) + rStartInterior,
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
	    }
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===

	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
	  int nConstraints = (degree * (degree - 1)) / 2;
	  for (int i = 0; i < nConstraints; i++) {
	    VecSetValue(vec_scale_lagrange,
			index + i,
			1.0 / static_cast<double>(degree),
			INSERT_VALUES);
	  }
	}
      }
    }

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


#ifndef NDEBUG