feti.h

Go to the documentation of this file.

/*
 * Copyright (c) 2010-2015:  G-CSC, Goethe University Frankfurt
 * Authors: Ingo Heppner, Andreas Vogel
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
/* TODO:
 * If finished: Check for unnecessary "set zero" operations!
 *
 */
 
#ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__FETI__
#define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__FETI__
 
 
#ifdef UG_PARALLEL
 
#include <iostream>
#include <sstream>
#include <string>
#include "lib_algebra/operator/interface/linear_operator.h"
#include "lib_algebra/operator/interface/linear_operator_inverse.h"
#include "lib_algebra/operator/interface/matrix_operator.h"
#include "lib_algebra/operator/interface/matrix_operator_inverse.h"
#include "lib_algebra/parallelization/parallelization.h"
#include "lib_algebra/operator/debug_writer.h"
#include "pcl/pcl.h"
 
/* 
 * Outline of FETI-DP algorithm (after Klawonn, A., Widlund, O., Dryja, M.:
 * "Dual-Primal FETI methods for three-dimensional elliptic problems with
 * heterogeneous coefficients", SIAM J. Numer. Anal., 40/2002a, 159--179.):
 *
 * 1. eliminate the primal variables associated with the interior nodes,
 *    the vertex nodes designated as primal (as well as the Lagrange
 *    multipliers related to the optional constraints - if any)
 *    ==> reduced system with Schur complement \tilde{S} related to the subspace
 *        \tilde{W}_{\Delta}, and a reduced right hand side \tilde{f}_{\Delta}.
 *        
 * 2. Reformulate as a variational problem, with a vector of Lagrange multipliers
 *    ==> saddle point problem.
 *
 * 3. Eliminate subvectors u_{\Delta}
 *    ==> system for the dual variables: F \lambda = d
 *
 * 4. Instead of solving ... action of the inverse of the Schur complement can be
 *    obtained by solving the entire linear system from which it originates after
 *    augmenting the given right hand side with zeros.
 *    Group all the interior and dual variables of each subdomain together
 *    and factor the resulting blocks in parallel across the subdomains using a
 *     good ordering algorithm.
 *
 *    The contributions of the remaining Schur complement, of the primal variables,
 *    can also be computed locally prior to subassembly and factorization of this
 *    final, global part of the linear system of equations.
 */
 
namespace ug{
 
 
template <typename TAlgebra>
class FetiLayouts
{
  public:
  //  Algebra type
    typedef TAlgebra algebra_type;
 
  //  Vector type
    typedef typename TAlgebra::vector_type vector_type;
 
  //  Matrix type
    typedef typename TAlgebra::matrix_type matrix_type;
 
  public:
    FetiLayouts() {}
 
  //  assigns standard layouts, creates the feti layouts
    void create_layouts(ConstSmartPtr<AlgebraLayouts> stdLayouts,
                        size_t numIndices,
                        pcl::IDomainDecompositionInfo& DDInfo,
                        bool bDebug = false)
    {
    //  if no indices, nothing to do
      if(numIndices == 0) return;
 
      m_spStdLayouts = stdLayouts;
      m_spInnerLayouts = SmartPtr<AlgebraLayouts>(new AlgebraLayouts());
 
    //  create FETI Layouts:
 //   \todo: For some documentation info see mail by S. Reiter, 30. Januar 2011 16:10:52 MEZ
//             (if the information therein is not already outdated ...)
      BuildDomainDecompositionLayouts(
              m_masterDualLayout, m_slaveDualLayout,
              m_spInnerLayouts->master(), m_spInnerLayouts->slave(),
              m_masterDualNbrLayout, m_slaveDualNbrLayout,
              m_masterPrimalLayout, m_slavePrimalLayout,
              stdLayouts->master(), stdLayouts->slave(),
              (int)(numIndices - 1), DDInfo);
      UG_LOG("[BuildDomainDecompositionLayouts done]");
 
    //  create intra feti subdomain communicator
      int localSubdomID = DDInfo.map_proc_id_to_subdomain_id(pcl::ProcRank());
      pcl::ProcessCommunicator worldComm;
      for(int i = 0; i < DDInfo.get_num_subdomains(); ++i){
        if(localSubdomID == i)
          m_spInnerLayouts->proc_comm() = worldComm.create_sub_communicator(true);
        else
          worldComm.create_sub_communicator(false);
      }
      UG_LOG("[intra feti sd comms created]");
 
    //  create set of unique indices (to avoid doubles due to several interfaces)
    //  collect from dual layouts
      CollectElements(m_vUniqueMasterDualIndex, m_masterDualLayout);
      sort_and_remove_doubles(m_vUniqueMasterDualIndex);
 
      CollectElements(m_vUniqueSlaveDualIndex, m_slaveDualLayout);
      sort_and_remove_doubles(m_vUniqueSlaveDualIndex);
 
      CollectElements(m_vUniqueMasterDualNbrIndex, m_masterDualNbrLayout);
      sort_and_remove_doubles(m_vUniqueMasterDualNbrIndex);
 
      CollectElements(m_vUniqueSlaveDualNbrIndex, m_slaveDualNbrLayout);
      sort_and_remove_doubles(m_vUniqueSlaveDualNbrIndex);
      UG_LOG("[CollectElements for duals done]");
 
    //  create union of all dual unknowns
      add_merge_sort_remove_doubles(m_vUniqueDualIndex, m_vUniqueMasterDualIndex);
      add_merge_sort_remove_doubles(m_vUniqueDualIndex, m_vUniqueSlaveDualIndex);
      add_merge_sort_remove_doubles(m_vUniqueDualIndex, m_vUniqueMasterDualNbrIndex);
      add_merge_sort_remove_doubles(m_vUniqueDualIndex, m_vUniqueSlaveDualNbrIndex);
      UG_LOG("[union of duals created]");
 
    //  collect from primal layouts
      CollectElements(m_vUniquePrimalMasterIndex, m_masterPrimalLayout);
      sort_and_remove_doubles(m_vUniquePrimalMasterIndex);
 
      CollectElements(m_vUniquePrimalSlaveIndex, m_slavePrimalLayout);
      sort_and_remove_doubles(m_vUniquePrimalSlaveIndex);
 
    //  create union of all primal unknowns
      add_merge_sort_remove_doubles(m_vUniquePrimalIndex, m_vUniquePrimalMasterIndex);
      add_merge_sort_remove_doubles(m_vUniquePrimalIndex, m_vUniquePrimalSlaveIndex);
      UG_LOG("[CollectElements for primals done] ");
 
    //  test output
      if(bDebug && false)
      {
        pcl::InterfaceCommunicator<IndexLayout> comTmp;
        UG_LOG("STANDARD LAYOUTS:\n");
        PrintLayout(m_spStdLayouts->proc_comm(), comTmp, m_spStdLayouts->master(), m_spStdLayouts->slave());
        UG_LOG("INNER LAYOUTS:\n");
        PrintLayout(m_spStdLayouts->proc_comm(), comTmp, m_spInnerLayouts->master(), m_spInnerLayouts->slave());
        UG_LOG("PRIMAL LAYOUTS:\n");
        PrintLayout(m_spStdLayouts->proc_comm(), comTmp, m_masterPrimalLayout, m_slavePrimalLayout);
        UG_LOG("DUAL LAYOUTS:\n");
        PrintLayout(m_spStdLayouts->proc_comm(), comTmp, m_masterDualLayout, m_slaveDualLayout);
        UG_LOG("DUAL NBR LAYOUTS:\n");
        PrintLayout(m_spStdLayouts->proc_comm(), comTmp, m_masterDualNbrLayout, m_slaveDualNbrLayout);
      }
    }
 
  //  tests the previously created layouts:
    void test_layouts(bool bPrint)
    {
      if (m_spInnerLayouts.invalid() || m_spInnerLayouts.invalid()) {
        UG_LOG("LAYOUTS not yet created!\n");
        return;
      }
      pcl::InterfaceCommunicator<IndexLayout> comTmp;
 
      UG_LOG("TEST STANDARD LAYOUTS:\n");
      if (TestLayout(m_spStdLayouts->proc_comm(), comTmp, m_spStdLayouts->master(), m_spStdLayouts->slave(), bPrint) != true) {
        UG_LOG("STANDARD LAYOUTS inconsistent!\n");
      } else {
        UG_LOG("STANDARD LAYOUTS are consistent!\n");
      }
 
      UG_LOG("TEST INNER LAYOUTS:\n");
      if (TestLayout(m_spStdLayouts->proc_comm(), comTmp, m_spInnerLayouts->master(), m_spInnerLayouts->slave(), bPrint) != true) {
        UG_LOG("INNER LAYOUTS inconsistent!\n");
      } else {
        UG_LOG("INNER LAYOUTS are consistent!\n");
      }
 
      UG_LOG("TEST PRIMAL LAYOUTS:\n");
      if (TestLayout(m_spStdLayouts->proc_comm(), comTmp, m_masterPrimalLayout, m_slavePrimalLayout, bPrint) != true) {
        UG_LOG("PRIMAL LAYOUTS inconsistent!\n");
      } else {
        UG_LOG("PRIMAL LAYOUTS are consistent!\n");
      }
 
      UG_LOG("TEST DUAL LAYOUTS:\n");
      if (TestLayout(m_spStdLayouts->proc_comm(), comTmp, m_masterDualLayout, m_slaveDualLayout, bPrint) != true) {
        UG_LOG("DUAL LAYOUTS inconsistent!\n");
      } else {
        UG_LOG("DUAL LAYOUTS are consistent!\n");
      }
 
      UG_LOG("TEST DUAL NBR LAYOUTS:\n");
      if (TestLayout(m_spStdLayouts->proc_comm(), comTmp, m_masterDualNbrLayout, m_slaveDualNbrLayout, bPrint) != true) {
        UG_LOG("DUAL NBR LAYOUTS inconsistent!\n");
      } else {
        UG_LOG("DUAL NBR LAYOUTS are consistent!\n");
      }
    }
 
  //  standard layouts
    const IndexLayout& get_std_master_layout() {return m_spStdLayouts->master();}
    const IndexLayout& get_std_slave_layout() {return m_spStdLayouts->slave();}
    const pcl::ProcessCommunicator& get_std_process_communicator() {return m_spStdLayouts->proc_comm();}
 
  //  intra subdomain layouts
    IndexLayout& get_intra_sd_master_layout() {return m_spInnerLayouts->master();}
    IndexLayout& get_intra_sd_slave_layout() {return m_spInnerLayouts->slave();}
    pcl::ProcessCommunicator& get_intra_sd_process_communicator() {return m_spInnerLayouts->proc_comm();}
 
  //  dual layouts
    IndexLayout& get_dual_master_layout() {return m_masterDualLayout;}
    IndexLayout& get_dual_slave_layout() {return m_slaveDualLayout;}
    const std::vector<IndexLayout::Element>& get_dual_master_indices() const {return m_vUniqueMasterDualIndex;}
    const std::vector<IndexLayout::Element>& get_dual_slave_indices() const {return m_vUniqueSlaveDualIndex;}
 
  //  dual nbr layouts
    IndexLayout& get_dual_nbr_master_layout() {return m_masterDualNbrLayout;}
    IndexLayout& get_dual_nbr_slave_layout() {return m_slaveDualNbrLayout;}
    const std::vector<IndexLayout::Element>& get_dual_nbr_master_indices() const {return m_vUniqueMasterDualNbrIndex;}
    const std::vector<IndexLayout::Element>& get_dual_nbr_slave_indices() const {return m_vUniqueSlaveDualNbrIndex;}
 
  //  primal layouts
    IndexLayout& get_primal_master_layout() {return m_masterPrimalLayout;}
    IndexLayout& get_primal_slave_layout() {return m_slavePrimalLayout;}
    const std::vector<IndexLayout::Element>& get_primal_master_indices() const {return m_vUniquePrimalMasterIndex;}
    const std::vector<IndexLayout::Element>& get_primal_slave_indices() const {return m_vUniquePrimalSlaveIndex;}
 
  //  union of indices
    const std::vector<IndexLayout::Element>& get_primal_indices() const {return m_vUniquePrimalIndex;}
    const std::vector<IndexLayout::Element>& get_dual_indices() const {return m_vUniqueDualIndex;}
 
  public:
 
  //  sets standard communication layouts and communicators
    void vec_use_std_communication(vector_type& vec)
    {
      vec.set_layouts(m_spStdLayouts);
    }
 
  //  sets intra subdomain communication layouts and communicators
    void vec_use_intra_sd_communication(vector_type& vec)
    {
      vec.set_layouts(m_spInnerLayouts);
    }
 
  public:
    void vec_scale_add_on_dual(vector_type& vecDest,
                               number alpha1, const vector_type& vecSrc1,
                               number alpha2, const vector_type& vecSrc2)
    {
      VecScaleAdd(vecDest, alpha1, vecSrc1, alpha2, vecSrc2, m_vUniqueDualIndex);
    }
 
    void vec_scale_add_on_primal(vector_type& vecDest,
                               number alpha1, const vector_type& vecSrc1,
                               number alpha2, const vector_type& vecSrc2)
    {
      VecScaleAdd(vecDest, alpha1, vecSrc1, alpha2, vecSrc2, m_vUniquePrimalIndex);
    }
 
    void vec_scale_append_on_dual(vector_type& vecInOut,
                                  const vector_type& vecSrc1, number alpha1)
    {
      VecScaleAdd(vecInOut, alpha1, vecSrc1, 1.0, vecInOut, m_vUniqueDualIndex);
    }
 
    void vec_scale_append_on_primal(vector_type& vecInOut,
                                  const vector_type& vecSrc1, number alpha1)
    {
      VecScaleAdd(vecInOut, alpha1, vecSrc1, 1.0, vecInOut, m_vUniquePrimalIndex);
    }
 
    void vec_set_on_dual(vector_type& vecDest, number alpha)
    {
      VecSet(vecDest, alpha, m_vUniqueDualIndex);
    }
 
    void vec_set_on_primal(vector_type& vecDest, number alpha)
    {
      VecSet(vecDest, alpha, m_vUniquePrimalIndex);
    }
 
    void vec_scale_assign_on_dual(vector_type& vecDest, const vector_type& vecSrc, number alpha)
    {
      VecScaleAssign(vecDest, alpha, vecSrc, m_vUniqueDualIndex);
    }
 
    void vec_scale_assign_on_primal(vector_type& vecDest, const vector_type& vecSrc, number alpha)
    {
      VecScaleAssign(vecDest, alpha, vecSrc, m_vUniquePrimalIndex);
    }
 
    number vec_norm_on_identified_dual(vector_type& vec)
    {
    //  forward to VecProc
      return sqrt(vec_prod_on_identified_dual(vec, vec));
    }
 
    number vec_prod_on_identified_dual(const vector_type& vecSrc1, const vector_type& vecSrc2)
    {
    //  compute norm on dual master indices
      double tProdLocal = VecProd(vecSrc1, vecSrc2, m_vUniqueMasterDualIndex);
 
    //  global value
      double tProdGlobal;
 
    //  sum up values
      m_spStdLayouts->proc_comm().allreduce(&tProdLocal, &tProdGlobal, 1,
                      PCL_DT_DOUBLE, PCL_RO_SUM);
 
    //  return result
      return tProdGlobal;
    }
 
  public:
    void mat_set_dirichlet_on_dual(matrix_type& mat)
    {
      SetDirichletRow(mat, m_vUniqueDualIndex);
    }
 
    void mat_set_dirichlet_on_primal(matrix_type& mat)
    {
      SetDirichletRow(mat, m_vUniquePrimalIndex);
    }
 
  //  sets standard communication layouts and communicators
    void mat_use_std_communication(matrix_type& mat)
    {
      mat.set_slave_layout(get_std_slave_layout());
      mat.set_master_layout(get_std_master_layout());
      mat.set_process_communicator(get_std_process_communicator());
    }
 
  //  sets intra subdomain communication layouts and communicators
    void mat_use_intra_sd_communication(matrix_type& mat)
    {
      mat.set_layouts(m_spInnerLayouts);
    }
 
  protected:
    void sort_and_remove_doubles(std::vector<IndexLayout::Element>& vIndex)
    {
    //  sort vector
      std::sort(vIndex.begin(), vIndex.end());
 
    //  remove doubles
      vIndex.erase(std::unique(vIndex.begin(), vIndex.end()),
                               vIndex.end());
    }
 
    void add_merge_sort_remove_doubles(std::vector<IndexLayout::Element>& vOut,
                                       const std::vector<IndexLayout::Element>& v)
    {
    //  on entry, vOut and v are supposed to be sorted
 
    //  original size of vOut
      size_t mid = vOut.size();
 
    //  add second vector
      vOut.insert(vOut.end(), v.begin(), v.end() );
 
    //  merge both
      std::inplace_merge(vOut.begin(), vOut.begin()+mid, vOut.end());
 
    //  remove doubles
      vOut.erase(std::unique(vOut.begin(), vOut.end()),
                               vOut.end());
    }
 
  protected:
  //  Standard Layouts
    ConstSmartPtr<AlgebraLayouts> m_spStdLayouts;
 
  //  Layouts and Communicator for Inner variables
    SmartPtr<AlgebraLayouts> m_spInnerLayouts;
 
  //  Layouts for Dual variables
    IndexLayout m_masterDualLayout;
    IndexLayout m_slaveDualLayout;
    std::vector<IndexLayout::Element> m_vUniqueMasterDualIndex;
    std::vector<IndexLayout::Element> m_vUniqueSlaveDualIndex;
 
  //  Layouts for Dual Neighbour variables
    IndexLayout m_masterDualNbrLayout;
    IndexLayout m_slaveDualNbrLayout;
    std::vector<IndexLayout::Element> m_vUniqueMasterDualNbrIndex;
    std::vector<IndexLayout::Element> m_vUniqueSlaveDualNbrIndex;
 
  //  Layouts for Primal variables
    IndexLayout m_masterPrimalLayout;
    IndexLayout m_slavePrimalLayout;
    std::vector<IndexLayout::Element> m_vUniquePrimalMasterIndex;
    std::vector<IndexLayout::Element> m_vUniquePrimalSlaveIndex;
 
  //  unique indices for layouts
    std::vector<IndexLayout::Element> m_vUniqueDualIndex;
    std::vector<IndexLayout::Element> m_vUniquePrimalIndex;
 
}; /* end class 'FetiLayouts' */
 
 
template <typename TVector>
void ComputeDifferenceOnDelta(TVector& diff, const TVector& u,
                IndexLayout&    dualMasterLayoutIn,
                IndexLayout&     dualSlaveLayoutIn,
                IndexLayout& dualNbrMasterLayoutIn,
                IndexLayout&  dualNbrSlaveLayoutIn)
{
  // Copy all values
  diff = u;
 
  // Communicate values:
  // (a) Slaves send their u values, every master subtracts the value of only
  //     one of his slaves ...
  VecSubtractOneSlaveFromMaster(&diff, dualMasterLayoutIn, dualSlaveLayoutIn);
 
  // (b) masters send their diff values to all slaves (but not vice versa)
  VecCopy(&diff, dualMasterLayoutIn, dualSlaveLayoutIn);
 
  // (c) and also to the additional slaves living in "Dual neighbour layout"
  VecCopy(&diff, dualNbrMasterLayoutIn, dualNbrSlaveLayoutIn);
 
  // now the vector 'diff' should be consistent!
  return;
};
 
 
template <typename TVector>
void ComputeDifferenceOnDeltaTransposed(TVector& f, const TVector& diff,
                    const std::vector<IndexLayout::Element>& vDualMasterIndex,
                    const std::vector<IndexLayout::Element>& vDualSlaveIndex,
                    const std::vector<IndexLayout::Element>& vDualNbrSlaveIndex)
{
  // Copy values (no communication is performed):
  // (a) set f = +1 * d on masters ...
  VecScaleAssign(f, 1.0, diff, vDualMasterIndex);
 
  // (b) set f = -1 * d on slaves
  VecScaleAssign(f, -1.0, diff, vDualSlaveIndex);
 
  // (c) set f = +1 * d on slaves living in "Dual neighbour layout"
  //     ("+1" since in this context these dofs play the role of masters!)
  VecScaleAssign(f, +1.0, diff, vDualNbrSlaveIndex);
  return;
 
};
 
 
template <typename TAlgebra>
class LocalSchurComplement
  : public ILinearOperator< typename TAlgebra::vector_type,
                              typename TAlgebra::vector_type>,
    public DebugWritingObject<TAlgebra>
{
  public:
  //  Algebra type
    typedef TAlgebra algebra_type;
 
  //  Vector type
    typedef typename TAlgebra::vector_type vector_type;
 
  //  Matrix type
    typedef typename TAlgebra::matrix_type matrix_type;
 
  protected:
    using DebugWritingObject<TAlgebra>::write_debug;
    using DebugWritingObject<TAlgebra>::debug_writer;
 
  public:
    LocalSchurComplement();
 
    virtual const char* name() const {return "Local Schur Complement Solver";}
 
 
    void set_matrix(SmartPtr<MatrixOperator<matrix_type, vector_type> > A)
    {
    //  save current operator
      m_spOperator = A;
    }
 
 
    void set_dirichlet_solver(SmartPtr<ILinearOperatorInverse<vector_type> > dirichletSolver)
    {
    //  remember the Dirichlet Solver
      m_spDirichletSolver = dirichletSolver;
    }
 
    void set_feti_layouts(FetiLayouts<algebra_type>& fetiLayouts)
    {
      m_pFetiLayouts = &fetiLayouts;
    }
 
    void init(const vector_type& u) {init();}
 
 
    virtual void init();
 
    virtual void apply(vector_type& f, const vector_type& u);
 
    virtual void apply_sub(vector_type& f, const vector_type& u);
 
    void set_statistic_type(std::string type) {m_statType = type;}
 
    void print_statistic_of_inner_solver(bool bPrintOnlyAverages); // const;
 
    void clear_total_itercnt_of_inner_solvers() {m_totalIterCntOfInnerSolvers = 0;}
    int get_total_itercnt_of_inner_solvers() {return m_totalIterCntOfInnerSolvers;}
 
  // destructor
    virtual ~LocalSchurComplement() {};
 
  protected:
  //  Operator that is inverted by this Inverse Operator
    SmartPtr<MatrixOperator<matrix_type,vector_type> > m_spOperator;
 
  //  Parallel Matrix
    matrix_type* m_pMatrix;
 
  //  Feti Layouts
    FetiLayouts<algebra_type>* m_pFetiLayouts;
 
  //  Copy of matrix
    SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spDirichletOperator;
 
  //  Parallel Dirichlet Matrix
    matrix_type* m_pDirichletMatrix;
 
  //  Linear Solver to invert the local Dirichlet problems
    SmartPtr<ILinearOperatorInverse<vector_type> > m_spDirichletSolver;
 
  //  Convergence history
    std::string m_statType;
 
    struct StepConv
    {
      int numIter3b;
      number lastDef3b;
    };
 
    std::map<std::string, std::vector<StepConv> > m_mvStepConv;
 
    int m_applyCnt;
 
    int m_totalIterCntOfInnerSolvers;
 
}; /* end class 'LocalSchurComplement' */
 
/* 1.7 Application of \f${\tilde{S}_{\Delta \Delta}}^{-1}\f$ */ 
 
template <typename TAlgebra>
class PrimalSubassembledMatrixInverse
  : public ILinearOperatorInverse<  typename TAlgebra::vector_type,
                                  typename TAlgebra::vector_type>,
    public DebugWritingObject<TAlgebra>
{
  public:
  //  Algebra type
    typedef TAlgebra algebra_type;
 
  //  Vector type
    typedef typename TAlgebra::vector_type vector_type;
 
  //  Matrix type
    typedef typename TAlgebra::matrix_type matrix_type;
 
    typedef ILinearOperatorInverse<vector_type> base_type;
 
  protected:
    using base_type::convergence_check;
    using DebugWritingObject<TAlgebra>::write_debug;
    using DebugWritingObject<TAlgebra>::debug_writer;
 
  public:
    PrimalSubassembledMatrixInverse();
 
    virtual const char* name() const {return "Schur Complement Inverse";}
 
    virtual bool supports_parallel() const
    {
      bool bRet = true;
      if(m_spNeumannSolver.valid() || !m_spNeumannSolver->supports_parallel())
        bRet = false;
      if(m_spCoarseProblemSolver.valid() || !m_spCoarseProblemSolver->supports_parallel())
        bRet = false;
      return bRet;
    }
 
    void set_neumann_solver(SmartPtr<ILinearOperatorInverse<vector_type> > neumannSolver)
    {
    //  remember the Dirichlet Solver
      m_spNeumannSolver = neumannSolver;
    }
 
    void set_coarse_problem_solver(SmartPtr<ILinearOperatorInverse<vector_type> > coarseProblemSolver)
    {
    //  remember the coarse problem Solver
      m_spCoarseProblemSolver = coarseProblemSolver;
    }
 
    void set_feti_layouts(FetiLayouts<algebra_type>& fetiLayouts)
    {
      m_pFetiLayouts = &fetiLayouts;
    }
 
  //  Init for Linear Operator L
    virtual bool init(SmartPtr<ILinearOperator<vector_type> > L);
 
 
  //  Init for Linear Operator J and Linearization point (current solution)
    virtual bool init(SmartPtr<ILinearOperator<vector_type> > J, const vector_type& u)
    {
      return init(J);
    }
 
  //  Solve A*u = f, such that u = A^{-1} f
    virtual bool apply(vector_type& u, const vector_type& f);
 
  //  Solve A*u = f, such that u = A^{-1} f
  //  This is done by iterating: u := u + B(f - A*u)
  //  In f the last defect f := f - A*u is returned
    virtual bool apply_return_defect(vector_type& u, vector_type& f);
 
    void set_statistic_type(std::string type) {m_statType = type;}
 
    void print_statistic_of_inner_solver(bool bPrintOnlyAverages); // const;
 
    void clear_total_itercnt_of_inner_solvers() {m_totalIterCntOfInnerSolvers = 0;}
    int get_total_itercnt_of_inner_solvers() {return m_totalIterCntOfInnerSolvers;}
 
    void set_test_one_to_many_layouts(bool bTest) {m_bTestOneToManyLayouts = bTest;}
 
  //  destructor
    virtual ~PrimalSubassembledMatrixInverse() {};
 
  protected:
  //  Operator that is inverted by this Inverse Operator ==> from which SC is built (05022011)
    SmartPtr<MatrixOperator<matrix_type,vector_type> > m_spOperator;
 
  //  Parallel Matrix to invert ==> from which SC is built (05022011)
    matrix_type* m_pMatrix;
 
  //  Feti Layouts
    FetiLayouts<algebra_type>* m_pFetiLayouts;
 
  //  Copy of matrix
    SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spNeumannOperator;
 
  //  Parallel Neumann Matrix
    matrix_type* m_pNeumannMatrix;
 
  //  Neumann Solver
    SmartPtr<ILinearOperatorInverse<vector_type> > m_spNeumannSolver;
 
  //  Coarse problem Solver on root.
    SmartPtr<ILinearOperatorInverse<vector_type> > m_spCoarseProblemSolver;
 
  //  Process gathering the Schur Complement w.r.t. Primal unknowns
    int m_primalRootProc;
 
  //  Layouts for all to one communication
    IndexLayout m_masterAllToOneLayout;
    IndexLayout m_slaveAllToOneLayout;
    pcl::ProcessCommunicator m_allToOneProcessComm;
 
  //  Schur Complement operator for gathered matrix
    SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spRootSchurComplementOp;
 
  //  Matrix for one proc Schur complement
    matrix_type* m_pRootSchurComplementMatrix;
 
  //  Convergence history
    std::string m_statType;
 
    struct StepConv
    {
      int numIter2a;
      int numIter7;
      int numIterSC;
      number lastDef2a;
      number lastDef7;
      number lastDefSC;
    };
 
    std::map<std::string, std::vector<StepConv> > m_mvStepConv;
 
  //  testing of 'one-to-many layouts'
    bool m_bTestOneToManyLayouts;
 
    int m_totalIterCntOfInnerSolvers;
 
}; /* end class 'PrimalSubassembledMatrixInverse' */
 
 
template <typename TAlgebra>
class FETISolver : public IMatrixOperatorInverse< typename TAlgebra::matrix_type,
                          typename TAlgebra::vector_type>,
  public DebugWritingObject<TAlgebra>
{
  public:
  //  Algebra type
    typedef TAlgebra algebra_type;
 
  //  Vector type
    typedef typename TAlgebra::vector_type vector_type;
 
  //  Matrix type
    typedef typename TAlgebra::matrix_type matrix_type;
 
    typedef ILinearOperatorInverse<vector_type> base_type;
 
  protected:
    using base_type::convergence_check;
    using DebugWritingObject<TAlgebra>::write_debug;
    using DebugWritingObject<TAlgebra>::debug_writer;
 
  public:
    FETISolver();
 
    virtual const char* name() const {return "FETI Solver";}
 
    virtual bool supports_parallel() const
    {
      if(m_spDirichletSolver.valid())
        if(!m_spDirichletSolver->supports_parallel())
          return false;
      if(m_spNeumannSolver.valid())
        if(!m_spNeumannSolver->supports_parallel())
          return false;
      if(m_spCoarseProblemSolver.valid())
        if(!m_spCoarseProblemSolver->supports_parallel())
          return false;
      return true;
    }
 
    void set_dirichlet_solver(SmartPtr<ILinearOperatorInverse<vector_type> > dirichletSolver)
    {
    //  remember the Dirichlet Solver
      m_spDirichletSolver = dirichletSolver;
    }
 
    void set_neumann_solver(SmartPtr<ILinearOperatorInverse<vector_type> > neumannSolver)
    {
    //  remember the Dirichlet Solver
      m_spNeumannSolver = neumannSolver;
    }
 
    void set_coarse_problem_solver(SmartPtr<ILinearOperatorInverse<vector_type> > coarseProblemSolver)
    {
    //  remember the coarse problem Solver
      m_spCoarseProblemSolver = coarseProblemSolver;
    }
 
  //  set debug output
    void set_debug(SmartPtr<IDebugWriter<algebra_type> > spDebugWriter)
    {
      m_LocalSchurComplement.set_debug(spDebugWriter);
      m_PrimalSubassembledMatrixInverse.set_debug(spDebugWriter);
      VectorDebugWritingObject<vector_type>::set_debug(spDebugWriter);
    }
 
    virtual bool init(SmartPtr<MatrixOperator<matrix_type, vector_type> > A);
 
 
    bool apply_F(vector_type& f, const vector_type& v);
 
 
    bool compute_d(vector_type& d, const vector_type& f);
 
    bool apply_scaling_matrix(vector_type& s, const vector_type& v) // maybe restrict to layout
    {
    //  copy values
      s = v;
 
    //  scale by 1/2
      s *= 1./2.;
 
      // \todo: more general: m_Ddelta.apply(s,v), with additional member 'matrix_type m_Ddelta;'
 
    //  we're done
      return true;
    }
 
 
    bool apply_M_inverse(vector_type& z, const vector_type& r);
 
  //  tests layouts:
    void test_layouts(bool bPrint);
 
  //  set member 'm_bTestOneToManyLayouts' of class 'PrimalSubassembledMatrixInverse':
    void set_test_one_to_many_layouts(bool bTest)
    {
      m_PrimalSubassembledMatrixInverse.set_test_one_to_many_layouts(bTest);
    }
 
 
    virtual bool apply_return_defect(vector_type& lambda, vector_type& d);
 
    virtual bool apply(vector_type& x, const vector_type& b)
    {
    //  copy defect
      vector_type d; d.resize(b.size());
      d = b;
 
    //  solve on copy of defect
      return apply_return_defect(x, d);
    }
 
    void print_statistic_of_inner_solver(bool bPrintOnlyAverages) //const
    {
      // sum over all calls of inner solvers is done in the following print functions!
      m_PrimalSubassembledMatrixInverse.clear_total_itercnt_of_inner_solvers();
      m_PrimalSubassembledMatrixInverse.print_statistic_of_inner_solver(bPrintOnlyAverages);
      UG_LOG("\n");
 
      m_LocalSchurComplement.clear_total_itercnt_of_inner_solvers();
      m_LocalSchurComplement.print_statistic_of_inner_solver(bPrintOnlyAverages);
 
      int totalIterCntOfInnerSolvers = m_PrimalSubassembledMatrixInverse.get_total_itercnt_of_inner_solvers();
      totalIterCntOfInnerSolvers     += m_LocalSchurComplement.get_total_itercnt_of_inner_solvers();
 
      UG_LOG("Total number of calls of sub problem solvers: " << totalIterCntOfInnerSolvers);
      UG_LOG(" in " << std::setw(5) << m_iterCnt << " FETI iterations ");
      UG_LOG(" on " << std::setw(5) << m_pDDInfo->get_num_subdomains() << " FETI subdomains, "
                  << std::fixed << (double)totalIterCntOfInnerSolvers/m_iterCnt
                  << std::scientific << " per FETI iteration, "
                  << std::fixed << (double)totalIterCntOfInnerSolvers/m_pDDInfo->get_num_subdomains()
                  << std::scientific << " per FETI subdomain.\n");
 
    }
 
    // destructor
    virtual ~FETISolver() {};
 
  protected:
  //  Prepare the convergence check
    void prepare_conv_check()
    {
      convergence_check()->set_name(name());
      convergence_check()->set_symbol('%');
 
      //  set preconditioner string
      std::stringstream ss;
      ss << " (Inherent Preconditioner) ";
      convergence_check()->set_info(ss.str());
    }
 
  protected:
    virtual void write_debug(const vector_type& vec, const char* filename)
    {
    //  add iter count to name
      std::string name(filename);
      char ext[20]; snprintf(ext, 20, "_iter%03d.vec", m_iterCnt);
      name.append(ext);
 
    //  if no debug writer set, we're done
      if(debug_writer() == SPNULL) return;
 
    //  write
      debug_writer()->write_vector(vec, name.c_str());
    }
 
    int m_iterCnt;
 
  protected:
  //  Operator that is inverted by this Inverse Operator
    SmartPtr<MatrixOperator<matrix_type,vector_type> > m_spOperator;
 
  //  Parallel Matrix to invert
    matrix_type* m_pMatrix;
 
  //  Feti Layouts
    FetiLayouts<algebra_type> m_fetiLayouts;
 
  //  Local Schur complement for each feti subdomain
    LocalSchurComplement<algebra_type> m_LocalSchurComplement;
 
  //  Dirichlet solver for inverse of \f$A_{II}\f$ in local Schur complement
    SmartPtr<ILinearOperatorInverse<vector_type> > m_spDirichletSolver;
 
  //  PrimalSubassembledMatrixInverse
    PrimalSubassembledMatrixInverse<algebra_type> m_PrimalSubassembledMatrixInverse;
 
  //  Neumann solver for inverse of \f$A_{\{I,\Delta\}, \{I,\Delta\}}\f$ in the
  //  creation of the S_{\Pi \Pi} Schur complement in PrimalSubassembledMatrixInverse
    SmartPtr<ILinearOperatorInverse<vector_type> > m_spNeumannSolver;
 
  //  Solver used in solving coarse problem on root.
  //  It solves \f$S_{\Pi \Pi} u_{\Pi} = \tilde{f}_{\Pi}\f$ 
    SmartPtr<ILinearOperatorInverse<vector_type> > m_spCoarseProblemSolver;
 
  public:
    void set_domain_decomp_info(pcl::IDomainDecompositionInfo& ddInfo)
    {
      m_pDDInfo = &ddInfo;
    }
 
  private:
  //  pointer to Domain decomposition info object
    pcl::IDomainDecompositionInfo* m_pDDInfo;
 
}; /* end class 'FETISolver' */
 
} // end namespace ug
 
#endif /* UG_PARALLEL */
 
#endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__FETI__ */