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