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


#include "AMDiS.h"
#include "MatrixVector.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"

namespace AMDiS {

  using namespace std;

  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
      schurPrimalSolver(0),
      multiLevelTest(false),
      subdomain(NULL),
      massMatrixSolver(NULL),
      meshLevel(0),
      rStartInterior(0),
      nGlobalOverallInterior(0),
      printTimings(false),
      stokesMode(false),
      augmentedLagrange(false),
      pressureComponent(-1)
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
      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());
    }

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

    Parameters::get("parallel->multi level test", multiLevelTest);
    if (multiLevelTest)
      meshLevel = 1;

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

    Parameters::get("parallel->feti->augmented lagrange", augmentedLagrange);
  }


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

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

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


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

      pressureFeSpace = feSpaces[pressureComponent];
    }
			   
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();

      if (meshLevel == 0) {
	subdomain->setMeshDistributor(meshDistributor, 
				      mpiCommGlobal, mpiCommLocal);
      } else {
	subdomain->setMeshDistributor(meshDistributor, 
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
	subdomain->setLevel(meshLevel);
      }
    }

    primalDofMap.init(levelData, feSpaces, uniqueFe);   
    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
    lagrangeMap.init(levelData, feSpaces, uniqueFe);

    if (stokesMode)
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);

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

      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
    }
  }


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

    bool removeDirichletRows = false;
    Parameters::get("parallel->feti->remove dirichlet", removeDirichletRows);
    if (!removeDirichletRows)
      return;

    int nComponents = mat.getSize();
    for (int i = 0; i < nComponents; i++) {
      DOFMatrix* dofMat = mat[i][i];
      if (!dofMat)
	continue;
      
      const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
      std::set<DegreeOfFreedom>& dRows = dofMat->getDirichletRows();
      if (dirichletRows.count(feSpace)) {
	WARNING("Implement test that for all components dirichlet rows are equal!\n");
      } else {
	dirichletRows[feSpace] = dRows;
      }
    }
  }


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

    double timeCounter = MPI::Wtime();

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

    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

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

    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());

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

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

    if (stokesMode) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

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

      createDuals(feSpace);

      createInterfaceNodes(feSpace);

      createIndexB(feSpace);     
    }

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

    if (stokesMode)
      interfaceDofMap.update();

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

    lagrangeMap.update();


    // === ===

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

      mpi::getDofNumbering(mpiCommGlobal, 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("FETI-DP data created on mesh level %d\n", meshLevel);
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      MSG("FETI-DP data for %d-ith FE space: %p\n", i, feSpace);

      if (feSpace == pressureFeSpace) {
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
	    primalDofMap[feSpace].nRankDofs, 
	    primalDofMap[feSpace].nLocalDofs,
	    primalDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
	    dualDofMap[feSpace].nRankDofs, 
	    dualDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
	    lagrangeMap[feSpace].nRankDofs, 
	    lagrangeMap[feSpace].nOverallDofs);

// 	TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
// 		      static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
// 	  ("Should not happen!\n");	
      }
    }

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


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

    if (feSpace == pressureFeSpace)
      return;

    // === 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[feSpace].count(**it))
	continue;

      if (meshLevel == 0) {
	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]) < 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);
	}
      }
    }


    // === 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 (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
      else
  	primalDofMap[feSpace].insertNonRankDof(*it);
  }


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

    if (feSpace == pressureFeSpace)
      return;

    // === 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[feSpace].count(**it))
	continue;

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

      if (meshLevel == 0) {
	dualDofMap[feSpace].insertRankDof(**it);
      } else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
	  dualDofMap[feSpace].insertRankDof(**it);
      }	  
    }
  }

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

    if (feSpace != pressureFeSpace)
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
	interfaceDofMap[feSpace].insertRankDof(**it);
      else
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
  }


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

    if (feSpace == pressureFeSpace)
      return;

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

    if (meshLevel > 0) {
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);

      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()) {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
	    subdomainRankDofs.push_back(1);
	  else
	    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); 
	   !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());
    }

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

    int mpiRank = meshDistributor->getMpiRank();
    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());	  
	}
      }
    }


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

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

    stdMpi.updateSendDataSize();

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

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

    stdMpi.startCommunication();

    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      int i = 0;
      for (; !it.endDofIter(); it.nextDof())
	if (!isPrimal(feSpace, it.getDofIndex()))
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
	  else
	    lagrangeMap[feSpace].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[feSpace].getMap();
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
	int degree = boundaryDofRanks[feSpace][it->first].size();
	nRankLagrange += (degree * (degree - 1)) / 2;
      } else {
	lagrangeMap[feSpace].insertNonRankDof(it->first);
      }
    }
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
  }


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

    if (!augmentedLagrange)
      return;
  }


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

    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(feSpace, i) ||
	  isDual(feSpace, i) ||
	  isInterface(feSpace, i) ||
	  dirichletRows[feSpace].count(i))
	continue;      

      if (meshLevel == 0) {
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
	
	if (fetiPreconditioner == FETI_DIRICHLET)
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	
	nLocalInterior++;	
      } else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	  localDofMap[feSpace].insertRankDof(i);
	else
	  localDofMap[feSpace].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[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);
      }
    }
  }


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

    double wtime = MPI::Wtime();
    int mpiRank = meshDistributor->getMpiRank();

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

    MatCreateAIJ(mpiCommGlobal,
		 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;
    lagrangeMap.createVec(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 k = 0; k < feSpaces.size(); k++) {
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();

      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
	int index = lagrangeMap.getMatIndex(k, it->first);

	// Copy set of all ranks that contain this dual node.
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][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(k, 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 (meshDistributor->getDofMap()[feSpaces[k]].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);


    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

    VecView(vec_scale_lagrange, PETSC_VIEWER_STDOUT_WORLD);

    if (fetiPreconditioner != FETI_NONE || stokesMode) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatDiagonalScale(mat_lagrange_scaled, vec_scale_lagrange, PETSC_NULL);
    }

    VecDestroy(&vec_scale_lagrange);


    // === Print final timings. ===

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


  void PetscSolverFeti::createMatAugmentedLagrange(vector<const FiniteElemSpace*> &feSpaces)
  {
    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

    nRankEdges = 0;
    nOverallEdges = 0;
    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
      if (it->rankObj.subObj == EDGE)
	nRankEdges++;
    
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

    MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);

    nRankEdges *= feSpaces.size();
    rStartEdges *= feSpaces.size();
    nOverallEdges *= feSpaces.size();

    MatCreateAIJ(mpiCommGlobal,
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
		 1, PETSC_NULL, 1, PETSC_NULL, 
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

    int rowCounter = rStartEdges;
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
      if (it->rankObj.subObj == EDGE) {
	for (int i = 0; i < feSpaces.size(); i++) {
	  DofContainer edgeDofs;
	  it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj, edgeDofs);

	  MSG("SIZE = %d\n", edgeDofs.size());

	  for (DofContainer::iterator it = edgeDofs.begin();
	       it != edgeDofs.end(); it++) {
	    TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false)
	      ("Should not be primal!\n");

	    int col = lagrangeMap.getMatIndex(i, **it);
	    double value = 1.0;
	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
	  }

	  rowCounter++;
	}
      }
    }

    MatAssemblyBegin(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);

    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 05a: %.5f seconds (creation of augmented lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
    }
  }


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

    if (schurPrimalSolver == 0) {
      MSG("Create iterative schur primal solver!\n");

      if (augmentedLagrange == false) {
	schurPrimalData.subSolver = subdomain;
	
	localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
	
	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getOverallDofs(), 
		       primalDofMap.getOverallDofs(),
		       &schurPrimalData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimal);	
      } else {
	schurPrimalAugmentedData.subSolver = subdomain;

	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b0, nGlobalOverallInterior);
	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b1, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalAugmentedData.tmp_vec_primal);
	lagrangeMap.createVec(schurPrimalAugmentedData.tmp_vec_lagrange);

	schurPrimalAugmentedData.mat_lagrange = &mat_lagrange;
	schurPrimalAugmentedData.mat_augmented_lagrange = &mat_augmented_lagrange;

	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges,
		       &schurPrimalAugmentedData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimalAugmented);
      }

      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
      KSPSetFromOptions(ksp_schur_primal);
    } else {
      MSG("Create direct schur primal solver!\n");

      TEST_EXIT(!augmentedLagrange)("Not yet supported!\n");

      double wtime = MPI::Wtime();

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

      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();

      Mat matBPi;
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);


      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;

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

      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);

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

	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
      }

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

      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);

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

	subdomain->solve(tmpVec, tmpVec);

	PetscScalar *tmpValues;
	VecGetArray(tmpVec, &tmpValues);
	for (int i  = 0; i < nRowsRankB; i++)
	  MatSetValue(matBPi, 
		      localDofMap.getStartDofs() + i,
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);

      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);

      Mat matPrimal;
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
      MatDestroy(&matBPi);

      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
		      SAME_NONZERO_PATTERN);
      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);
      KSPSetFromOptions(ksp_schur_primal);

      if (printTimings) {
	MPI::COMM_WORLD.Barrier();
	MatInfo minfo;
	MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
	MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);
	
	MSG("FETI-DP timing 06: %.5f seconds (creation of schur primal matrix)\n",
	    MPI::Wtime() - wtime);

	wtime = MPI::Wtime();
	KSPSetUp(ksp_schur_primal);
	KSPSetUpOnBlocks(ksp_schur_primal);
	MPI::COMM_WORLD.Barrier();
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
    }
  }


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

    if (schurPrimalSolver == 0) {
      if (augmentedLagrange == false) {
	schurPrimalData.subSolver = NULL;
	
	VecDestroy(&schurPrimalData.tmp_vec_b);
	VecDestroy(&schurPrimalData.tmp_vec_primal);
      } else {
	schurPrimalAugmentedData.subSolver = NULL;
	schurPrimalAugmentedData.mat_lagrange = NULL;
	schurPrimalAugmentedData.mat_augmented_lagrange = NULL;

	VecDestroy(&schurPrimalAugmentedData.tmp_vec_b0);
	VecDestroy(&schurPrimalAugmentedData.tmp_vec_b1);
	VecDestroy(&schurPrimalAugmentedData.tmp_vec_primal);
	VecDestroy(&schurPrimalAugmentedData.tmp_vec_lagrange);
      }
    }
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
  }


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

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

    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
    primalDofMap.createVec(fetiData.tmp_vec_primal0);

    if (stokesMode == false) {      
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      if (augmentedLagrange == false) {
	MatShellSetOperation(mat_feti, MATOP_MULT, 
			     (void(*)(void))petscMultMatFeti);
      } else {
	fetiData.mat_augmented_lagrange = &mat_augmented_lagrange;
	primalDofMap.createVec(fetiData.tmp_vec_primal1);
	MatShellSetOperation(mat_feti, MATOP_MULT, 
			     (void(*)(void))petscMultMatFetiAugmented);	
      }
    }  else {
      TEST_EXIT_DBG(!augmentedLagrange)("Not yet implemented!\n");

      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFetiInterface);      
    }

    KSPCreate(mpiCommGlobal, &ksp_feti);
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
    KSPSetFromOptions(ksp_feti);
    if (stokesMode)
      KSPMonitorSet(ksp_feti, KSPMonitorFetiStokes, &fetiKspData, PETSC_NULL);

    // === Create null space objects. ===
    
    createNullSpace();


    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:
      TEST_EXIT(meshLevel == 0)
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

      TEST_EXIT(!stokesMode)
	("Stokes mode does not yet support the Dirichlet precondition!\n");

      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
      KSPSetType(ksp_interior, KSPPREONLY);
      PC pc_interior;
      KSPGetPC(ksp_interior, &pc_interior);
      PCSetType(pc_interior, PCLU);
      PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK);
      KSPSetFromOptions(ksp_interior);
            
      fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
      fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
      fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
      fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
      fetiDirichletPreconData.ksp_interior = &ksp_interior;
      
      localDofMap.createVec(fetiDirichletPreconData.tmp_vec_b, 
			    nGlobalOverallInterior);
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_interior));

      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;	  
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);


      // For the case, that we want to print the timings, we force the LU 
      // factorization of the local problems to be done here explicitly.
      if (printTimings) {
	double wtime = MPI::Wtime();
	KSPSetUp(ksp_interior);
	KSPSetUpOnBlocks(ksp_interior);
	MPI::COMM_WORLD.Barrier();
	MSG("FETI-DP timing 08: %.5f seconds (factorization of Dirichlet preconditoner matrices)\n",
	    MPI::Wtime() - wtime);
      }
      
      break;

    case FETI_LUMPED:
      {
	FetiLumpedPreconData *lumpedData = 
	  (stokesMode ? &fetiInterfaceLumpedPreconData : &fetiLumpedPreconData);

	lumpedData->mat_lagrange_scaled = &mat_lagrange_scaled;
	lumpedData->mat_duals_duals = &mat_duals_duals;
	
	localDofMap.createVec(lumpedData->tmp_vec_b0);
	MatGetVecs(mat_duals_duals, PETSC_NULL, 
		   &(lumpedData->tmp_vec_duals0));
	MatGetVecs(mat_duals_duals, PETSC_NULL, 
		   &(lumpedData->tmp_vec_duals1));
	
	for (unsigned int i = 0; i < feSpaces.size(); i++) {
	  if (stokesMode && i == pressureComponent)
	    continue;
	  
	  DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
	  for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
	    DegreeOfFreedom d = it->first;
	    int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	    int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	    lumpedData->localToDualMap[matIndexLocal] = matIndexDual;
	  }
	}

	if (stokesMode) {
	  // === Create H2 vec ===

	  DOFVector<double> tmpvec(pressureFeSpace, "tmpvec");
	  createH2vec(tmpvec);
	  interfaceDofMap.createVec(fetiInterfaceLumpedPreconData.h2vec);

	  DofMap& m = interfaceDofMap[pressureFeSpace].getMap();
	  for (DofMap::iterator it = m.begin(); it != m.end(); ++it) {
	    if (meshDistributor->getDofMap()[pressureFeSpace].isRankDof(it->first)) {
	      int index = interfaceDofMap.getMatIndex(pressureComponent, it->first);
	      VecSetValue(fetiInterfaceLumpedPreconData.h2vec, 
			  index, tmpvec[it->first], INSERT_VALUES);
	    }
	  }  

	  VecAssemblyBegin(fetiInterfaceLumpedPreconData.h2vec);
	  VecAssemblyEnd(fetiInterfaceLumpedPreconData.h2vec);

	  // === Create mass matrix solver ===
	  
	  if (!massMatrixSolver) {
	    MSG("START CREATE MASS MATRIX!\n");
	    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
	    Operator op(pressureFeSpace, pressureFeSpace);
	    Simple_ZOT zot;
	    op.addTerm(&zot);
	    massMatrix.assembleOperator(op);

	    ParallelDofMapping massMatMapping = interfaceDofMap;
	    massMatrixSolver = new PetscSolverGlobalMatrix;
	    massMatrixSolver->setKspPrefix("mass_");
	    massMatrixSolver->setMeshDistributor(meshDistributor,
						 mpiCommGlobal,
						 mpiCommLocal);
	    massMatrixSolver->setDofMapping(&massMatMapping);
	    massMatrixSolver->fillPetscMatrix(&massMatrix);

	    int r, c;
	    MatGetSize(massMatrixSolver->getMatInterior(), &r, &c);
	    MatInfo info;
	    MatGetInfo(massMatrixSolver->getMatInterior(), MAT_GLOBAL_SUM, &info);
	    MSG("MASS MAT INFO:  size = %d x %d   nnz = %d\n",
		r, c,
		static_cast<int>(info.nz_used));
	    MSG("END CREATE MASS MATRIX!\n");

	    fetiInterfaceLumpedPreconData.ksp_mass = 
	      massMatrixSolver->getSolver();
	  }

	  // === Create tmp vectors ===
 
	  localDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_vec_b1);
	  primalDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_primal);
	  fetiInterfaceLumpedPreconData.subSolver = subdomain;	
	}
	
	KSPGetPC(ksp_feti, &precon_feti);
	PCSetType(precon_feti, PCSHELL);

	if (stokesMode) {
	  PCShellSetContext(precon_feti, static_cast<void*>(&fetiInterfaceLumpedPreconData));
	  PCShellSetApply(precon_feti, petscApplyFetiInterfaceLumpedPrecon);
	} else {
	  PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
	  PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
	}
      }

      break;
    default:
      break;
    }
  }
  

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

    // === Destroy FETI-DP solver object. ===

    fetiData.mat_lagrange = PETSC_NULL;
    fetiData.subSolver = NULL;
    fetiData.ksp_schur_primal = PETSC_NULL;

    VecDestroy(&fetiData.tmp_vec_b0);
    VecDestroy(&fetiData.tmp_vec_lagrange);
    VecDestroy(&fetiData.tmp_vec_primal0);

    if (augmentedLagrange) {
      fetiData.mat_augmented_lagrange = PETSC_NULL;
      VecDestroy(&fetiData.tmp_vec_primal1);
    }

    if (stokesMode) {
      VecDestroy(&fetiData.tmp_vec_b1);
      VecDestroy(&fetiData.tmp_vec_primal1);
      VecDestroy(&fetiData.tmp_vec_interface);
    }

    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);


    // === Destroy FETI-DP preconditioner object. ===

    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPDestroy(&ksp_interior);

      fetiDirichletPreconData.mat_lagrange_scaled = NULL;
      fetiDirichletPreconData.mat_interior_interior = NULL;
      fetiDirichletPreconData.mat_duals_duals = NULL;
      fetiDirichletPreconData.mat_interior_duals = NULL;
      fetiDirichletPreconData.mat_duals_interior = NULL;
      fetiDirichletPreconData.ksp_interior = NULL;
      
      VecDestroy(&fetiDirichletPreconData.tmp_vec_b);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals0);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals1);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_interior);
      MatDestroy(&mat_lagrange_scaled);
      break;

    case FETI_LUMPED:
      {
	FetiLumpedPreconData &lumpedData = 
	  (stokesMode ? fetiInterfaceLumpedPreconData : fetiLumpedPreconData);
	
	lumpedData.mat_lagrange_scaled = NULL;
	lumpedData.mat_duals_duals = NULL;
	
	VecDestroy(&lumpedData.tmp_vec_b0);
	VecDestroy(&lumpedData.tmp_vec_duals0);
	VecDestroy(&lumpedData.tmp_vec_duals1);
      }
      break;
    default:
      break;
    }
  }


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

    if (!stokesMode)
      return;

    Vec ktest0, ktest1;
    localDofMap.createLocalVec(ktest0);
    localDofMap.createLocalVec(ktest1);
    DofMap& m = localDofMap[pressureFeSpace].getMap();
    for (DofMap::iterator it = m.begin(); it != m.end(); ++it) {
      if (meshDistributor->getDofMap()[pressureFeSpace].isRankDof(it->first)) {
	int index = localDofMap.getLocalMatIndex(pressureComponent, it->first);
	VecSetValue(ktest0, index, 1.0, INSERT_VALUES);
      }
    }    
    VecAssemblyBegin(ktest0);
    VecAssemblyEnd(ktest0);
    MatMult(subdomain->getMatInterior(), ktest0, ktest1);
    
    PetscScalar *valarray;
    Vec ktest2, ktest3;
    VecGetArray(ktest1, &valarray);
    VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, 
			  localDofMap.getRankDofs(), nGlobalOverallInterior, 
			  valarray, &ktest2);
    localDofMap.createVec(ktest3, nGlobalOverallInterior);
    
    Vec vecArray[2];
    interfaceDofMap.createVec(vecArray[0]);
    lagrangeMap.createVec(vecArray[1]);

    VecSet(vecArray[0], 1.0);   
    MatMult(subdomain->getMatInteriorCoarse(1), vecArray[0], ktest3);
    VecAXPY(ktest3, 1.0, ktest2);
    MatMult(mat_lagrange_scaled, ktest3, vecArray[1]);
    VecScale(vecArray[1], -1.0);
        
    Vec nullSpaceBasis;
    VecCreateNest(mpiCommGlobal, 2, PETSC_NULL, vecArray, &nullSpaceBasis); 

#if (DEBUG != 0)
    PetscSolverFetiDebug::writeNullSpace(*this, nullSpaceBasis);
#endif

    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, 
		       &matNullSpace);
    MatSetNullSpace(mat_feti, matNullSpace);
    KSPSetNullSpace(ksp_feti, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);

    VecDestroy(&ktest0);
    VecDestroy(&ktest1);
    VecDestroy(&ktest2);
    VecDestroy(&ktest3);

    VecDestroy(&(vecArray[0]));
    VecDestroy(&(vecArray[1]));
    VecDestroy(&nullSpaceBasis);
  }


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

#if (DEBUG != 0)
    PetscInt nRow, nCol;
    MatGetSize(subdomain->getMatInterior(), &nRow, &nCol);
    mpi::globalAdd(nRow);
    mpi::globalAdd(nCol);

    MatInfo minfo;
    MatGetInfo(subdomain->getMatInterior(), MAT_GLOBAL_SUM, &minfo);
    int nnz = static_cast<int>(minfo.nz_used);
    mpi::globalAdd(nnz);
    MSG("Interior matrices [%d x %d] nnz = %d\n", nRow, nCol, nnz);

    MatGetSize(subdomain->getMatCoarse(), &nRow, &nCol);
    MatGetInfo(subdomain->getMatCoarse(), MAT_GLOBAL_SUM, &minfo);
    MSG("Primal matrix [%d x %d] nnz = %d\n", nRow, nCol, 
	static_cast<int>(minfo.nz_used));

    MatGetSize(subdomain->getMatCoarseInterior(), &nRow, &nCol);
    MatGetInfo(subdomain->getMatCoarseInterior(), MAT_GLOBAL_SUM, &minfo);
    MSG("Primal-Interior matrix [%d x %d] nnz = %d\n", nRow, nCol, 
	static_cast<int>(minfo.nz_used));

    MatGetSize(subdomain->getMatInteriorCoarse(), &nRow, &nCol);    
    MatGetInfo(subdomain->getMatInteriorCoarse(), MAT_GLOBAL_SUM, &minfo);
    MSG("Interior-Primal matrix [%d x %d] nnz = %d\n", nRow, nCol, 
	static_cast<int>(minfo.nz_used));
    
    if (stokesMode) {
      MatGetSize(subdomain->getMatCoarse(1, 1), &nRow, &nCol);   
      MatGetInfo(subdomain->getMatCoarse(1, 1), MAT_GLOBAL_SUM, &minfo);
      MSG("Interface matrix [%d x %d] nnz = %d\n", nRow, nCol, 
	  static_cast<int>(minfo.nz_used));

      MatGetSize(subdomain->getMatCoarseInterior(1), &nRow, &nCol);         
      MatGetInfo(subdomain->getMatCoarseInterior(1), MAT_GLOBAL_SUM, &minfo);
      MSG("Interface-Interior matrix [%d x %d] nnz = %d\n", nRow, nCol, 
	  static_cast<int>(minfo.nz_used));
      
      MatGetSize(subdomain->getMatInteriorCoarse(1), &nRow, &nCol);   
      MatGetInfo(subdomain->getMatInteriorCoarse(1), MAT_GLOBAL_SUM, &minfo);
      MSG("Interior-Interface matrix [%d x %d] nnz = %d\n", nRow, nCol, 
	  static_cast<int>(minfo.nz_used));

      MatGetSize(subdomain->getMatCoarse(1, 0), &nRow, &nCol);   
      MatGetInfo(subdomain->getMatCoarse(1, 0), MAT_GLOBAL_SUM, &minfo);
      MSG("Interface-Primal matrix [%d x %d] nnz = %d\n", nRow, nCol, 
	  static_cast<int>(minfo.nz_used));

      MatGetSize(subdomain->getMatCoarse(0, 1), &nRow, &nCol);
      MatGetInfo(subdomain->getMatCoarse(0, 1), MAT_GLOBAL_SUM, &minfo);
      MSG("Primal-Interface matrix [%d x %d] nnz = %d\n", nRow, nCol, 
	  static_cast<int>(minfo.nz_used));
    }
#endif

    int writeInteriorMatrix = -1;
    Parameters::get("parallel->debug->write interior matrix", 
		    writeInteriorMatrix);

    if (writeInteriorMatrix >= 0 &&
	writeInteriorMatrix == MPI::COMM_WORLD.Get_rank()) {
      PetscViewer petscView;
      PetscViewerBinaryOpen(PETSC_COMM_SELF, "interior.mat", 
			    FILE_MODE_WRITE, &petscView);
      MatView(subdomain->getMatInterior(), petscView);
      PetscViewerDestroy(&petscView);
    }


    int writeCoarseMatrix = 0;
    Parameters::get("parallel->debug->write coarse matrix", 
		    writeCoarseMatrix);
    
    if (writeCoarseMatrix > 0) {
      PetscViewer petscView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "coarse.mat", 
			    FILE_MODE_WRITE, &petscView);
      MatView(subdomain->getMatCoarse(), petscView);
      PetscViewerDestroy(&petscView);
    }

    
    int writeSchurPrimalMatrix = 0;
    Parameters::get("parallel->debug->write schur primal matrix",
		    writeSchurPrimalMatrix);
    if (writeSchurPrimalMatrix) {
      PetscViewer petscView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "schurprimal.mat", 
			    FILE_MODE_WRITE, &petscView);
      MatView(mat_schur_primal, petscView);
      PetscViewerDestroy(&petscView);
    }
  }


  void PetscSolverFeti::recoverSolution(Vec &vec_sol_b,
					Vec &vec_sol_primal,
					SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::recoverSolution()");

    // === Get local part of the solution for B variables. ===

    PetscScalar *localSolB;
    VecGetArray(vec_sol_b, &localSolB);

    PetscInt bsize;
    VecGetLocalSize(vec_sol_b, &bsize);

    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
    
    vector<const FiniteElemSpace*> feSpaces = vec.getComponentFeSpaces();

    vector<PetscInt> globalIsIndex, localIsIndex;
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());

    {
      int cnt = 0;
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	DofMap& dofMap = primalDofMap[feSpaces[i]].getMap();
	for (DofMap::iterator it = dofMap.begin();it != dofMap.end(); ++it) {
	  globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
	  localIsIndex.push_back(cnt++);
	}
      }
      
      TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs())
	("Should not happen!\n");
    }
    
    IS globalIs, localIs;
    ISCreateGeneral(PETSC_COMM_SELF, 
		    globalIsIndex.size(), 
		    &(globalIsIndex[0]),
		    PETSC_USE_POINTER,
		    &globalIs);

    ISCreateGeneral(PETSC_COMM_SELF, 
		    localIsIndex.size(), 
		    &(localIsIndex[0]),
		    PETSC_USE_POINTER,
		    &localIs);

    Vec local_sol_primal;
    VecCreateSeq(PETSC_COMM_SELF, localIsIndex.size(), &local_sol_primal);

    VecScatter primalScatter;
    VecScatterCreate(vec_sol_primal, globalIs, local_sol_primal, localIs, &primalScatter);
    VecScatterBegin(primalScatter, vec_sol_primal, local_sol_primal, 
		    INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(primalScatter, vec_sol_primal, local_sol_primal, 
		  INSERT_VALUES, SCATTER_FORWARD);

    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    

    PetscScalar *localSolPrimal;
    VecGetArray(local_sol_primal, &localSolPrimal);

    // === And copy from PETSc local vectors to the DOF vectors. ===
    
    int cnt = 0;
    for (int i = 0; i < vec.getSize(); i++) {
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));
      
      for (DofMap::iterator it = localDofMap[feSpaces[i]].getMap().begin();
	   it != localDofMap[feSpaces[i]].getMap().end(); ++it) {
	if (meshLevel == 0) {
	  int petscIndex = localDofMap.getLocalMatIndex(i, it->first);
	  dofVec[it->first] = localSolB[petscIndex];
	} else {
	  if (meshDistributor->getDofMapSd()[feSpaces[i]].isRankDof(it->first)) {
	    int petscIndex = localDofMap.getLocalMatIndex(i, it->first);
	    TEST_EXIT(petscIndex < bsize)("Should not happen!\n");
	    dofVec[it->first] = localSolB[petscIndex];	 
	  }
	}
      }

      for (DofMap::iterator it = primalDofMap[feSpaces[i]].getMap().begin();
	   it != primalDofMap[feSpaces[i]].getMap().end(); ++it)
	dofVec[it->first] = localSolPrimal[cnt++];      
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
    VecDestroy(&local_sol_primal);
  }    


  void PetscSolverFeti::recoverInterfaceSolution(Vec& vecInterface, SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::recoverInterfaceSolution()");

    if (!stokesMode)
      return;

    vector<PetscInt> globalIsIndex, localIsIndex;
    globalIsIndex.reserve(interfaceDofMap.getLocalDofs());
    localIsIndex.reserve(interfaceDofMap.getLocalDofs());
    
    const FiniteElemSpace *pressureFeSpace = 
      vec.getDOFVector(pressureComponent)->getFeSpace();
    int cnt = 0;
    DofMap& dofMap = interfaceDofMap[pressureFeSpace].getMap();
    for (DofMap::iterator it = dofMap.begin();it != dofMap.end(); ++it) {
      globalIsIndex.push_back(interfaceDofMap.getMatIndex(pressureComponent, 
							  it->first));
      localIsIndex.push_back(cnt++);	
    }
    
    IS globalIs, localIs;
    ISCreateGeneral(PETSC_COMM_SELF, 
		    globalIsIndex.size(), 
		    &(globalIsIndex[0]),
		    PETSC_USE_POINTER,
		    &globalIs);
    
    ISCreateGeneral(PETSC_COMM_SELF, 
		    localIsIndex.size(), 
		    &(localIsIndex[0]),
		    PETSC_USE_POINTER,
		    &localIs);
    
    Vec local_sol_interface;
    VecCreateSeq(PETSC_COMM_SELF, localIsIndex.size(), &local_sol_interface);
    
    VecScatter interfaceScatter;
    VecScatterCreate(vecInterface, globalIs, local_sol_interface, localIs, &interfaceScatter);
    VecScatterBegin(interfaceScatter, vecInterface, local_sol_interface, 
		    INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(interfaceScatter, vecInterface, local_sol_interface, 
		  INSERT_VALUES, SCATTER_FORWARD);
    
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&interfaceScatter);    
    
    
    PetscScalar *localSolInterface;
    VecGetArray(local_sol_interface, &localSolInterface);
    
    // === And copy from PETSc local vectors to the DOF vectors. ===
    
    cnt = 0;
    DOFVector<double>& dofVec = *(vec.getDOFVector(pressureComponent));
    for (DofMap::iterator it = interfaceDofMap[pressureFeSpace].getMap().begin();
	 it != interfaceDofMap[pressureFeSpace].getMap().end(); ++it) {
      dofVec[it->first] = localSolInterface[cnt++];
    }
    
    VecRestoreArray(local_sol_interface, &localSolInterface);
    VecDestroy(&local_sol_interface);  
  }


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

    // === Create all sets and indices. ===
    
    vector<const FiniteElemSpace*> feSpaces = AMDiS::getComponentFeSpaces(*mat);

    initialize(feSpaces);

    createDirichletData(*mat);

    createFetiData();

    dirichletRows.clear();

   
    // === Create matrices for the FETI-DP method. ===

    if (printTimings)
      MPI::COMM_WORLD.Barrier();
    double wtime = MPI::Wtime();
  
    subdomain->fillPetscMatrix(mat);

    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 02: %.5f seconds (creation of interior matrices)\n",
	  MPI::Wtime() - wtime);

      // For the case, that we want to print the timings, we force the LU 
      // factorization of the local problems to be done here explicitly.
      wtime = MPI::Wtime();
      KSPSetUp(subdomain->getSolver());
      KSPSetUpOnBlocks(subdomain->getSolver());
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 04: %.5f seconds (factorization of subdomain matrices)\n",
	  MPI::Wtime() - wtime);
    }
   
    
    // === Create matrices for FETI-DP preconditioner. ===    
    createPreconditionerMatrix(mat);
   
    // === Create and fill PETSc matrix for Lagrange constraints. ===
    createMatLagrange(feSpaces);

    // === ===
    createMatAugmentedLagrange(feSpaces);

    // === Create PETSc solver for the Schur complement on primal variables. ===   
    createSchurPrimalKsp(feSpaces);

    // === Create PETSc solver for the FETI-DP operator. ===
    createFetiKsp(feSpaces);

    // === If required, run debug tests. ===
    dbgMatrix();
  }


  void PetscSolverFeti::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverFeti::fillPetscRhs()");

    subdomain->fillPetscRhs(vec);
  }


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

    if (fetiPreconditioner == FETI_NONE)
      return;

    double wtime = MPI::Wtime();
    vector<const FiniteElemSpace*> feSpaces = AMDiS::getComponentFeSpaces(*mat);
    int nRowsDual = dualDofMap.getRankDofs();
    
    MatCreateSeqAIJ(PETSC_COMM_SELF,
		    nRowsDual, nRowsDual, 100, PETSC_NULL,
		    &mat_duals_duals);
    MatSetOption(mat_duals_duals, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
    
    if (fetiPreconditioner == FETI_DIRICHLET) {
      int nRowsInterior = interiorDofMap.getRankDofs();
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsInterior, nRowsInterior, 100, PETSC_NULL,
		      &mat_interior_interior);
      MatSetOption(mat_interior_interior, 
		   MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsInterior, nRowsDual, 100, PETSC_NULL,
		      &mat_interior_duals);
      MatSetOption(mat_interior_duals, 
		   MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsDual, nRowsInterior, 100, PETSC_NULL,
		      &mat_duals_interior);
      MatSetOption(mat_duals_interior, 
		   MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
    }
    
    // === Prepare traverse of sequentially created matrices. ===
    
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;
    
    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;
    
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);
    
    
    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===
    
    int nComponents = mat->getSize();
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*mat)[rowComponent][colComponent];

	if (!dofMat)
	  continue;
	
	if (stokesMode && 
	    (rowComponent == pressureComponent || 
	     colComponent == pressureComponent))
	  continue;

	const FiniteElemSpace *rowFeSpace = feSpaces[rowComponent];
	const FiniteElemSpace *colFeSpace = feSpaces[colComponent];

	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
	
	// Traverse all rows.
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
	  
	  if (isPrimal(rowFeSpace, *cursor))
	    continue;


	  switch (fetiPreconditioner) {
	  case FETI_DIRICHLET:
	    colsLocal.clear();
	    colsLocalOther.clear();
	    valuesLocal.clear();
	    valuesLocalOther.clear();	    

	    // Traverse all columns.
	    for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
		 icursor != icend; ++icursor) {
	      
	      if (isPrimal(colFeSpace, col(*icursor)))
		continue;
	      
	      if (!isDual(rowFeSpace, *cursor)) {
		if (!isDual(colFeSpace, col(*icursor))) {
		  int colIndex = 
		    interiorDofMap.getLocalMatIndex(colComponent, col(*icursor));
		  colsLocal.push_back(colIndex);
		  valuesLocal.push_back(value(*icursor));
		} else {
		  int colIndex = 
		    dualDofMap.getLocalMatIndex(colComponent, col(*icursor));
		  colsLocalOther.push_back(colIndex);
		  valuesLocalOther.push_back(value(*icursor));
		}
	      } else {
		if (!isDual(colFeSpace, col(*icursor))) {
		  int colIndex = 
		    interiorDofMap.getLocalMatIndex(colComponent, col(*icursor));
		  colsLocalOther.push_back(colIndex);
		  valuesLocalOther.push_back(value(*icursor));
		} else {
		  int colIndex = 
		    dualDofMap.getLocalMatIndex(colComponent, col(*icursor));
		  colsLocal.push_back(colIndex);
		  valuesLocal.push_back(value(*icursor));
		}
	      }	 	    
	    }  // for each nnz in row

	    if (!isDual(rowFeSpace, *cursor)) {
	      int rowIndex = 
		interiorDofMap.getLocalMatIndex(rowComponent, *cursor);
	      MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(),
			   &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
	      
	      if (colsLocalOther.size()) 
		MatSetValues(mat_interior_duals, 1, &rowIndex, colsLocalOther.size(),
			     &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
	    } else {
	      int rowIndex = 
		dualDofMap.getLocalMatIndex(rowComponent, *cursor);
	      MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
			   &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);  
	      if (colsLocalOther.size()) 
		MatSetValues(mat_duals_interior, 1, &rowIndex, colsLocalOther.size(),
			     &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
	    }

	    break;


	  case FETI_LUMPED:
	    if (!isDual(rowFeSpace, *cursor))
	      continue;

	    colsLocal.clear();
	    valuesLocal.clear();

	    // Traverse all columns.
	    for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
		 icursor != icend; ++icursor) {
	      if (!isDual(colFeSpace, col(*icursor)))
		continue;

	      int colIndex = 
		dualDofMap.getLocalMatIndex(colComponent, col(*icursor));
	      colsLocal.push_back(colIndex);
	      valuesLocal.push_back(value(*icursor));
	    }

	    {
	      int rowIndex = dualDofMap.getLocalMatIndex(rowComponent, *cursor);		
	      MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
			   &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
	    }
	    break;


	  default:
	    break;
	  }	   
	} 
      }
    }
    
    // === Start global assembly procedure for preconditioner matrices. ===
    
    if (fetiPreconditioner == FETI_LUMPED) {
      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
    }
    
    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
      
      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
      
      MatAssemblyBegin(mat_interior_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_interior_duals, MAT_FINAL_ASSEMBLY);
      
      MatAssemblyBegin(mat_duals_interior, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_interior, MAT_FINAL_ASSEMBLY);
    }
    
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 03: %.5f seconds (creation of preconditioner matrices)\n",
	  MPI::Wtime() - wtime);
    }   
  }
  

  void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveFetiMatrix()");

    // The function was removed when FETI-DP was reformulated for mixed finite
    // elements. To reconstruct this function, check for revision number 2024.
    ERROR_EXIT("This function was removed!\n");
  }


  void PetscSolverFeti::solveReducedFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");

#if (DEBUG != 0)
    if (stokesMode)
      PetscSolverFetiDebug::debugNullSpace(*this, vec);
#endif

    // === Some temporary vectors. ===

    Vec tmp_b0, tmp_b1;
    localDofMap.createVec(tmp_b0, nGlobalOverallInterior);
    localDofMap.createVec(tmp_b1, nGlobalOverallInterior);

    Vec tmp_primal0, tmp_primal1;
    primalDofMap.createVec(tmp_primal0);
    primalDofMap.createVec(tmp_primal1);

    Vec tmp_lagrange, tmp_interface;
    MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange);

    // === Create RHS and solution vectors. ===

    Vec vecRhs, vecSol;
    Vec vecRhsLagrange, vecSolLagrange, vecRhsInterface, vecSolInterface;
    if (!stokesMode) {
      MatGetVecs(mat_lagrange, PETSC_NULL, &vecRhsLagrange);
      MatGetVecs(mat_lagrange, PETSC_NULL, &vecSolLagrange);

      vecRhs = vecRhsLagrange;
      vecSol = vecSolLagrange;
    } else {
      MatGetVecs(mat_lagrange, PETSC_NULL, &vecRhsLagrange);
      MatGetVecs(mat_lagrange, PETSC_NULL, &vecSolLagrange);


      interfaceDofMap.createVec(tmp_interface);
      interfaceDofMap.createVec(vecRhsInterface);
      interfaceDofMap.createVec(vecSolInterface);

      Vec vecRhsArray[2] = {vecRhsInterface, vecRhsLagrange}; 
      VecCreateNest(mpiCommGlobal, 2, PETSC_NULL, vecRhsArray, &vecRhs);

      Vec vecSolArray[2] = {vecSolInterface, vecSolLagrange}; 
      VecCreateNest(mpiCommGlobal, 2, PETSC_NULL, vecSolArray, &vecSol);

      VecAssemblyBegin(vecRhs);
      VecAssemblyEnd(vecRhs);

      VecAssemblyBegin(vecSol);
      VecAssemblyEnd(vecSol);
    }

    VecDuplicate(vecSol, &fetiKspData.draft);

    // === Create reduced RHS ===

    double wtime = MPI::Wtime();

    if (!stokesMode) {
      //    d = L inv(K_BB) f_B - L inv(K_BB) K_BPi inv(S_PiPi) [f_Pi - K_PiB inv(K_BB) f_B]
      // vecRhs = L * inv(K_BB) * f_B
      subdomain->solveGlobal(subdomain->getVecRhsInterior(), tmp_b0);
      MatMult(mat_lagrange, tmp_b0, vecRhsLagrange);
	
      // tmp_primal0 = M_PiB * inv(K_BB) * f_B
      MatMult(subdomain->getMatCoarseInterior(), tmp_b0, tmp_primal0);
      
      // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_B
      VecAXPBY(tmp_primal0, 1.0, -1.0, subdomain->getVecRhsCoarse());
      
      if (!augmentedLagrange) {
	double wtime2 = MPI::Wtime();
	// tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_B)
	KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
	if (printTimings) {
	  MSG("FETI-DP timing 09a: %.5f seconds (create rhs vector)\n", 
	      MPI::Wtime() - wtime2);
	}
	
	MatMult(subdomain->getMatInteriorCoarse(), tmp_primal0, tmp_b0);
	subdomain->solveGlobal(tmp_b0, tmp_b0);
	MatMult(mat_lagrange, tmp_b0, tmp_lagrange);
      } else {
	Vec tmp_mu;
	MatGetVecs(mat_augmented_lagrange, PETSC_NULL, &tmp_mu);
	
	MatMult(mat_augmented_lagrange, vecRhsLagrange, tmp_mu);
	VecScale(tmp_mu, -1.0);

	Vec vec_array[2] = {tmp_primal0, tmp_mu};
	Vec vec_nest;
	VecCreateNest(PETSC_COMM_WORLD, 2, PETSC_NULL, vec_array, &vec_nest);	

	KSPSolve(ksp_schur_primal, vec_nest, vec_nest);

	// 1 Step
	MatMult(subdomain->getMatInteriorCoarse(), tmp_primal0, tmp_b0);
	subdomain->solveGlobal(tmp_b0, tmp_b0);
	MatMult(mat_lagrange, tmp_b0, tmp_lagrange);

	// 2 Step
	Vec tmp_lagrange1;
	VecDuplicate(tmp_lagrange, &tmp_lagrange1);
	MatMultTranspose(mat_augmented_lagrange, tmp_mu, tmp_lagrange1);
	MatMultTranspose(mat_lagrange, tmp_lagrange1, tmp_b0);
	subdomain->solveGlobal(tmp_b0, tmp_b0);
	MatMult(mat_lagrange, tmp_b0, tmp_lagrange1);
	VecAXPY(tmp_lagrange, 1.0, tmp_lagrange1);
	VecDestroy(&tmp_lagrange1);


	VecDestroy(&tmp_mu);
	VecDestroy(&vec_nest);
      }

      VecAXPY(vecRhsLagrange, -1.0, tmp_lagrange);
    } else {
      TEST_EXIT(augmentedLagrange)("Not yet implemented!\n");

      subdomain->solveGlobal(subdomain->getVecRhsInterior(), tmp_b0);  

      MatMult(subdomain->getMatCoarseInterior(), tmp_b0, tmp_primal0);
      VecAXPBY(tmp_primal0, 1.0, -1.0, subdomain->getVecRhsCoarse());
      double wtime2 = MPI::Wtime();
      // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_B)
      KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
      if (printTimings) {
	MSG("FETI-DP timing 09a: %.5f seconds (create rhs vector)\n", 
	    MPI::Wtime() - wtime2);
      }
      
      MatMult(subdomain->getMatInteriorCoarse(), tmp_primal0, tmp_b1);
      subdomain->solveGlobal(tmp_b1, tmp_b1);  
      VecAXPY(tmp_b0, -1.0, tmp_b1);
      
      MatMult(mat_lagrange, tmp_b0, vecRhsLagrange);
      
      MatMult(subdomain->getMatCoarseInterior(1), tmp_b0, vecRhsInterface);
      MatMult(subdomain->getMatCoarse(1, 0), tmp_primal0, tmp_interface);
      VecAXPY(vecRhsInterface, 1.0, tmp_interface);
    }

    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 09: %.5f seconds (create rhs vector)\n", 
	  MPI::Wtime() - wtime);
      wtime = MPI::Wtime();
      FetiTimings::reset();
    }

    // === Optionally run some debug procedures on the FETI-DP system. ===
    PetscSolverFetiDebug::debugFeti(*this, vecRhs);

    // === Solve with FETI-DP operator. ===
    KSPSolve(ksp_feti, vecRhs, vecSol);


    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 10: %.5f seconds (application of FETI-DP operator)\n", 
	  MPI::Wtime() - wtime);
      wtime = MPI::Wtime();

      MSG("FETI-DP timing 10a: %.5f [ %.5f %.5f ] seconds (FETI-DP KSP Solve)\n", 
	  FetiTimings::fetiSolve, FetiTimings::fetiSolve01, FetiTimings::fetiSolve02);
      MSG("FETI-DP timing 10b: %.5f seconds (FETI-DP KSP Solve)\n", 
	  FetiTimings::fetiPreconditioner);
    }


    // === Solve for u_primals. ===

    if (!stokesMode) {
      MatMultTranspose(mat_lagrange, vecSol, tmp_b0);
      VecAYPX(tmp_b0, -1.0, subdomain->getVecRhsInterior());  
    } else {
      MatMultTranspose(mat_lagrange, vecSolLagrange, tmp_b0);
      VecAYPX(tmp_b0, -1.0, subdomain->getVecRhsInterior());  

      MatMult(subdomain->getMatInteriorCoarse(1), vecSolInterface, tmp_b1);
      VecAXPY(tmp_b0, -1.0, tmp_b1);
    }
  
    subdomain->solveGlobal(tmp_b0, tmp_b1);
    MatMult(subdomain->getMatCoarseInterior(), tmp_b1, tmp_primal0);
    VecAYPX(tmp_primal0, -1.0, subdomain->getVecRhsCoarse());

    if (stokesMode) {
      MatMult(subdomain->getMatCoarse(0, 1), vecSolInterface, tmp_primal1);
      VecAXPY(tmp_primal0, -1.0, tmp_primal1);
    }

    if (augmentedLagrange == false) {
      KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
    } else {
      Vec tmp_mu;
      MatGetVecs(mat_augmented_lagrange, PETSC_NULL, &tmp_mu);
      
      MatMult(mat_lagrange, tmp_b1, tmp_lagrange);
      MatMult(mat_augmented_lagrange, tmp_lagrange, tmp_mu);
      VecScale(tmp_mu, -1.0);
      
      Vec vec_array[2] = {tmp_primal0, tmp_mu};
      Vec vec_nest;
      VecCreateNest(PETSC_COMM_WORLD, 2, PETSC_NULL, vec_array, &vec_nest);	

      KSPSolve(ksp_schur_primal, vec_nest, vec_nest);

      MatMultTranspose(mat_augmented_lagrange, tmp_mu, tmp_lagrange);
      MatMultTranspose(mat_lagrange, tmp_lagrange, tmp_b1);
      VecAXPY(tmp_b0, -1.0, tmp_b1);

      VecDestroy(&tmp_mu);
      VecDestroy(&vec_nest);
    }

    // === Solve for u_b. ===

    MatMult(subdomain->getMatInteriorCoarse(), tmp_primal0, tmp_b1);
    VecAXPY(tmp_b0, -1.0, tmp_b1);
    subdomain->solveGlobal(tmp_b0, tmp_b0);


    // === And recover AMDiS solution vectors. ===

    recoverSolution(tmp_b0, tmp_primal0, vec);
    if (stokesMode)
      recoverInterfaceSolution(vecSolInterface, vec);
    
    
    // === Print timings and delete data. ===

    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 11: %.5f seconds (Inner solve and solution recovery)\n", 
	  MPI::Wtime() - wtime);
    }

    VecDestroy(&vecRhs);
    VecDestroy(&vecSol);
    VecDestroy(&tmp_b0);
    VecDestroy(&tmp_b1);
    VecDestroy(&tmp_lagrange);
    VecDestroy(&tmp_primal0);
    VecDestroy(&tmp_primal1);

    if (stokesMode) {
      VecDestroy(&tmp_interface);
      VecDestroy(&vecRhsInterface);
      VecDestroy(&vecSolInterface);      
    }
  }


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

    MatDestroy(&mat_lagrange);

    if (augmentedLagrange)
      MatDestroy(&mat_augmented_lagrange);

    // === Destroy preconditioner data structures. ===

    if (fetiPreconditioner != FETI_NONE)
      MatDestroy(&mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatDestroy(&mat_interior_interior);
      MatDestroy(&mat_interior_duals);
      MatDestroy(&mat_duals_interior);
    }

    destroySchurPrimalKsp();

    destroyFetiKsp();

    subdomain->destroyMatrixData();
  }


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

    subdomain->destroyVectorData();
  }


  void PetscSolverFeti::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverFeti::solvePetscMatrix()");

    int debug = 0;
    Parameters::get("parallel->debug feti", debug);

    if (debug)
      solveFetiMatrix(vec);
    else
      solveReducedFetiMatrix(vec);

    MeshDistributor::globalMeshDistributor->synchVector(vec);
  }


  void PetscSolverFeti::createH2vec(DOFVector<double> &vec)
  {
    FUNCNAME("PetscSolverFeti::createH2vec()");
    
    vec.set(0.0);
    
    Mesh *mesh = vec.getFeSpace()->getMesh();    
    TEST_EXIT(mesh->getDim() == 2)("This works only in 2D!\n");
    int n0 = vec.getFeSpace()->getAdmin()->getNumberOfPreDofs(VERTEX);

    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      Element *el = elInfo->getElement();
      std::vector<double> length(3, 0.0);
      double sumLength = 0.0;

      for (int i = 0; i < 3; i++) {
	int idx0 = el->getVertexOfEdge(i, 0);
	int idx1 = el->getVertexOfEdge(i, 1);
	DegreeOfFreedom dof0 = el->getDof(idx0, n0);
	DegreeOfFreedom dof1 = el->getDof(idx1, n0);

	WorldVector<double> c0 = elInfo->getCoord(idx0);
	WorldVector<double> c1 = elInfo->getCoord(idx1);
	c0 -= c1;
	length[i] = norm(c0);
	sumLength += length[i];
      }

      sumLength *= 0.5;
      double area = 
	sqrt(sumLength * (sumLength - length[0]) * (sumLength - length[1]) * (sumLength - length[2]));

      for (int i = 0; i < 3; i++) {
	DegreeOfFreedom d = el->getDof(i, n0);
	vec[d] += area;
      }
      
      elInfo = stack.traverseNext(elInfo);
    }

    meshDistributor->synchAddVector(vec);

    double scalePower = 2.0;
    Parameters::get("lp->scale power", scalePower);
    DOFIterator<double> dofIt(&vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
      DegreeOfFreedom d = dofIt.getDOFIndex();
      vec[d] = 1.0 / pow(vec[d], scalePower);      
    }
  }

}