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);
    SchurPrimalData* data = static_cast<SchurPrimalData*>(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)
    // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)

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

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

    // tmp_vec_b1 = inv(K_BB) K_BPi  inv(S_PiPi) K_PiB tmp_vec_b0
    MatMult(*(data->mat_primal_b), data->tmp_vec_b0, data->tmp_vec_primal);
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b1);

    // y = L (tmp_vec_b0 + tmp_vec_b1)
    VecAXPBY(data->tmp_vec_b0, 1.0, 1.0, data->tmp_vec_b1);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y);

    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);

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


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

    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);

    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];

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


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

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

    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

    VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);


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

    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];

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


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

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


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

    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);

    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];

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


    // === K_DD ===

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


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

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

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];

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


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
      nComponents(-1)
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
      fetiPreconditioner = FETI_LUMPED;
    } else {
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", preconditionerName.c_str());
    }
  }


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

    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("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("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("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. ===

    nLocalInterior = 0;
    for (DofMapping::iterator it = globalIndexB.begin(); 
	 it != globalIndexB.end(); ++it) {
      it->second = nLocalInterior + rStartB;
      nLocalInterior++;
    }
    nLocalDuals = duals.size();

    TEST_EXIT_DBG(nLocalInterior + 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()");

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

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

    MatCreateShell(PETSC_COMM_WORLD,
		   nRankPrimals * nComponents, nRankPrimals * nComponents,
		   nOverallPrimals * nComponents, nOverallPrimals * nComponents,
		   &schurPrimalData, 
		   &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()");

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

    VecDestroy(&schurPrimalData.tmp_vec_b);
    VecDestroy(&schurPrimalData.tmp_vec_primal);

    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
  }


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

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

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

    VecDuplicate(f_b, &(fetiData.tmp_vec_b0));
    VecDuplicate(f_b, &(fetiData.tmp_vec_b1));
    VecDuplicate(f_primal, &(fetiData.tmp_vec_primal));

    MatCreateShell(PETSC_COMM_WORLD,
		   nRankLagrange * nComponents, nRankLagrange * nComponents,
		   nOverallLagrange * nComponents, nOverallLagrange * nComponents,
		   &fetiData, &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);


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

    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }

    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_interior, "solver_interior_");
      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;
      
      VecDuplicate(f_b, &(fetiDirichletPreconData.tmp_vec_b));      
      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));
      
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
      
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;

      VecDuplicate(f_b, &(fetiLumpedPreconData.tmp_vec_b));
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals1));

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
      break;
    }
  }
  

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

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

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

    VecDestroy(&fetiData.tmp_vec_b0);
    VecDestroy(&fetiData.tmp_vec_b1);
    VecDestroy(&fetiData.tmp_vec_primal);
    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.mat_lagrange_scaled = NULL;
      fetiLumpedPreconData.mat_duals_duals = NULL;

      VecDestroy(&fetiLumpedPreconData.tmp_vec_b);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1);
      break;
    }
  }


  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;
    int nRowsInterior = nLocalInterior * nComponents;
    int nRowsDual = nLocalDuals * 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);


    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsDual, nRowsDual, 100, PETSC_NULL,
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsInterior, nRowsInterior, 100, PETSC_NULL,
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsInterior, nRowsDual, 100, PETSC_NULL,
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
		      nRowsDual, nRowsInterior, 100, PETSC_NULL,
		      &mat_duals_interior);
    }

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

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

    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();
	  colsOther.clear();
	  values.clear();	  
	  valuesOther.clear();

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

	    bool colPrimal = primals.count(col(*icursor)) != 0;

	    if (colPrimal) {
	      // 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));
	      }
	    }



	    // === For preconditioner ===

	    if (!rowPrimal && !colPrimal) {
	      int rowIndex = globalIndexB[*cursor] - rStartB;
	      int colIndex = globalIndexB[col(*icursor)] - rStartB;
		
	      if (rowIndex < nLocalInterior) {
		if (colIndex < nLocalInterior) {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB) * nComponents + j;

		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		} else {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		}
	      } else {
		if (colIndex < nLocalInterior) {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB) * nComponents + j;

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		} else {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;

		  colsLocal.push_back(colIndex2);
		  valuesLocal.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);
	  }

	  // === For preconditioner ===

	  if (!rowPrimal) {
	    switch (fetiPreconditioner) {
	    case FETI_DIRICHLET:
	      {
		int rowIndex = globalIndexB[*cursor] - rStartB;
		
		if (rowIndex < nLocalInterior) {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB) * nComponents + i;
		  
		  MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
		  
		  if (colsLocalOther.size()) 
		    MatSetValues(mat_interior_duals, 1, &rowIndex2, colsLocalOther.size(),
				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
		} else {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
		  
		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
		  
		  if (colsLocalOther.size()) 
		    MatSetValues(mat_duals_interior, 1, &rowIndex2, colsLocalOther.size(),
				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
		  
		}
	      }
	      break;

	    case FETI_LUMPED:
	      {
		int rowIndex = globalIndexB[*cursor] - rStartB;
		
		if (rowIndex >= nLocalInterior) {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
		  
		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);		
		}
	      }		
	      break;
	    }	  
	  }
	} 
      }
    }
    

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


    // === Start global assembly procedure for preconditioner matrices. ===

    if (fetiPreconditioner != FETI_NONE) {
      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_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);
    }


    // === 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 + nRankLagrange) * nComponents;
    int nOverallNest = (nOverallB + nOverallPrimals + nOverallLagrange) * nComponents;
    int rStartNest = (rStartB + rStartPrimals + rStartLagrange) * nComponents;

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