1 /* 2 * Copyright © 2004-2011 Ondra Kamenik 3 * Copyright © 2019 Dynare Team 4 * 5 * This file is part of Dynare. 6 * 7 * Dynare is free software: you can redistribute it and/or modify 8 * it under the terms of the GNU General Public License as published by 9 * the Free Software Foundation, either version 3 of the License, or 10 * (at your option) any later version. 11 * 12 * Dynare is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU General Public License 18 * along with Dynare. If not, see <http://www.gnu.org/licenses/>. 19 */ 20 21 #ifndef QUASI_TRIANGULAR_H 22 #define QUASI_TRIANGULAR_H 23 24 #include "Vector.hh" 25 #include "KronVector.hh" 26 #include "SylvMatrix.hh" 27 28 #include <list> 29 #include <memory> 30 31 class DiagonalBlock; 32 class Diagonal; 33 class DiagPair 34 { 35 private: 36 double *a1; 37 double *a2; 38 public: 39 DiagPair() = default; DiagPair(double * aa1,double * aa2)40 DiagPair(double *aa1, double *aa2) : a1{aa1}, a2{aa2} 41 { 42 } 43 DiagPair(const DiagPair &p) = default; 44 DiagPair &operator=(const DiagPair &p) = default; 45 DiagPair & operator =(double v)46 operator=(double v) 47 { 48 *a1 = v; 49 *a2 = v; 50 return *this; 51 } 52 const double & operator *() const53 operator*() const 54 { 55 return *a1; 56 } 57 /* Here we must not define double& operator*(), since it wouldn't 58 rewrite both values, we use operator=() for this */ 59 friend class Diagonal; 60 friend class DiagonalBlock; 61 }; 62 63 /* Stores a diagonal block of a quasi-triangular real matrix: 64 – either a 1×1 block, i.e. a real scalar, stored in α₁ 65 ⎛α₁ β₁⎞ 66 – or a 2×2 block, stored as ⎝β₂ α₂⎠ 67 */ 68 class DiagonalBlock 69 { 70 private: 71 int jbar; // Index of block in the diagonal 72 bool real; 73 DiagPair alpha; 74 double *beta1; 75 double *beta2; 76 77 public: 78 DiagonalBlock() = default; DiagonalBlock(int jb,bool r,double * a1,double * a2,double * b1,double * b2)79 DiagonalBlock(int jb, bool r, double *a1, double *a2, 80 double *b1, double *b2) 81 : jbar{jb}, real{r}, alpha{a1, a2}, beta1{b1}, beta2{b2} 82 { 83 } 84 // Construct a complex 2×2 block 85 /* β₁ and β₂ will be deduced from pointers to α₁ and α₂ */ DiagonalBlock(int jb,double * a1,double * a2)86 DiagonalBlock(int jb, double *a1, double *a2) 87 : jbar{jb}, real{false}, alpha{a1, a2}, beta1{a2-1}, beta2{a1+1} 88 { 89 } 90 // Construct a real 1×1 block DiagonalBlock(int jb,double * a1)91 DiagonalBlock(int jb, double *a1) 92 : jbar{jb}, real{true}, alpha{a1, a1}, beta1{nullptr}, beta2{nullptr} 93 { 94 } 95 DiagonalBlock(const DiagonalBlock &b) = default; 96 DiagonalBlock &operator=(const DiagonalBlock &b) = default; 97 int getIndex() const98 getIndex() const 99 { 100 return jbar; 101 } 102 bool isReal() const103 isReal() const 104 { 105 return real; 106 } 107 const DiagPair & getAlpha() const108 getAlpha() const 109 { 110 return alpha; 111 } 112 DiagPair & getAlpha()113 getAlpha() 114 { 115 return alpha; 116 } 117 double & getBeta1() const118 getBeta1() const 119 { 120 return *beta1; 121 } 122 double & getBeta2() const123 getBeta2() const 124 { 125 return *beta2; 126 } 127 // Returns determinant of this block (assuming it is 2×2) 128 double getDeterminant() const; 129 // Returns −β₁β₂ 130 double getSBeta() const; 131 // Returns the modulus of the eigenvalue(s) contained in this block 132 double getSize() const; 133 // Transforms this block into a real one 134 void setReal(); 135 // Verifies that the block information is consistent with the matrix d (for debugging) 136 void checkBlock(const double *d, int d_size); 137 friend class Diagonal; 138 }; 139 140 // Stores the diagonal blocks of a quasi-triangular real matrix 141 class Diagonal 142 { 143 public: 144 using const_diag_iter = std::list<DiagonalBlock>::const_iterator; 145 using diag_iter = std::list<DiagonalBlock>::iterator; 146 private: 147 int num_all{0}; // Total number of blocks 148 std::list<DiagonalBlock> blocks; 149 int num_real{0}; // Number of 1×1 (real) blocks 150 public: 151 Diagonal() = default; 152 // Construct the diagonal blocks of (quasi-triangular) matrix ‘data’ 153 Diagonal(double *data, int d_size); 154 /* Construct the diagonal blocks of (quasi-triangular) matrix ‘data’, 155 assuming it has the same shape as ‘d’ */ 156 Diagonal(double *data, const Diagonal &d); 157 Diagonal(const Diagonal &d) = default; 158 Diagonal &operator=(const Diagonal &d) = default; 159 virtual ~Diagonal() = default; 160 161 // Returns number of 2×2 blocks on the diagonal 162 int getNumComplex() const163 getNumComplex() const 164 { 165 return num_all - num_real; 166 } 167 // Returns number of 1×1 blocks on the diagonal 168 int getNumReal() const169 getNumReal() const 170 { 171 return num_real; 172 } 173 // Returns number of scalar elements on the diagonal 174 int getSize() const175 getSize() const 176 { 177 return getNumReal() + 2*getNumComplex(); 178 } 179 // Returns total number of blocks on the diagonal 180 int getNumBlocks() const181 getNumBlocks() const 182 { 183 return num_all; 184 } 185 void getEigenValues(Vector &eig) const; 186 void swapLogically(diag_iter it); 187 void checkConsistency(diag_iter it); 188 double getAverageSize(diag_iter start, diag_iter end); 189 diag_iter findClosestBlock(diag_iter start, diag_iter end, double a); 190 diag_iter findNextLargerBlock(diag_iter start, diag_iter end, double a); 191 void print() const; 192 193 diag_iter begin()194 begin() 195 { 196 return blocks.begin(); 197 } 198 const_diag_iter begin() const199 begin() const 200 { 201 return blocks.begin(); 202 } 203 diag_iter end()204 end() 205 { 206 return blocks.end(); 207 } 208 const_diag_iter end() const209 end() const 210 { 211 return blocks.end(); 212 } 213 214 /* redefine pointers as data start at p */ 215 void changeBase(double *p); 216 private: 217 constexpr static double EPS = 1.0e-300; 218 /* Computes number of 2×2 diagonal blocks on the quasi-triangular matrix 219 represented by data (of size d_size×d_size) */ 220 static int getNumComplex(const double *data, int d_size); 221 // Checks whether |p|<EPS 222 static bool isZero(double p); 223 }; 224 225 template<class _TRef, class _TPtr> 226 struct _matrix_iter 227 { 228 using _Self = _matrix_iter<_TRef, _TPtr>; 229 int d_size; 230 bool real; 231 _TPtr ptr; 232 public: _matrix_iter_matrix_iter233 _matrix_iter(_TPtr base, int ds, bool r) 234 { 235 ptr = base; 236 d_size = ds; 237 real = r; 238 } 239 virtual ~_matrix_iter() = default; 240 bool operator ==_matrix_iter241 operator==(const _Self &it) const 242 { 243 return ptr == it.ptr; 244 } 245 bool operator !=_matrix_iter246 operator!=(const _Self &it) const 247 { 248 return ptr != it.ptr; 249 } 250 _TRef operator *_matrix_iter251 operator*() const 252 { 253 return *ptr; 254 } 255 _TRef a_matrix_iter256 a() const 257 { 258 return *ptr; 259 } 260 virtual _Self &operator++() = 0; 261 }; 262 263 template<class _TRef, class _TPtr> 264 class _column_iter : public _matrix_iter<_TRef, _TPtr> 265 { 266 using _Tparent = _matrix_iter<_TRef, _TPtr>; 267 using _Self = _column_iter<_TRef, _TPtr>; 268 int row; 269 public: _column_iter(_TPtr base,int ds,bool r,int rw)270 _column_iter(_TPtr base, int ds, bool r, int rw) 271 : _matrix_iter<_TRef, _TPtr>(base, ds, r), row(rw) 272 { 273 }; 274 _Self & operator ++()275 operator++() override 276 { 277 _Tparent::ptr++; 278 row++; 279 return *this; 280 } 281 _TRef b() const282 b() const 283 { 284 if (_Tparent::real) 285 return *(_Tparent::ptr); 286 else 287 return *(_Tparent::ptr+_Tparent::d_size); 288 } 289 int getRow() const290 getRow() const 291 { 292 return row; 293 } 294 }; 295 296 template<class _TRef, class _TPtr> 297 class _row_iter : public _matrix_iter<_TRef, _TPtr> 298 { 299 using _Tparent = _matrix_iter<_TRef, _TPtr>; 300 using _Self = _row_iter<_TRef, _TPtr>; 301 int col; 302 public: _row_iter(_TPtr base,int ds,bool r,int cl)303 _row_iter(_TPtr base, int ds, bool r, int cl) 304 : _matrix_iter<_TRef, _TPtr>(base, ds, r), col(cl) 305 { 306 }; 307 _Self & operator ++()308 operator++() override 309 { 310 _Tparent::ptr += _Tparent::d_size; 311 col++; 312 return *this; 313 } 314 virtual _TRef b() const315 b() const 316 { 317 if (_Tparent::real) 318 return *(_Tparent::ptr); 319 else 320 return *(_Tparent::ptr+1); 321 } 322 int getCol() const323 getCol() const 324 { 325 return col; 326 } 327 }; 328 329 class SchurDecomp; 330 class SchurDecompZero; 331 332 /* Represents an upper quasi-triangular matrix. 333 All the elements are stored in the SqSylvMatrix super-class. 334 Additionally, a list of the diagonal blocks (1×1 or 2×2), is stored in the 335 “diagonal” member, in order to optimize some operations (where the matrix is 336 seen as an upper-triangular matrix, plus sub-diagonal elements of the 2×2 337 diagonal blocks) */ 338 class QuasiTriangular : public SqSylvMatrix 339 { 340 public: 341 using const_col_iter = _column_iter<const double &, const double *>; 342 using col_iter = _column_iter<double &, double *>; 343 using const_row_iter = _row_iter<const double &, const double *>; 344 using row_iter = _row_iter<double &, double *>; 345 using const_diag_iter = Diagonal::const_diag_iter; 346 using diag_iter = Diagonal::diag_iter; 347 protected: 348 Diagonal diagonal; 349 public: 350 QuasiTriangular(const ConstVector &d, int d_size); 351 // Initializes with r·t 352 QuasiTriangular(double r, const QuasiTriangular &t); 353 // Initializes with r·t+r₂·t₂ 354 QuasiTriangular(double r, const QuasiTriangular &t, 355 double r2, const QuasiTriangular &t2); 356 // Initializes with t² 357 QuasiTriangular(const std::string &dummy, const QuasiTriangular &t); 358 explicit QuasiTriangular(const SchurDecomp &decomp); 359 explicit QuasiTriangular(const SchurDecompZero &decomp); 360 QuasiTriangular(const QuasiTriangular &t); 361 362 ~QuasiTriangular() override = default; 363 const Diagonal & getDiagonal() const364 getDiagonal() const 365 { 366 return diagonal; 367 } 368 int getNumOffdiagonal() const; 369 void swapDiagLogically(diag_iter it); 370 void checkDiagConsistency(diag_iter it); 371 double getAverageDiagSize(diag_iter start, diag_iter end); 372 diag_iter findClosestDiagBlock(diag_iter start, diag_iter end, double a); 373 diag_iter findNextLargerBlock(diag_iter start, diag_iter end, double a); 374 375 /* (I+this)·y = x, y→x */ 376 virtual void solvePre(Vector &x, double &eig_min); 377 /* (I+thisᵀ)·y = x, y→x */ 378 virtual void solvePreTrans(Vector &x, double &eig_min); 379 /* (I+this)·x = b */ 380 virtual void solve(Vector &x, const ConstVector &b, double &eig_min); 381 /* (I+thisᵀ)·x = b */ 382 virtual void solveTrans(Vector &x, const ConstVector &b, double &eig_min); 383 /* x = this·b */ 384 virtual void multVec(Vector &x, const ConstVector &b) const; 385 /* x = thisᵀ·b */ 386 virtual void multVecTrans(Vector &x, const ConstVector &b) const; 387 /* x = x + this·b */ 388 virtual void multaVec(Vector &x, const ConstVector &b) const; 389 /* x = x + thisᵀ·b */ 390 virtual void multaVecTrans(Vector &x, const ConstVector &b) const; 391 /* x = (this⊗I)·x */ 392 virtual void multKron(KronVector &x) const; 393 /* x = (thisᵀ⊗I)·x */ 394 virtual void multKronTrans(KronVector &x) const; 395 /* A = this·A */ 396 virtual void multLeftOther(GeneralMatrix &a) const; 397 /* A = thisᵀ·A */ 398 virtual void multLeftOtherTrans(GeneralMatrix &a) const; 399 400 const_diag_iter diag_begin() const401 diag_begin() const 402 { 403 return diagonal.begin(); 404 } 405 diag_iter diag_begin()406 diag_begin() 407 { 408 return diagonal.begin(); 409 } 410 const_diag_iter diag_end() const411 diag_end() const 412 { 413 return diagonal.end(); 414 } 415 diag_iter diag_end()416 diag_end() 417 { 418 return diagonal.end(); 419 } 420 421 /* iterators for off diagonal elements */ 422 virtual const_col_iter col_begin(const DiagonalBlock &b) const; 423 virtual col_iter col_begin(const DiagonalBlock &b); 424 virtual const_row_iter row_begin(const DiagonalBlock &b) const; 425 virtual row_iter row_begin(const DiagonalBlock &b); 426 virtual const_col_iter col_end(const DiagonalBlock &b) const; 427 virtual col_iter col_end(const DiagonalBlock &b); 428 virtual const_row_iter row_end(const DiagonalBlock &b) const; 429 virtual row_iter row_end(const DiagonalBlock &b); 430 431 virtual std::unique_ptr<QuasiTriangular> clone() const432 clone() const 433 { 434 return std::make_unique<QuasiTriangular>(*this); 435 } 436 // Returns this² 437 virtual std::unique_ptr<QuasiTriangular> square() const438 square() const 439 { 440 return std::make_unique<QuasiTriangular>("square", *this); 441 } 442 // Returns r·this 443 virtual std::unique_ptr<QuasiTriangular> scale(double r) const444 scale(double r) const 445 { 446 return std::make_unique<QuasiTriangular>(r, *this); 447 } 448 // Returns r·this + r₂·t₂ 449 virtual std::unique_ptr<QuasiTriangular> linearlyCombine(double r,double r2,const QuasiTriangular & t2) const450 linearlyCombine(double r, double r2, const QuasiTriangular &t2) const 451 { 452 return std::make_unique<QuasiTriangular>(r, *this, r2, t2); 453 } 454 protected: 455 // this = r·t 456 void setMatrix(double r, const QuasiTriangular &t); 457 // this = this + r·t 458 void addMatrix(double r, const QuasiTriangular &t); 459 private: 460 // this = this + I 461 void addUnit(); 462 /* x = x + (this⊗I)·b */ 463 void multaKron(KronVector &x, const ConstKronVector &b) const; 464 /* x = x + (thisᵀ⊗I)·b */ 465 void multaKronTrans(KronVector &x, const ConstKronVector &b) const; 466 /* hide noneffective implementations of parents */ 467 void multsVec(Vector &x, const ConstVector &d) const; 468 void multsVecTrans(Vector &x, const ConstVector &d) const; 469 }; 470 471 #endif /* QUASI_TRIANGULAR_H */ 472