PetscSolverFeti.cc

//
// 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 "parallel/PetscSolverFeti.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

  using namespace std;


#ifdef HAVE_PETSC_DEV 
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
    PetscSchurPrimalData* data = static_cast<PetscSchurPrimalData*>(ctx);

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);

    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(*(data->mat_primal_primal), x, y);
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
    // F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)

    void *ctx;
    MatShellGetContext(mat, &ctx);
    PetscFetiData* data = static_cast<PetscFetiData*>(ctx);

    // y = L inv(K_BB) trans(L) x
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);

    // tmp_vec_primal = inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);

    // tmp_vec_lagrange = L inv(K_BB) K_BPi tmp_vec_primal
    //                  = L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b);
    KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);

    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);

    return 0;
  }


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

    TEST_EXIT(meshDistributor->getMesh()->getDim() == 2)
      ("Works for 2D problems only!");

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
   
    createPrimals();

    createDuals();

    createLagrange();

    createIndexB();
  }


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

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

    primals.clear();
    DofContainerSet& vertices = 
      meshDistributor->getBoundaryDofInfo().geoDofs[VERTEX];
    TEST_EXIT_DBG(vertices.size())("No primal vertices on this rank!\n");
    for (DofContainerSet::iterator it = vertices.begin(); 
	 it != vertices.end(); ++it)
      primals.insert(**it);


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

    globalPrimalIndex.clear();
    nRankPrimals = 0;
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it)) {
	globalPrimalIndex[*it] = nRankPrimals;
	nRankPrimals++;
      }


    // === Get overall number of primals and rank's displacement in the ===
    // === numbering of the primals.                                    ===

    nOverallPrimals = 0;
    rStartPrimals = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankPrimals, rStartPrimals, nOverallPrimals);


    // === Create global primal index for all primals. ===

    for (DofMapping::iterator it = globalPrimalIndex.begin();
	 it != globalPrimalIndex.end(); ++it)
      it->second += rStartPrimals;

    MSG_DBG("nRankPrimals = %d   nOverallPrimals = %d\n",
	    nRankPrimals, nOverallPrimals);


    // === Communicate primal's global index from ranks that own the     ===
    // === primals to ranks that contain this primals but are not owning ===
    // === them.                                                         ===

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (globalPrimalIndex.count(**dofIt))
	  stdMpi.getSendData(it->first).push_back(globalPrimalIndex[**dofIt]);
    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvFromRank = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false) {
	  recvFromRank = true;
	  break;
	}

      if (recvFromRank) 
	stdMpi.recv(it->first);
    }
    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int i = 0;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false)
	  globalPrimalIndex[**dofIt] = stdMpi.getRecvData(it->first)[i++];
      }
    }

    TEST_EXIT_DBG(primals.size() == globalPrimalIndex.size())
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
       primals.size(), globalPrimalIndex.size());


    TEST_EXIT_DBG(nOverallPrimals > 0)
      ("There are zero primal nodes in domain!\n");
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
    
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===

    boundaryDofRanks.clear();

    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it) {
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	// If DOF is not primal, i.e., its a dual node
	if (primals.count(**dofIt) == 0) {
	  boundaryDofRanks[**dofIt].insert(mpiRank);
	  boundaryDofRanks[**dofIt].insert(it->first);
	}
      }
    }


    // === 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 (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) == 0)
	  stdMpi.getSendData(it->first).push_back(boundaryDofRanks[**dofIt]);

    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvFromRank = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) == 0) {
	  recvFromRank = true;
	  break;
	}

      if (recvFromRank)
	stdMpi.recv(it->first);
    }
    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int i = 0;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)	
	if (primals.count(**dofIt) == 0)
	  boundaryDofRanks[**dofIt] = stdMpi.getRecvData(it->first)[i++];	      
    }


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

    duals.clear();
    globalDualIndex.clear();

    int nRankAllDofs = meshDistributor->getFeSpace()->getAdmin()->getUsedDofs();
    nRankB = nRankAllDofs - primals.size();
    nOverallB = 0;
    rStartB = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankB, rStartB, nOverallB);
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(allBoundaryDofs);
    int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size();

    int nRankDuals = 0;
    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it) {
      if (primals.count(**it) == 0) {
	duals.insert(**it);
	globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals;
	nRankDuals++;
      }
    }

    int nOverallDuals = nRankDuals;
    mpi::globalAdd(nOverallDuals);

    MSG_DBG("nRankDuals = %d   nOverallDuals = %d\n",
	    nRankDuals, nOverallDuals);
  }

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

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

    nRankLagrange = 0;
    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
      if (meshDistributor->getIsRankDof(*it)) {
	dofFirstLagrange[*it] = nRankLagrange;
	int degree = boundaryDofRanks[*it].size();
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }


    // === Get the overall number of Lagrange constraints and create the ===
    // === mapping dofFirstLagrange, that defines for each dual boundary ===
    // === node the first Lagrange constraint global index.              ===

    nOverallLagrange = 0;
    rStartLagrange = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankLagrange, rStartLagrange, nOverallLagrange);

    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it))
	dofFirstLagrange[*it] += rStartLagrange;

    MSG_DBG("nRankLagrange = %d  nOverallLagrange = %d\n",
	    nRankLagrange, nOverallLagrange);


    // === Communicate dofFirstLagrange to all other ranks. ===

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (primals.count(**dofIt) == 0) {
	  TEST_EXIT_DBG(dofFirstLagrange.count(**dofIt))("Should not happen!\n");
	  stdMpi.getSendData(it->first).push_back(dofFirstLagrange[**dofIt]);
	}
      }
    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvData = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) == 0) {
	  recvData = true;
	  break;
	}
	  
      if (recvData)
	stdMpi.recv(it->first);
    }

    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int counter = 0;
      for (unsigned int i = 0; i < it->second.size(); i++)
	if (primals.count(*(it->second[i])) == 0)
	  dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++];
    }     
  }


  void PetscSolverFeti::createIndexB()
  {
    FUNCNAME("PetscSolverFeti::createIndeB()");

    globalIndexB.clear();
    DOFAdmin* admin = meshDistributor->getFeSpace()->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.                                 ===
    
    for (int i = 0; i < admin->getUsedSize(); i++)
      if (admin->isDofFree(i) == false && primals.count(i) == 0)
	if (duals.count(i) == 0 && primals.count(i) == 0)
	  globalIndexB[i] = -1;


    // === Get correct index for all interior nodes. ===

    int nInterior = 0;
    for (DofMapping::iterator it = globalIndexB.begin(); 
	 it != globalIndexB.end(); ++it) {
      it->second = nInterior + rStartB;
      nInterior++;
    }

    TEST_EXIT_DBG(nInterior + primals.size() + duals.size() == 
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");


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

    for (DofIndexSet::iterator it = duals.begin();
	 it != duals.end(); ++it)
      globalIndexB[*it] = globalDualIndex[*it];
  }


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

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

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_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 (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
      TEST_EXIT_DBG(dofFirstLagrange.count(*it))("Should not happen!\n");
      TEST_EXIT_DBG(boundaryDofRanks.count(*it))("Should not happen!\n");

      // Global index of the first Lagrange constriant for this node.
      int index = dofFirstLagrange[*it];
      // Copy set of all ranks that contain this dual node.
      vector<int> W(boundaryDofRanks[*it].begin(), boundaryDofRanks[*it].end());
      // Number of ranks that contain this dual node.
      int degree = W.size();

      TEST_EXIT_DBG(globalDualIndex.count(*it))("Should not happen!\n");
      int dualCol = globalDualIndex[*it];

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
	      int colIndex = dualCol * nComponents + k;
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
	    }
	  }

	  index++;
	}
      }
    }

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


  void PetscSolverFeti::createSchurPrimalKsp()
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

    petscSchurPrimalData.mat_primal_primal = &mat_primal_primal;
    petscSchurPrimalData.mat_primal_b = &mat_primal_b;
    petscSchurPrimalData.mat_b_primal = &mat_b_primal;
    petscSchurPrimalData.ksp_b = &ksp_b;

    VecDuplicate(f_b, &(petscSchurPrimalData.tmp_vec_b));
    VecDuplicate(f_primal, &(petscSchurPrimalData.tmp_vec_primal));

    MatCreateShell(PETSC_COMM_WORLD,
		   nRankPrimals * nComponents, nRankPrimals * nComponents,
		   nOverallPrimals * nComponents, nOverallPrimals * nComponents,
		   &petscSchurPrimalData, 
		   &mat_schur_primal);
    MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			 (void(*)(void))petscMultMatSchurPrimal);

    KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
    KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
    KSPSetFromOptions(ksp_schur_primal);
  }


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

    petscSchurPrimalData.mat_primal_primal = PETSC_NULL;
    petscSchurPrimalData.mat_primal_b = PETSC_NULL;
    petscSchurPrimalData.mat_b_primal = PETSC_NULL;
    petscSchurPrimalData.ksp_b = PETSC_NULL;

    VecDestroy(petscSchurPrimalData.tmp_vec_b);
    VecDestroy(petscSchurPrimalData.tmp_vec_primal);

    MatDestroy(mat_schur_primal);
    KSPDestroy(ksp_schur_primal);
  }


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

    petscFetiData.mat_primal_primal = &mat_primal_primal;
    petscFetiData.mat_primal_b = &mat_primal_b;
    petscFetiData.mat_b_primal = &mat_b_primal;
    petscFetiData.mat_lagrange = &mat_lagrange;
    petscFetiData.ksp_b = &ksp_b;
    petscFetiData.ksp_schur_primal = &ksp_schur_primal;


    VecDuplicate(f_b, &(petscFetiData.tmp_vec_b));
    VecDuplicate(f_primal, &(petscFetiData.tmp_vec_primal));
    MatGetVecs(mat_lagrange, PETSC_NULL, &(petscFetiData.tmp_vec_lagrange));


    MatCreateShell(PETSC_COMM_WORLD,
		   nRankLagrange, nRankLagrange,
		   nOverallLagrange, nOverallLagrange,
		   &petscFetiData, &mat_feti);
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


    KSPCreate(PETSC_COMM_WORLD, &ksp_feti);
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_feti, "solver_feti_");
    KSPSetFromOptions(ksp_feti);
  }
  

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

    petscFetiData.mat_primal_primal = PETSC_NULL;
    petscFetiData.mat_primal_b = PETSC_NULL;
    petscFetiData.mat_b_primal = PETSC_NULL;
    petscFetiData.mat_lagrange = PETSC_NULL;
    petscFetiData.ksp_b = PETSC_NULL;
    petscFetiData.ksp_schur_primal = PETSC_NULL;

    VecDestroy(petscFetiData.tmp_vec_b);
    VecDestroy(petscFetiData.tmp_vec_primal);
    VecDestroy(petscFetiData.tmp_vec_lagrange);

    MatDestroy(mat_feti);
    KSPDestroy(ksp_feti);
  }


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


    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
    
    vector<PetscInt> globalIsIndex, localIsIndex;
    globalIsIndex.reserve(globalPrimalIndex.size() * nComponents);
    localIsIndex.reserve(globalPrimalIndex.size() * nComponents);

    {
      int counter = 0;
      for (DofMapping::iterator it = globalPrimalIndex.begin();
	   it != globalPrimalIndex.end(); ++it) {
	for (int i = 0; i < nComponents; i++) {
	  globalIsIndex.push_back(it->second * nComponents + i);
	  localIsIndex.push_back(counter++);
	}
      }
    }
    
    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. ===

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));

      for (DofMapping::iterator it = globalIndexB.begin();
	   it != globalIndexB.end(); ++it) {
	int petscIndex = (it->second - rStartB) * nComponents + i;
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
      for (DofMapping::iterator it = globalPrimalIndex.begin();
	   it != globalPrimalIndex.end(); ++it) {
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }



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


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

    nComponents = vec->getSize();

    // === Create all sets and indices. ===

    updateDofData();


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

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
		    nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_b_b);

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
		    10, PETSC_NULL, 10, PETSC_NULL, &mat_primal_primal);

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_b_primal);

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
		    100, PETSC_NULL, 100, PETSC_NULL, &mat_primal_b);

    
    // === 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> cols, colsOther;
    vector<double> values, valuesOther;
    cols.reserve(300);
    colsOther.reserve(300);
    values.reserve(300);
    valuesOther.reserve(300);


    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	if (!(*mat)[i][j])
	  continue;

	traits::col<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
	traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
	
	// Traverse all rows.
	for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()), 
	       cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
	  bool rowPrimal = primals.count(*cursor) != 0;
	  
	  cols.clear();
	  values.clear();
	  
	  colsOther.clear();
	  valuesOther.clear();
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

	    if (primals.count(col(*icursor)) != 0) {
	      // Column is a primal variable.

	      TEST_EXIT_DBG(globalPrimalIndex.count(col(*icursor)))
		("No global primal index for DOF %d!\n", col(*icursor));
	      
	      int colIndex = globalPrimalIndex[col(*icursor)] * nComponents + j;
	      
	      if (rowPrimal) {
		cols.push_back(colIndex);
		values.push_back(value(*icursor));
	      } else {
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      // Column is not a primal variable.

	      TEST_EXIT_DBG(globalIndexB.count(col(*icursor)))
		("No global B index for DOF %d!\n", col(*icursor));
	      
	      int colIndex = globalIndexB[col(*icursor)] * nComponents + j;

	      if (rowPrimal) {
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(colIndex);
		values.push_back(value(*icursor));
	      }
	    }
	  }

	  if (rowPrimal) {
	    TEST_EXIT_DBG(globalPrimalIndex.count(*cursor))
	      ("Should not happen!\n");

	    int rowIndex = globalPrimalIndex[*cursor] * nComponents + i;
	    MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size())
	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  } else {
	    TEST_EXIT_DBG(globalIndexB.count(*cursor))
	      ("Should not happen!\n");

	    int rowIndex = globalIndexB[*cursor] * nComponents + i;
	    MatSetValues(mat_b_b, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size())
	      MatSetValues(mat_b_primal, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  }
	} 
      }
    }
    

    // === Start global assembly procedure. ===

    MatAssemblyBegin(mat_b_b, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_b_b, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_primal_primal, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_primal_primal, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_b_primal, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_b_primal, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_primal_b, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_primal_b, MAT_FINAL_ASSEMBLY);
	  

    // === Create and fill PETSc's right hand side vectors. ===

    VecCreate(PETSC_COMM_WORLD, &f_b);
    VecSetSizes(f_b, nRankB * nComponents, nOverallB * nComponents);
    VecSetType(f_b, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &f_primal);
    VecSetSizes(f_primal, nRankPrimals * nComponents, 
		nOverallPrimals * nComponents);
    VecSetType(f_primal, VECMPI);
    
    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
      for (dofIt.reset(); !dofIt.end(); ++dofIt) {
	int index = dofIt.getDOFIndex();
	if (primals.count(index)) {
	  TEST_EXIT_DBG(globalPrimalIndex.count(index))
	    ("Should not happen!\n");

	  index = globalPrimalIndex[index] * nComponents + i;
	  double value = *dofIt;
	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
	} else {
	  TEST_EXIT_DBG(globalIndexB.count(index))
	    ("Should not happen!\n");

	  index = globalIndexB[index] * nComponents + i;
	  double value = *dofIt;
	  VecSetValues(f_b, 1, &index, &value, ADD_VALUES);
	}      
      }
    }

    VecAssemblyBegin(f_b);
    VecAssemblyEnd(f_b);

    VecAssemblyBegin(f_primal);
    VecAssemblyEnd(f_primal);


    // === Create and fill PETSc matrix for Lagrange constraints. ===

    createMatLagrange();

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


    // === Create PETSc solver for the FETI-DP operator. ===

    createFetiKsp();
  }


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

    // Create transpose of Lagrange matrix.
    Mat mat_lagrange_transpose;
    MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);


    // === Create nested matrix which will contain the overall FETI system. ===

    Mat A;
    Mat nestedA[3][3];
    nestedA[0][0] = mat_b_b;
    nestedA[0][1] = mat_b_primal;
    nestedA[0][2] = mat_lagrange_transpose;
    nestedA[1][0] = mat_primal_b;
    nestedA[1][1] = mat_primal_primal;
    nestedA[1][2] = PETSC_NULL;
    nestedA[2][0] = mat_lagrange;
    nestedA[2][1] = PETSC_NULL;
    nestedA[2][2] = PETSC_NULL;

    MatCreateNest(PETSC_COMM_WORLD, 3, PETSC_NULL, 3, PETSC_NULL, &(nestedA[0][0]), &A);

    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  


    int nRankNest = (nRankB + nRankPrimals) * nComponents + nRankLagrange;    
    int nOverallNest = (nOverallB + nOverallPrimals) * nComponents + nOverallLagrange;
    int rStartNest = (rStartB + rStartPrimals) * nComponents + rStartLagrange;

    {
      // === Test some matrix sizes. ===

      int matRow, matCol;
      MatGetLocalSize(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == nRankNest)("Should not happen!\n");
      mpi::globalAdd(matRow);
      TEST_EXIT_DBG(matRow == nOverallNest)("Should not happen!\n");

      MatGetOwnershipRange(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == rStartNest)("Should not happen!\n");
    }

    // === Create rhs and solution vectors for the overall FETI system. ===

    Vec f;
    VecCreate(PETSC_COMM_WORLD, &f);
    VecSetSizes(f, nRankNest, nOverallNest);
    VecSetType(f, VECMPI);

    Vec b;
    VecDuplicate(f, &b);

    
    // === Fill rhs vector by coping the primal and non primal PETSc vectors. ===

    PetscScalar *local_f_b;
    VecGetArray(f_b, &local_f_b);

    PetscScalar *local_f_primal;
    VecGetArray(f_primal, &local_f_primal);

    {
      int size;
      VecGetLocalSize(f_b, &size);
      TEST_EXIT_DBG(size == nRankB * nComponents)("Should not happen!\n");
      VecGetLocalSize(f_primal, &size);
      TEST_EXIT_DBG(size == nRankPrimals * nComponents)("Should not happen!\n");
    }

    PetscScalar *local_f;
    VecGetArray(f, &local_f);

    int index = 0;
    for (int i = 0; i < nRankB * nComponents; i++)
      local_f[index++] = local_f_b[i];
    for (int i = 0; i < nRankPrimals * nComponents; i++)
      local_f[index++] = local_f_primal[i];

    VecRestoreArray(f, &local_f);  
    VecRestoreArray(f_b, &local_f_b);
    VecRestoreArray(f_primal, &local_f_primal);

    
    // === Create solver and solve the overall FETI system. ===

    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, A, A, SAME_NONZERO_PATTERN);
    KSPSetFromOptions(ksp);


    KSPSolve(ksp, f, b);


    // === Reconstruct FETI solution vectors. ===
    
    Vec u_b, u_primal;
    VecDuplicate(f_b, &u_b);
    VecDuplicate(f_primal, &u_primal);
    

    PetscScalar *local_b;
    VecGetArray(b, &local_b);

    PetscScalar *local_u_b;
    VecGetArray(u_b, &local_u_b);

    PetscScalar *local_u_primal;
    VecGetArray(u_primal, &local_u_primal);

    index = 0;
    for (int i = 0; i < nRankB * nComponents; i++)
      local_u_b[i] = local_b[index++];
    for (int i = 0; i < nRankPrimals * nComponents; i++)
      local_u_primal[i] = local_b[index++];

    VecRestoreArray(f, &local_b);
    VecRestoreArray(u_b, &local_u_b);
    VecRestoreArray(u_primal, &local_u_primal);

    recoverSolution(u_b, u_primal, vec);

    VecDestroy(u_b);
    VecDestroy(u_primal);
    VecDestroy(b);
    VecDestroy(f);

    KSPDestroy(ksp);
  }


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

    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(PETSC_COMM_WORLD, &ksp_b);
    KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_b, "solver_b_");
    KSPSetFromOptions(ksp_b);

    // RHS and solution vector.
    Vec vec_rhs;

    // Some temporary vectors.
    Vec tmp_b0, tmp_b1, tmp_lagrange0, tmp_primal0, tmp_primal1;
    MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange0);
    MatGetVecs(mat_lagrange, PETSC_NULL, &vec_rhs);
    MatGetVecs(mat_b_b, PETSC_NULL, &tmp_b0);
    MatGetVecs(mat_b_b, PETSC_NULL, &tmp_b1);
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal0);
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal1);


    // === Create new rhs ===

    // vec_rhs = L * inv(K_BB) * f_b
    KSPSolve(ksp_b, f_b, tmp_b0);
    MatMult(mat_lagrange, tmp_b0, vec_rhs);

    // tmp_primal0 = M_PiB * inv(K_BB) * f_b
    MatMult(mat_primal_b, tmp_b0, tmp_primal0);

    // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_b
    VecAXPBY(tmp_primal0, -1.0, 1.0, f_primal);

    // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_b)
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);

    //
    MatMult(mat_b_primal, tmp_primal0, tmp_b0);
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
    MatMult(mat_lagrange, tmp_b0, tmp_lagrange0);

    //
    VecAXPBY(vec_rhs, 1.0, 1.0, tmp_lagrange0);


    // === Solve with FETI-DP operator. ===

    KSPSolve(ksp_feti, vec_rhs, vec_rhs);

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

    VecCopy(f_primal, tmp_primal0);

    KSPSolve(ksp_b, f_b, tmp_b0);
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);

    VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1);

    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b0);
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);

    VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1);
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);

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

    VecCopy(f_b, tmp_b0);
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);

    MatMult(mat_b_primal, tmp_primal0, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);

    KSPSolve(ksp_b, tmp_b0, tmp_b0);


    // === And recover AMDiS solution vectors. ===
    
    recoverSolution(tmp_b0, tmp_primal0, vec);


    // === Destroy all data structures. ===
    
    VecDestroy(vec_rhs);
    VecDestroy(tmp_b0);
    VecDestroy(tmp_b1);
    VecDestroy(tmp_lagrange0);
    VecDestroy(tmp_primal0);
    VecDestroy(tmp_primal1);
	    

    KSPDestroy(ksp_b);

    MatDestroy(mat_b_b);
    MatDestroy(mat_primal_primal);
    MatDestroy(mat_b_primal);
    MatDestroy(mat_primal_b);
    MatDestroy(mat_lagrange);

    VecDestroy(f_b);
    VecDestroy(f_primal);

    destroySchurPrimalKsp();

    destroyFetiKsp();   
  }


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

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

    if (debug) {
      WARNING("FETI matrix is solved globally, thus without reducing to the lagrange multipliers!\n");

      solveFetiMatrix(vec);
    } else {
      solveReducedFetiMatrix(vec);
    }
      
  }
#endif

}