base.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/util/exception.hpp>
#include <kernel/util/math.hpp>
#include <kernel/util/property_map.hpp>
#include <kernel/util/string.hpp>
 
// includes, system
#include <memory>
 
namespace FEAT
{
  namespace Solver
  {
 
    namespace Intern
    {
    } // namespace Intern
 
    enum class Status : std::uint8_t
    {
      undefined = 0u,
      progress = 1u,
      success = 2u,
      aborted = 3u,
      diverged = 4u,
      max_iter = 5u,
      stagnated = 6u
    };
 
    inline std::ostream& operator<<(std::ostream& os, Status status)
    {
      switch(status)
      {
      case Status::undefined:
        return os << "undefined";
      case Status::progress:
        return os << "progress";
      case Status::success:
        return os << "success";
      case Status::aborted:
        return os << "aborted";
      case Status::diverged:
        return os << "diverged";
      case Status::max_iter:
        return os << "max-iter";
      case Status::stagnated:
        return os << "stagnated";
      default:
        return os << "-unknown-";
      }
    }
 
    inline bool status_success(Status status)
    {
      switch(status)
      {
      case Status::success:
      case Status::max_iter:
      case Status::stagnated:
        return true;
 
      default:
        return false;
      }
    }
 
    class SolverException :
      public Exception
    {
    public:
      explicit SolverException(const String& msg) :
        Exception(msg)
      {
      }
    };
 
    class InvalidMatrixStructureException :
      public SolverException
    {
    public:
      explicit InvalidMatrixStructureException(const String& msg = "invalid matrix structure") :
        SolverException(msg)
      {
      }
    };
 
    class SingularMatrixException :
      public SolverException
    {
    public:
      explicit SingularMatrixException(const String& msg = "singular matrix") :
        SolverException(msg)
      {
      }
    };
 
    template<typename Vector_>
    class SolverBase
    {
    public:
      typedef Vector_ VectorType;
 
    public:
      SolverBase()
      {
      }
 
      explicit SolverBase(const String& DOXY(section_name), const PropertyMap* DOXY(config_section))
      {
      }
 
      virtual ~SolverBase()
      {
      }
 
      virtual void init_symbolic()
      {
      }
 
      virtual void init_numeric()
      {
      }
 
      virtual void done_numeric()
      {
      }
 
      virtual void done_symbolic()
      {
      }
 
      virtual void init()
      {
        init_symbolic();
        init_numeric();
      }
 
      virtual void done()
      {
        done_numeric();
        done_symbolic();
      }
 
      virtual String name() const = 0;
 
      virtual Status apply(Vector_& vec_cor, const Vector_& vec_def) = 0;
    }; // class SolverBase<...>
 
    template<
      typename Vector_,
      typename Matrix_,
      typename Filter_>
    inline Status solve(
      SolverBase<Vector_>& solver,
      Vector_& vec_sol,
      const Vector_& vec_rhs,
      const Matrix_& matrix,
      const Filter_& filter)
    {
      typedef typename Vector_::DataType DataType;
 
      // create two temporary vectors
      Vector_ vec_def(vec_rhs.clone(LAFEM::CloneMode::Layout));
      Vector_ vec_cor(vec_sol.clone(LAFEM::CloneMode::Layout));
 
      // compute defect
      matrix.apply(vec_def, vec_sol, vec_rhs, -DataType(1));
 
      // apply defect filter
      filter.filter_def(vec_def);
 
      // apply solver
      Status status = solver.apply(vec_cor, vec_def);
      if(status_success(status))
      {
        // apply correction filter
        filter.filter_cor(vec_cor);
 
        // update solution
        vec_sol.axpy(vec_cor);
      }
 
      // return status
      return status;
    }
 
    class IterationStats
    {
    private:
      const String _solver_name;
      TimeStamp _at;
      double _mpi_execute_reduction_start;
      double _mpi_execute_reduction_stop;
      double _mpi_execute_blas2_start;
      double _mpi_execute_blas2_stop;
      double _mpi_execute_blas3_start;
      double _mpi_execute_blas3_stop;
      double _mpi_execute_collective_start;
      double _mpi_execute_collective_stop;
      double _mpi_wait_start_reduction;
      double _mpi_wait_start_blas2;
      double _mpi_wait_start_blas3;
      double _mpi_wait_start_collective;
      double _mpi_wait_stop_reduction;
      double _mpi_wait_stop_blas2;
      double _mpi_wait_stop_blas3;
      double _mpi_wait_stop_collective;
      bool _destroyed;
 
    public:
      template<typename Vector_>
      explicit IterationStats(const SolverBase<Vector_>& solver) :
        _solver_name(solver.name()),
        _destroyed(false)
      {
        _mpi_execute_reduction_start = Statistics::get_time_mpi_execute_reduction();
        _mpi_execute_blas2_start = Statistics::get_time_mpi_execute_blas2();
        _mpi_execute_blas3_start = Statistics::get_time_mpi_execute_blas3();
        _mpi_execute_collective_start = Statistics::get_time_mpi_execute_collective();
        _mpi_wait_start_reduction    = Statistics::get_time_mpi_wait_reduction();
        _mpi_wait_start_blas2    = Statistics::get_time_mpi_wait_blas2();
        _mpi_wait_start_blas3    = Statistics::get_time_mpi_wait_blas3();
        _mpi_wait_start_collective    = Statistics::get_time_mpi_wait_collective();
        _mpi_execute_reduction_stop = _mpi_execute_reduction_start;
        _mpi_execute_blas2_stop = _mpi_execute_blas2_start;
        _mpi_execute_blas3_stop = _mpi_execute_blas3_start;
        _mpi_execute_collective_stop = _mpi_execute_collective_start;
        _mpi_wait_stop_reduction    = _mpi_wait_start_reduction;
        _mpi_wait_stop_blas2    = _mpi_wait_start_blas2;
        _mpi_wait_stop_blas3    = _mpi_wait_start_blas3;
        _mpi_wait_stop_collective    = _mpi_wait_start_collective;
      }
 
      // delete copy-ctor and assign operator
      IterationStats(const IterationStats&) = delete;
      IterationStats& operator=(const IterationStats&) = delete;
 
      ~IterationStats()
      {
        if (!_destroyed)
          destroy();
      }
 
      void destroy()
      {
        XASSERTM(!_destroyed, "IterationStats::destroy() was already called before!");
 
        _mpi_execute_reduction_stop = Statistics::get_time_mpi_execute_reduction();
        _mpi_execute_blas2_stop = Statistics::get_time_mpi_execute_blas2();
        _mpi_execute_blas3_stop = Statistics::get_time_mpi_execute_blas3();
        _mpi_execute_collective_stop = Statistics::get_time_mpi_execute_collective();
        _mpi_wait_stop_reduction    = Statistics::get_time_mpi_wait_reduction();
        _mpi_wait_stop_blas2    = Statistics::get_time_mpi_wait_blas2();
        _mpi_wait_stop_blas3    = Statistics::get_time_mpi_wait_blas3();
        _mpi_wait_stop_collective    = Statistics::get_time_mpi_wait_collective();
        Statistics::add_solver_expression(std::make_shared<ExpressionTimings>(_solver_name, _at.elapsed_now(),
          _mpi_execute_reduction_stop - _mpi_execute_reduction_start,
          _mpi_execute_blas2_stop - _mpi_execute_blas2_start,
          _mpi_execute_blas3_stop - _mpi_execute_blas3_start,
          _mpi_execute_collective_stop - _mpi_execute_collective_start,
          _mpi_wait_stop_reduction - _mpi_wait_start_reduction,
          _mpi_wait_stop_blas2 - _mpi_wait_start_blas2,
          _mpi_wait_stop_blas3 - _mpi_wait_start_blas3,
          _mpi_wait_stop_collective - _mpi_wait_start_collective));
 
        _destroyed = true;
      }
    }; // class IterationStats
  } // namespace Solver
} // namespace FEAT