qpmg.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2017-2020 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef QPMG_HH
#define QPMG_HH
 
#include <memory>
#include <tuple>
 
#include "dune/istl/bvector.hh"
 
#include "linalg/dynamicMatrix.hh"
#include "linalg/threadedMatrix.hh"
#include "mg/prolongation.hh"
#include "utilities/enums.hh"
 
namespace Kaskade
{
 
  template <int d, class Real=double>
  class MGSolverStatistics
  {
  public:
    using VectorX = Dune::BlockVector<Dune::FieldVector<Real,d>>;
    using VectorB = Dune::BlockVector<Dune::FieldVector<Real,1>>;
 
    ~MGSolverStatistics() {}
 
    virtual void enterPreSmoothingCorrection(int level, VectorX const& dx) {}
    virtual void enterCoarseGridCorrection(int level, VectorX const& dx) {}
    virtual void enterPostSmoothingCorrection(int level, VectorX const& dx) {}
 
    virtual void enterVCycleCorrection(int iter, VectorX const& dx, Real energy, 
                                      VectorX const& grad, VectorX const& grad0, VectorX const& grad1,
                                      VectorB const& cons, int activeSetChanges, 
                                      std::vector<VectorX> &correctionStack,
                                      std::vector<std::tuple<double,double,double>> &coarseRes) 
                                      {}
  };
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
 
 
  template <int d, class Prolongation, class Smoother, class CoarseSmoother=Smoother, class Real=double>
  class QPMultiGridBase
  {
    using EntryA = Dune::FieldMatrix<Real,d,d>;
  public:
    using VectorX = Dune::BlockVector<Dune::FieldVector<Real,d>>;
    using VectorB = Dune::BlockVector<Dune::FieldVector<Real,1>>;
    using MatrixA = NumaBCRSMatrix<EntryA>;
    using MatrixB = NumaBCRSMatrix<Dune::FieldMatrix<Real,1,d>>;
    using Self = QPMultiGridBase<d,Prolongation,Smoother,CoarseSmoother,Real>;
    
    // check if coarse smoother is of equal type as rest of smoothers
    // if so, SmootherType = Smoother, else use std::variant<Smoother,CoarseSmoother>
    static constexpr bool smoothersEqual = std::is_same<Smoother,CoarseSmoother>::value;
    using SmootherType = std::conditional_t<smoothersEqual,
                                            Smoother,
                                            std::variant<Smoother,CoarseSmoother>>;
 
    QPMultiGridBase(MatrixA A, MatrixB const& B, std::vector<Prolongation>&& prolongations,
                    Real smootherRegularization=0, bool blocks=true, bool directOnCoarse=true);
 
    virtual ~QPMultiGridBase();
 
    VectorX step(VectorX const& x, VectorX c, VectorB b) const;
 
    std::tuple<VectorX,int> solve(VectorX x, VectorX const& c, VectorB const& b, double tol, int vcycles) const;
 
 
    Self& setCoarseLevel(int coarselevel);
 
    Self& setSmoothings(int pre, int post);
 
    Self& setCoarseCorrections(int n);
 
    Self& setBulkMode(ParallelMode m);
 
    Self& setLogger(MGSolverStatistics<d,Real>* logger);
 
  protected:
    virtual double computeEnergy(MatrixA const& A, MatrixB const& B,
                                 VectorX const& c, VectorB const& b, VectorX const& x) const = 0;
    virtual std::tuple<VectorX,VectorB,int> gradient(MatrixA const& A, MatrixB const& B,
                                                     VectorX const& c, VectorB const& b, VectorX const& x) const = 0;
                                                     
    virtual std::tuple<VectorX,VectorX,VectorX,VectorB,int> gradient_extended(MatrixA const& A, MatrixB const& B,
                                                                              VectorX const& c, VectorB const& b, VectorX const& x) const = 0;
                                                         
    virtual double qpLinesearch(MatrixA const& A, MatrixB const& B, VectorX const& c,
                                VectorB const& b, VectorX const& dx) const = 0;
 
    
    // vector of smoothers, one for each level
    // note: SmootherType might be either of Type
    // Smoother or std::variant<Smoother,CoarseSmoother>
   std::vector<SmootherType> smoothers;
 
    MatrixB const& B() const;
 
  private:
    VectorX mgRecursion(int level, int lpre, int lpost, VectorX c, VectorB b) const;
 
    MultiGridStack<Prolongation,EntryA,size_t> mgStackA;
    std::vector<MatrixB> mgStackB;
    mutable std::vector<VectorX> correctionStack;                                 // vector with coarse grid corrections on each level (interpolated to fine grid)
    mutable std::vector<std::tuple<double,double,double>> coarseResidual;         // vector with norm of residual on each level (after computing correction)
 
    mutable Real gamma;
    int coarseLevel = 0;
    int pre = 3;
    int post = 3;
    int mid = 1;
 
    std::unique_ptr<MGSolverStatistics<d,Real>> dummyLogger; // this dummy will be used if no other logger is given
    MGSolverStatistics<d,Real>* logger;
  };
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  // forward declaration
  template <int d, class Real, QPStructure sparsity=QPStructure::DENSE>
  class QPPenaltySmoother;
 
 
 
  template <int d, class Prolongation, class Real=double, class Smoother=QPPenaltySmoother<d,Real>, class CoarseSmoother=QPPenaltySmoother<d,Real>>
  class QPMultiGrid
  : public QPMultiGridBase<d,Prolongation,Smoother,CoarseSmoother,Real>
  {
    using EntryA = Dune::FieldMatrix<Real,d,d>;
    using Base = QPMultiGridBase<d,Prolongation,Smoother,CoarseSmoother,Real>;
 
  public:
    using VectorX = Dune::BlockVector<Dune::FieldVector<Real,d>>;
    using VectorB = Dune::BlockVector<Dune::FieldVector<Real,1>>;
    using MatrixA = NumaBCRSMatrix<EntryA>;
    using MatrixB = NumaBCRSMatrix<Dune::FieldMatrix<Real,1,d>>;
    using Self = QPMultiGrid<d,Prolongation,Real,Smoother,CoarseSmoother>;
 
    QPMultiGrid(MatrixA const& A, MatrixB const& B, std::vector<Prolongation>&& prolongations,
                 Real smootherRegularization=0, bool blocks=true, bool directOnCoarse=true);
 
    ~QPMultiGrid();
 
    std::tuple<VectorX,VectorB,int> solve(VectorX x, VectorX const& c, VectorB const& b, VectorB const& lambda, double tol, int vcycles) const;
 
//     /**
//      * \brief A least-squares multiplier update.
//      *
//      * This can be used if a good initial guess for \f$ x \f$ is available, but not for \f$ \lambda \f$.
//      * A common situation where this is typically the case is finite strain contact mechanics, where the
//      * constraints and hence the multiplier changes with deformation, after having solved a QP linearization.
//      */
//     void leastSquaresMultiplierUpdate(VectorX const& c, VectorB const& b, VectorX const& x, VectorB& lambda) const;
//
 
    VectorB firstOrderUpdate(VectorX const& x, VectorB const& b, VectorB const& lambda) const;
 
 
    Self& setPenalty(Real gamma, bool levelDependent=false);
 
    Self& setConstraintsMode(ParallelMode m);
    Self& setGlobalMode(ParallelMode m);
    
    double getPenalty();
 
  protected:
    virtual double computeEnergy(MatrixA const& A, MatrixB const& B,
                                 VectorX const& c, VectorB const& b, VectorX const& x) const;
    virtual std::tuple<VectorX,VectorB,int> gradient(MatrixA const& A, MatrixB const& B,
                                                     VectorX const& c, VectorB const& b, VectorX const& x) const;
    virtual std::tuple<VectorX,VectorX,VectorX,VectorB,int> gradient_extended(MatrixA const& A, MatrixB const& B,
                                                                              VectorX const& c, VectorB const& b, VectorX const& x) const;
    
    virtual double qpLinesearch(MatrixA const& A, MatrixB const& B, VectorX const& c,
                                VectorB const& b, VectorX const& dx) const;
 
 
  private:
    mutable Real gamma;
  };
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  // forward declaration
  template <int d, class Real, QPStructure sparsity=QPStructure::DENSE>
  class QPSmoother;
 
  template <int d, class Prolongation, class Real=double, class Smoother=QPSmoother<d,Real>, class CoarseSmoother=QPSmoother<d,Real>>
  class QPMultiGridStrict
  : public QPMultiGridBase<d,Prolongation,Smoother,CoarseSmoother,Real>
  {
    using Base = QPMultiGridBase<d,Prolongation,Smoother,CoarseSmoother,Real>;
 
  public:
    using Self = QPMultiGridStrict<d,Prolongation,Real,Smoother,CoarseSmoother>;
    using typename Base::MatrixA;
    using typename Base::MatrixB;
    using typename Base::VectorX;
    using typename Base::VectorB;
 
    QPMultiGridStrict(MatrixA const& A, MatrixB const& B, std::vector<Prolongation>&& prolongations,
                      Real smootherRegularization=0, bool blocks=true, bool directOnCoarse=true);
 
    ~QPMultiGridStrict();
    
    
    std::tuple<VectorX,int> solve(VectorX x, VectorX const& c, VectorB const& b, double tol, int vcycles) const;
 
 
  protected:
    virtual double computeEnergy(MatrixA const& A, MatrixB const& B,
                                 VectorX const& c, VectorB const& b, VectorX const& x) const;
                                 
    virtual std::tuple<VectorX,VectorB,int> gradient(MatrixA const& A, MatrixB const& B,
                                                     VectorX const& c, VectorB const& b, VectorX const& x) const;
                                                     
                                                     
    virtual std::tuple<VectorX,VectorX,VectorX,VectorB,int> gradient_extended(MatrixA const& A, MatrixB const& B,
                                                                              VectorX const& c, VectorB const& b, VectorX const& x) const;
 
    virtual double qpLinesearch(MatrixA const& A, MatrixB const& B, VectorX const& c,
                                VectorB const& b, VectorX const& dx) const;
  };
}
 
#endif