1 /* 2 * Copyright © 2004 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 // Kronecker product. 22 23 /* Here we define an abstraction for a Kronecker product of a sequence of 24 matrices, i.e. A₁⊗…⊗Aₙ. Obviously we do not store the product in memory. 25 First we need to represent a dimension of the Kronecker product. Then we 26 represent the Kronecker product, simply it is the Kronecker product 27 dimension with a vector of references to the matrices A₁,…,Aₙ. 28 29 The main task of this class is to calculate a matrix product B·(A₁⊗A₂⊗…⊗Aₙ), 30 which in our application has much more moderate dimensions than A₁⊗A₂⊗…⊗Aₙ. 31 We calculate it as B·(A₁⊗I)·…·(I⊗Aᵢ⊗I)·…·(I⊗Aₙ) where dimensions of identity 32 matrices differ and are given by the chosen order. One can naturally ask, 33 whether there is some optimal order minimizing maximum storage needed for 34 intermediate results. The optimal ordering is implemented by class 35 KronProdAllOptim. 36 37 For this multiplication, we also need to represent products of type A⊗I, 38 I⊗A⊗I and I⊗A. */ 39 40 #ifndef KRON_PROD_H 41 #define KRON_PROD_H 42 43 #include <utility> 44 #include <vector> 45 #include <memory> 46 47 #include "twod_matrix.hh" 48 #include "permutation.hh" 49 #include "int_sequence.hh" 50 51 class KronProdAll; 52 class KronProdAllOptim; 53 class KronProdIA; 54 class KronProdIAI; 55 class KronProdAI; 56 57 /* KronProdDimens maintains a dimension of the Kronecker product. So, it 58 maintains two sequences, one for rows, and one for columns. */ 59 60 class KronProdDimens 61 { 62 friend class KronProdAll; 63 friend class KronProdAllOptim; 64 friend class KronProdIA; 65 friend class KronProdIAI; 66 friend class KronProdAI; 67 private: 68 IntSequence rows; 69 IntSequence cols; 70 public: 71 // Initializes to a given dimension, and all rows and cols are set to zeros KronProdDimens(int dim)72 KronProdDimens(int dim) 73 : rows(dim, 0), cols(dim, 0) 74 { 75 } 76 KronProdDimens(const KronProdDimens &kd) = default; 77 KronProdDimens(KronProdDimens &&kd) = default; 78 79 /* Takes dimensions of A₁⊗…⊗Aₙ, and makes dimensions of A₁⊗I, I⊗Aᵢ⊗I or I⊗Aₙ 80 for a given i. The dimensions of identity matrices are such that 81 A₁⊗…⊗Aₙ=(A₁⊗I)·…·(I⊗Aᵢ⊗I)·…·(I⊗Aₙ). Note that the matrices on the right do 82 not commute only because sizes of identity matrices which are then given 83 by this ordering. */ 84 KronProdDimens(const KronProdDimens &kd, int i); 85 86 KronProdDimens &operator=(const KronProdDimens &kd) = default; 87 KronProdDimens &operator=(KronProdDimens &&kd) = default; 88 bool operator ==(const KronProdDimens & kd) const89 operator==(const KronProdDimens &kd) const 90 { 91 return rows == kd.rows && cols == kd.cols; 92 } 93 94 int dimen() const95 dimen() const 96 { 97 return rows.size(); 98 } 99 void setRC(int i,int r,int c)100 setRC(int i, int r, int c) 101 { 102 rows[i] = r; 103 cols[i] = c; 104 } 105 std::pair<int, int> getRC(int i) const106 getRC(int i) const 107 { 108 return { rows[i], cols[i] }; 109 } 110 std::pair<int, int> getRC() const111 getRC() const 112 { 113 return { rows.mult(), cols.mult() }; 114 } 115 int nrows() const116 nrows() const 117 { 118 return rows.mult(); 119 } 120 int ncols() const121 ncols() const 122 { 123 return cols.mult(); 124 } 125 int nrows(int i) const126 nrows(int i) const 127 { 128 return rows[i]; 129 } 130 int ncols(int i) const131 ncols(int i) const 132 { 133 return cols[i]; 134 } 135 }; 136 137 /* Here we define an abstract class for all Kronecker product classes, which 138 are KronProdAll (the most general), KronProdIA (for I⊗A), KronProdAI (for 139 A⊗I), and KronProdIAI (for I⊗A⊗I). The purpose of the super class is to only 140 define some common methods and common member ‘kpd’ for dimensions and 141 declare pure virtual mult() which is implemented by the subclasses. 142 143 The class also contains a static method kronMult(), which calculates a 144 Kronecker product of two vectors and stores it in the provided vector. It is 145 useful at a few points of the library. */ 146 147 class KronProd 148 { 149 protected: 150 KronProdDimens kpd; 151 public: KronProd(int dim)152 KronProd(int dim) 153 : kpd(dim) 154 { 155 } KronProd(const KronProdDimens & kd)156 KronProd(const KronProdDimens &kd) 157 : kpd(kd) 158 { 159 } 160 KronProd(const KronProd &kp) = default; 161 KronProd(KronProd &&kp) = default; 162 virtual ~KronProd() = default; 163 164 int dimen() const165 dimen() const 166 { 167 return kpd.dimen(); 168 } 169 170 virtual void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const = 0; 171 void mult(const TwoDMatrix & in,TwoDMatrix & out) const172 mult(const TwoDMatrix &in, TwoDMatrix &out) const 173 { 174 mult(ConstTwoDMatrix(in), out); 175 } 176 177 void checkDimForMult(const ConstTwoDMatrix &in, const TwoDMatrix &out) const; 178 void checkDimForMult(const TwoDMatrix & in,const TwoDMatrix & out) const179 checkDimForMult(const TwoDMatrix &in, const TwoDMatrix &out) const 180 { 181 checkDimForMult(ConstTwoDMatrix(in), out); 182 } 183 184 static void kronMult(const ConstVector &v1, const ConstVector &v2, 185 Vector &res); 186 187 int nrows() const188 nrows() const 189 { 190 return kpd.nrows(); 191 } 192 int ncols() const193 ncols() const 194 { 195 return kpd.ncols(); 196 } 197 int nrows(int i) const198 nrows(int i) const 199 { 200 return kpd.nrows(i); 201 } 202 int ncols(int i) const203 ncols(int i) const 204 { 205 return kpd.ncols(i); 206 } 207 }; 208 209 /* KronProdAll is the main class of this file. It represents the Kronecker 210 product A₁⊗A₂⊗…⊗Aₙ. Besides dimensions, it stores pointers to matrices in 211 ‘matlist’ array. If a pointer is null, then the matrix is considered to be 212 unit. The array is set by calls to setMat() method (for real matrices) or 213 setUnit() method (for unit matrices). 214 215 The object is constructed by a constructor, which allocates the ‘matlist’ 216 and initializes dimensions to zeros. Then a caller must feed the object with 217 matrices by calling setMat() and setUnit() repeatedly for different indices. 218 219 We implement the mult() method of KronProd, and a new method multRows(), 220 which creates a vector of kronecker product of all rows of matrices in the 221 object. The rows are given by the IntSequence. */ 222 223 class KronProdAll : public KronProd 224 { 225 friend class KronProdIA; 226 friend class KronProdIAI; 227 friend class KronProdAI; 228 protected: 229 std::vector<const TwoDMatrix *> matlist; 230 public: KronProdAll(int dim)231 KronProdAll(int dim) 232 : KronProd(dim), matlist(dim) 233 { 234 } 235 ~KronProdAll() override = default; 236 void setMat(int i, const TwoDMatrix &m); 237 void setUnit(int i, int n); 238 const TwoDMatrix & getMat(int i) const239 getMat(int i) const 240 { 241 return *(matlist[i]); 242 } 243 244 void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const override; 245 std::unique_ptr<Vector> multRows(const IntSequence &irows) const; 246 private: 247 bool isUnit() const; 248 }; 249 250 /* The class KronProdAllOptim minimizes memory consumption of the product 251 B·(A₁⊗A₂⊗…⊗Aₖ). The optimization is done by reordering of the matrices 252 A₁,…,Aₖ, in order to minimize a sum of all storages needed for intermediate 253 results. The optimal ordering is also nearly optimal with respect to number 254 of flops. 255 256 Let (mᵢ,nᵢ) be dimensions of Aᵢ. It is easy to observe, that for the i-th 257 step we need storage of r·(n₁·…·nᵢ·mᵢ₊₁·…·mₖ), where r is the number of rows 258 of B. To minimize the sum through all i over all permutations of matrices, 259 it is equivalent to minimizing the sum 260 ₖ 261 ∑ (mᵢ₊₁·…·mₖ)/(nᵢ₊₁·…·nₖ) 262 ⁱ⁼¹ 263 The optimal ordering will yield mₖ/nₖ ≤ mₖ₋₁/nₖ₋₁ ≤ … ≤ m₁/n₁ 264 265 Now observe, that the number of flops for the i-th step is 266 r·(n₁·…·nᵢ·mᵢ·…·mₖ). In order to 267 minimize the number of flops, it is equivalent to minimize 268 ₖ 269 ∑ mᵢ(mᵢ₊₁·…·mₖ)/(nᵢ₊₁·…·nₖ) 270 ⁱ⁼¹ 271 Note that, normally, mᵢ does not change as much as nᵢ₊₁,…,nₖ, so the 272 ordering minimizing the memory will be nearly optimal with respect to number 273 of flops. 274 275 The class KronProdAllOptim inherits from KronProdAll. A public method 276 optimizeOrder() does the reordering. The permutation is stored in ‘oper’. So 277 as long as optimizeOrder() is not called, the class is equivalent to 278 KronProdAll. */ 279 280 class KronProdAllOptim : public KronProdAll 281 { 282 protected: 283 Permutation oper; 284 public: KronProdAllOptim(int dim)285 KronProdAllOptim(int dim) 286 : KronProdAll(dim), oper(dim) 287 { 288 } 289 void optimizeOrder(); 290 const Permutation & getPer() const291 getPer() const 292 { 293 return oper; 294 } 295 }; 296 297 /* This class represents I⊗A. We have only one reference to the matrix, which 298 is set by constructor. */ 299 300 class KronProdIA : public KronProd 301 { 302 friend class KronProdAll; 303 const TwoDMatrix &mat; 304 public: KronProdIA(const KronProdAll & kpa)305 KronProdIA(const KronProdAll &kpa) 306 : KronProd(KronProdDimens(kpa.kpd, kpa.dimen()-1)), 307 mat(kpa.getMat(kpa.dimen()-1)) 308 { 309 } 310 void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const override; 311 }; 312 313 /* This class represents A⊗I. We have only one reference to the matrix, which 314 is set by constructor. */ 315 316 class KronProdAI : public KronProd 317 { 318 friend class KronProdIAI; 319 friend class KronProdAll; 320 const TwoDMatrix &mat; 321 public: KronProdAI(const KronProdAll & kpa)322 KronProdAI(const KronProdAll &kpa) 323 : KronProd(KronProdDimens(kpa.kpd, 0)), 324 mat(kpa.getMat(0)) 325 { 326 } 327 KronProdAI(const KronProdIAI &kpiai); 328 329 void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const override; 330 }; 331 332 /* This class represents I⊗A⊗I. We have only one reference to the matrix, which 333 is set by constructor. */ 334 335 class KronProdIAI : public KronProd 336 { 337 friend class KronProdAI; 338 friend class KronProdAll; 339 const TwoDMatrix &mat; 340 public: KronProdIAI(const KronProdAll & kpa,int i)341 KronProdIAI(const KronProdAll &kpa, int i) 342 : KronProd(KronProdDimens(kpa.kpd, i)), 343 mat(kpa.getMat(i)) 344 { 345 } 346 void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const override; 347 }; 348 349 #endif 350