117 lines
5.4 KiB
C++
117 lines
5.4 KiB
C++
/*
|
||
* Copyright © 2004-2011 Ondra Kamenik
|
||
* Copyright © 2019-2023 Dynare Team
|
||
*
|
||
* This file is part of Dynare.
|
||
*
|
||
* Dynare is free software: you can redistribute it and/or modify
|
||
* it under the terms of the GNU General Public License as published by
|
||
* the Free Software Foundation, either version 3 of the License, or
|
||
* (at your option) any later version.
|
||
*
|
||
* Dynare 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 General Public License for more details.
|
||
*
|
||
* You should have received a copy of the GNU General Public License
|
||
* along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||
*/
|
||
|
||
#ifndef TRIANGULAR_SYLVESTER_HH
|
||
#define TRIANGULAR_SYLVESTER_HH
|
||
|
||
#include "KronVector.hh"
|
||
#include "QuasiTriangular.hh"
|
||
#include "QuasiTriangularZero.hh"
|
||
#include "SimilarityDecomp.hh"
|
||
#include "SylvesterSolver.hh"
|
||
|
||
#include <memory>
|
||
|
||
class TriangularSylvester : public SylvesterSolver
|
||
{
|
||
const std::unique_ptr<const QuasiTriangular> matrixKK;
|
||
const std::unique_ptr<const QuasiTriangular> matrixFF;
|
||
|
||
public:
|
||
TriangularSylvester(const QuasiTriangular& k, const QuasiTriangular& f);
|
||
TriangularSylvester(const SchurDecompZero& kdecomp, const SchurDecomp& fdecomp);
|
||
TriangularSylvester(const SchurDecompZero& kdecomp, const SimilarityDecomp& fdecomp);
|
||
|
||
~TriangularSylvester() override = default;
|
||
void print() const;
|
||
void solve(SylvParams& pars, KronVector& d) const override;
|
||
|
||
void solvi(double r, KronVector& d, double& eig_min) const;
|
||
void solvii(double alpha, double beta1, double beta2, KronVector& d1, KronVector& d2,
|
||
double& eig_min) const;
|
||
void solviip(double alpha, double betas, KronVector& d, double& eig_min) const;
|
||
/* Computes:
|
||
⎛x₁⎞ ⎛d₁⎞ ⎛ α −β₁⎞ ⎛d₁⎞
|
||
⎢ ⎥=⎢ ⎥+⎢ ⎥⊗Fᵀ⊗Fᵀ⊗…⊗K·⎢ ⎥
|
||
⎝x₂⎠ ⎝d₂⎠ ⎝−β₂ α ⎠ ⎝d₂⎠
|
||
*/
|
||
void linEval(double alpha, double beta1, double beta2, KronVector& x1, KronVector& x2,
|
||
const ConstKronVector& d1, const ConstKronVector& d2) const;
|
||
void
|
||
linEval(double alpha, double beta1, double beta2, KronVector& x1, KronVector& x2,
|
||
const KronVector& d1, const KronVector& d2) const
|
||
{
|
||
linEval(alpha, beta1, beta2, x1, x2, ConstKronVector(d1), ConstKronVector(d2));
|
||
}
|
||
|
||
/* Computes:
|
||
⎛x₁⎞ ⎛d₁⎞ ⎛γ −δ₁⎞ ⎛d₁⎞ ⎛γ −δ₁⎞² ⎛d₁⎞
|
||
⎢ ⎥=⎢ ⎥+2α⎢ ⎥⊗Fᵀ⊗Fᵀ⊗…⊗K·⎢ ⎥+(α²+β)⎢ ⎥ ⊗Fᵀ²⊗Fᵀ²⊗…⊗K²·⎢ ⎥
|
||
⎝x₂⎠ ⎝d₂⎠ ⎝δ₂ γ ⎠ ⎝d₂⎠ ⎝δ₂ γ ⎠ ⎝d₂⎠
|
||
*/
|
||
void quaEval(double alpha, double betas, double gamma, double delta1, double delta2,
|
||
KronVector& x1, KronVector& x2, const ConstKronVector& d1,
|
||
const ConstKronVector& d2) const;
|
||
void
|
||
quaEval(double alpha, double betas, double gamma, double delta1, double delta2, KronVector& x1,
|
||
KronVector& x2, const KronVector& d1, const KronVector& d2) const
|
||
{
|
||
quaEval(alpha, betas, gamma, delta1, delta2, x1, x2, ConstKronVector(d1), ConstKronVector(d2));
|
||
}
|
||
|
||
private:
|
||
/* Returns square of size of minimal eigenvalue of the system solved,
|
||
now obsolete */
|
||
[[nodiscard]] double getEigSep(int depth) const;
|
||
// Recursively calculates kronecker product of complex vectors (used in getEigSep)
|
||
static void multEigVector(KronVector& eig, const Vector& feig, const Vector& keig);
|
||
|
||
using const_diag_iter = QuasiTriangular::const_diag_iter;
|
||
using const_row_iter = QuasiTriangular::const_row_iter;
|
||
|
||
// Called from solvi
|
||
void solviRealAndEliminate(double r, const_diag_iter di, KronVector& d, double& eig_min) const;
|
||
void solviComplexAndEliminate(double r, const_diag_iter di, KronVector& d, double& eig_min) const;
|
||
// Called from solviip
|
||
void solviipRealAndEliminate(double alpha, double betas, const_diag_iter di, const_diag_iter dsi,
|
||
KronVector& d, double& eig_min) const;
|
||
void solviipComplexAndEliminate(double alpha, double betas, const_diag_iter di,
|
||
const_diag_iter dsi, KronVector& d, double& eig_min) const;
|
||
// Eliminations
|
||
void solviEliminateReal(const_diag_iter di, KronVector& d, const KronVector& y,
|
||
double divisor) const;
|
||
void solviEliminateComplex(const_diag_iter di, KronVector& d, const KronVector& y1,
|
||
const KronVector& y2, double divisor) const;
|
||
void solviipEliminateReal(const_diag_iter di, const_diag_iter dsi, KronVector& d,
|
||
const KronVector& y1, const KronVector& y2, double divisor,
|
||
double divisor2) const;
|
||
void solviipEliminateComplex(const_diag_iter di, const_diag_iter dsi, KronVector& d,
|
||
const KronVector& y1, const KronVector& y11, const KronVector& y2,
|
||
const KronVector& y22, double divisor) const;
|
||
// Lemma 2
|
||
void solviipComplex(double alpha, double betas, double gamma, double delta1, double delta2,
|
||
KronVector& d1, KronVector& d2, double& eig_min) const;
|
||
// Norms for what we consider zero on diagonal of F
|
||
static constexpr double diag_zero = 1.e-15;
|
||
static constexpr double diag_zero_sq = diag_zero * diag_zero;
|
||
};
|
||
|
||
#endif
|