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 #include "kron_prod.hh"
22 #include "tl_exception.hh"
23
24 #include <tuple>
25
26 /* Here we construct Kronecker product dimensions from Kronecker product
27 dimensions by picking a given matrix and all other set to identity. The
28 constructor takes dimensions of A₁⊗A₂⊗…⊗Aₙ and makes dimensions of A₁⊗I,
29 I⊗Aᵢ⊗I or I⊗Aₙ for a given i. The identity matrices must fit into the
30 described order. See header file.
31
32 We first decide what is the length of the resulting dimensions. Possible
33 length is three for I⊗A⊗I, and two for I⊗A or A⊗I.
34
35 Then we fork according to i. */
36
KronProdDimens(const KronProdDimens & kd,int i)37 KronProdDimens::KronProdDimens(const KronProdDimens &kd, int i)
38 : rows((i == 0 || i == kd.dimen()-1) ? 2 : 3),
39 cols((i == 0 || i == kd.dimen()-1) ? 2 : 3)
40 {
41 TL_RAISE_IF(i < 0 || i >= kd.dimen(),
42 "Wrong index for pickup in KronProdDimens constructor");
43
44 int kdim = kd.dimen();
45 if (i == 0)
46 {
47 // set A⊗I dimensions
48 /* The first rows and cols are taken from ‘kd’. The dimension of the
49 identity matrix is the number of rows in A₂⊗…⊗Aₙ since the matrix A₁⊗I
50 is the first. */
51 rows[0] = kd.rows[0];
52 rows[1] = kd.rows.mult(1, kdim);
53 cols[0] = kd.cols[0];
54 cols[1] = rows[1];
55 }
56 else if (i == kdim-1)
57 {
58 // set I⊗A dimensions
59 /* The second dimension is taken from ‘kd’. The dimension of the identity
60 matrix is the number of columns of A₁⊗…⊗Aₙ₋₁, since the matrix I⊗Aₙ is
61 the last. */
62 rows[0] = kd.cols.mult(0, kdim-1);
63 rows[1] = kd.rows[kdim-1];
64 cols[0] = rows[0];
65 cols[1] = kd.cols[kdim-1];
66 }
67 else
68 {
69 // set I⊗A⊗I dimensions
70 /* The dimensions of the middle matrix are taken from ‘kd’. The dimension
71 of the first identity matrix is the number of columns of A₁⊗…⊗Aᵢ₋₁,
72 and the dimension of the last identity matrix is the number of rows of
73 Aᵢ₊₁⊗…⊗Aₙ. */
74 rows[0] = kd.cols.mult(0, i);
75 cols[0] = rows[0];
76 rows[1] = kd.rows[i];
77 cols[1] = kd.cols[i];
78 cols[2] = kd.rows.mult(i+1, kdim);
79 rows[2] = cols[2];
80 }
81 }
82
83 /* This raises an exception if dimensions are bad for multiplication
84 out = in·this. */
85
86 void
checkDimForMult(const ConstTwoDMatrix & in,const TwoDMatrix & out) const87 KronProd::checkDimForMult(const ConstTwoDMatrix &in, const TwoDMatrix &out) const
88 {
89 auto [my_rows, my_cols] = kpd.getRC();
90 TL_RAISE_IF(in.nrows() != out.nrows() || in.ncols() != my_rows,
91 "Wrong dimensions for KronProd in KronProd::checkDimForMult");
92 }
93
94 /* Here we Kronecker multiply two given vectors v₁ and v₂ and
95 store the result in preallocated ‘res’. */
96
97 void
kronMult(const ConstVector & v1,const ConstVector & v2,Vector & res)98 KronProd::kronMult(const ConstVector &v1, const ConstVector &v2,
99 Vector &res)
100 {
101 TL_RAISE_IF(res.length() != v1.length()*v2.length(),
102 "Wrong vector lengths in KronProd::kronMult");
103 res.zeros();
104 for (int i = 0; i < v1.length(); i++)
105 {
106 Vector sub(res, i *v2.length(), v2.length());
107 sub.add(v1[i], v2);
108 }
109 }
110
111 void
setMat(int i,const TwoDMatrix & m)112 KronProdAll::setMat(int i, const TwoDMatrix &m)
113 {
114 matlist[i] = &m;
115 kpd.setRC(i, m.nrows(), m.ncols());
116 }
117
118 void
setUnit(int i,int n)119 KronProdAll::setUnit(int i, int n)
120 {
121 matlist[i] = nullptr;
122 kpd.setRC(i, n, n);
123 }
124
125 bool
isUnit() const126 KronProdAll::isUnit() const
127 {
128 int i = 0;
129 while (i < dimen() && matlist[i] == nullptr)
130 i++;
131 return i == dimen();
132 }
133
134 /* Here we compute B·(I⊗A). If m is the dimension of the identity matrix, and B
135 is partitioned accordingly, then the result is (B₁·A B₂·A … Bₘ·A).
136
137 In the implementation, ‘outi’ are partitions of ‘out’, ‘ini’ are const
138 partitions of ‘in’, and ‘id_cols’ is m. We employ level-2 BLAS. */
139
140 void
mult(const ConstTwoDMatrix & in,TwoDMatrix & out) const141 KronProdIA::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
142 {
143 checkDimForMult(in, out);
144
145 int id_cols = kpd.cols[0];
146 ConstTwoDMatrix a(mat);
147
148 for (int i = 0; i < id_cols; i++)
149 {
150 TwoDMatrix outi(out, i * a.ncols(), a.ncols());
151 ConstTwoDMatrix ini(in, i * a.nrows(), a.nrows());
152 outi.mult(ini, a);
153 }
154 }
155
156 /* Here we construct KronProdAI from KronProdIAI. It is clear. */
KronProdAI(const KronProdIAI & kpiai)157 KronProdAI::KronProdAI(const KronProdIAI &kpiai)
158 : KronProd(KronProdDimens(2)), mat(kpiai.mat)
159 {
160 kpd.rows[0] = mat.nrows();
161 kpd.cols[0] = mat.ncols();
162 kpd.rows[1] = kpiai.kpd.rows[2];
163 kpd.cols[1] = kpiai.kpd.cols[2];
164 }
165
166 /* Here we compute B·(A⊗I). Let the dimension of the
167 matrix A be m×n, the dimension of I be p, and the number
168 of rows of B be q. We use the fact that:
169
170 B·(A⊗I)=reshape(reshape(B, q, mp)·A, q, np).
171
172 This works only for matrix B, whose storage has leading dimension equal to
173 number of rows.
174
175 For cases where the leading dimension is not equal to the number of
176 rows, we partition the matrix A⊗I into m×n square partitions aᵢⱼI.
177 Therefore, we partition B into m partitions (B₁ B₂ … Bₘ). Each partition of B has the same number of
178 columns as the identity matrix. If R denotes the resulting matrix,
179 then it can be partitioned into n partitions
180 (R₁ R₂ … Rₙ). Each partition of R has the same number of
181 columns as the identity matrix. Then we have Rᵢ=∑aⱼᵢBⱼ.
182
183 In the implementation, ‘outi’ is Rᵢ, ‘ini’ is Bⱼ, and ‘id_cols’ is the
184 dimension of the identity matrix. */
185
186 void
mult(const ConstTwoDMatrix & in,TwoDMatrix & out) const187 KronProdAI::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
188 {
189 checkDimForMult(in, out);
190
191 int id_cols = kpd.cols[1];
192 ConstTwoDMatrix a(mat);
193
194 if (in.getLD() == in.nrows())
195 {
196 ConstTwoDMatrix in_resh(in.nrows()*id_cols, a.nrows(), in.getData());
197 TwoDMatrix out_resh(in.nrows()*id_cols, a.ncols(), out.getData());
198 out_resh.mult(in_resh, a);
199 }
200 else
201 {
202 out.zeros();
203 for (int i = 0; i < a.ncols(); i++)
204 {
205 TwoDMatrix outi(out, i *id_cols, id_cols);
206 for (int j = 0; j < a.nrows(); j++)
207 {
208 ConstTwoDMatrix ini(in, j *id_cols, id_cols);
209 outi.add(a.get(j, i), ini);
210 }
211 }
212 }
213 }
214
215 /* Here we compute R=B·(I⊗A⊗I). If n is the dimension of the first identity
216 matrix, and B is partitioned accordingly, then the result is:
217 (B₁·(A⊗I) B₂·(A⊗I) … Bₘ·(A⊗I)).
218 Each Bᵢ·(A⊗I) is in fact KronProdAI::mult(). Note that the number of columns
219 of partitions of B is the number of rows of A⊗I, and the number of columns
220 of partitions of R is the number of columns of A⊗I.
221
222 In the implementation, ‘id_cols’ is n, ‘akronid’ is A⊗I, and ‘in_bl_width’
223 and ‘out_bl_width’ are the rows and cols of A⊗I. */
224
225 void
mult(const ConstTwoDMatrix & in,TwoDMatrix & out) const226 KronProdIAI::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
227 {
228 checkDimForMult(in, out);
229
230 int id_cols = kpd.cols[0];
231
232 KronProdAI akronid(*this);
233 auto [in_bl_width, out_bl_width] = akronid.kpd.getRC();
234
235 for (int i = 0; i < id_cols; i++)
236 {
237 TwoDMatrix outi(out, i *out_bl_width, out_bl_width);
238 ConstTwoDMatrix ini(in, i *in_bl_width, in_bl_width);
239 akronid.mult(ini, outi);
240 }
241 }
242
243 /* Here we compute B·(A₁⊗…⊗Aₙ). First we compute B·(A₁⊗I), then this is
244 multiplied by all I⊗Aᵢ⊗I, and finally by I⊗Aₙ.
245
246 If the dimension of the Kronecker product is only 1, then we multiply
247 two matrices in a straight manner and return.
248
249 The intermediate results are stored on heap pointed by ‘last’. A new
250 result is allocated, and then the former storage is deallocated.
251
252 We have to be careful in cases when last or first matrix is unit and
253 no calculations are performed in corresponding codes. The codes should
254 handle ‘last’ safely also if no calcs are done. */
255
256 void
mult(const ConstTwoDMatrix & in,TwoDMatrix & out) const257 KronProdAll::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
258 {
259 // quick copy if product is unit
260 if (isUnit())
261 {
262 out.zeros();
263 out.add(1.0, in);
264 return;
265 }
266
267 // quick zero if one of the matrices is zero
268 /* If one of the matrices is exactly zero or the ‘in’ matrix is zero,
269 set out to zero and return */
270 bool is_zero = false;
271 for (int i = 0; i < dimen() && !is_zero; i++)
272 is_zero = matlist[i] && matlist[i]->isZero();
273 if (is_zero || in.isZero())
274 {
275 out.zeros();
276 return;
277 }
278
279 // quick multiplication if dimension is 1
280 if (dimen() == 1)
281 {
282 if (matlist[0]) // always true
283 out.mult(in, ConstTwoDMatrix(*(matlist[0])));
284 return;
285 }
286
287 int c;
288 std::unique_ptr<TwoDMatrix> last;
289
290 // perform first multiplication by A₁⊗I
291 /* Here we have to construct A₁⊗I, allocate intermediate result ‘last’, and
292 perform the multiplication. */
293 if (matlist[0])
294 {
295 KronProdAI akronid(*this);
296 c = akronid.kpd.ncols();
297 last = std::make_unique<TwoDMatrix>(in.nrows(), c);
298 akronid.mult(in, *last);
299 }
300 else
301 last = std::make_unique<TwoDMatrix>(in.nrows(), in.ncols(), Vector{in.getData()});
302
303 // perform intermediate multiplications by I⊗Aᵢ⊗I
304 /* Here we go through all I⊗Aᵢ⊗I, construct the product, allocate new storage
305 for result ‘newlast’, perform the multiplication and set ‘last’ to
306 ‘newlast’. */
307 for (int i = 1; i < dimen()-1; i++)
308 if (matlist[i])
309 {
310 KronProdIAI interkron(*this, i);
311 c = interkron.kpd.ncols();
312 auto newlast = std::make_unique<TwoDMatrix>(in.nrows(), c);
313 interkron.mult(*last, *newlast);
314 last = std::move(newlast);
315 }
316
317 // perform last multiplication by I⊗Aₙ
318 if (matlist[dimen()-1])
319 {
320 KronProdIA idkrona(*this);
321 idkrona.mult(*last, out);
322 }
323 else
324 out = *last;
325 }
326
327 /* This calculates a Kornecker product of rows of matrices, the row
328 indices are given by the integer sequence. The result is allocated and
329 returned. */
330
331 std::unique_ptr<Vector>
multRows(const IntSequence & irows) const332 KronProdAll::multRows(const IntSequence &irows) const
333 {
334 TL_RAISE_IF(irows.size() != dimen(),
335 "Wrong length of row indices in KronProdAll::multRows");
336
337 std::unique_ptr<Vector> last;
338 std::unique_ptr<ConstVector> row;
339 std::vector<std::unique_ptr<Vector>> to_delete;
340 for (int i = 0; i < dimen(); i++)
341 {
342 int j = dimen()-1-i;
343
344 // set ‘row’ to the number of rows of j-th matrix
345 /* If the j-th matrix is a real matrix, then the row is constructed from
346 the matrix. It the matrix is unit, we construct a new vector, fill it
347 with zeros, than set the unit to appropriate place, and make the ‘row’
348 as ConstVector of this vector, which is sheduled for deallocation. */
349 if (matlist[j])
350 row = std::make_unique<ConstVector>(matlist[j]->getRow(irows[j]));
351 else
352 {
353 auto aux = std::make_unique<Vector>(ncols(j));
354 aux->zeros();
355 (*aux)[irows[j]] = 1.0;
356 row = std::make_unique<ConstVector>(*aux);
357 to_delete.emplace_back(std::move(aux));
358 }
359
360 // set ‘last’ to product of ‘row’ and ‘last’
361 /* If the ‘last’ exists, we allocate new storage and Kronecker multiply.
362 If the ‘last’ does not exist, then we only make ‘last’ equal to
363 ‘row’. */
364 if (last)
365 {
366 auto newlast = std::make_unique<Vector>(last->length()*row->length());
367 kronMult(*row, ConstVector(*last), *newlast);
368 last = std::move(newlast);
369 }
370 else
371 last = std::make_unique<Vector>(*row);
372 }
373
374 return last;
375 }
376
377 /* This permutes the matrices so that the new ordering would minimize memory
378 consumption. As shown in KronProdAllOptim class declaration, we want:
379
380 mₖ/nₖ ≤ mₖ₋₁/nₖ₋₁ ≤ … ≤ m₁/n₁
381
382 where mᵢ×nᵢ is the dimension of Aᵢ. So we implement the bubble sort. */
383
384 void
optimizeOrder()385 KronProdAllOptim::optimizeOrder()
386 {
387 for (int i = 0; i < dimen(); i++)
388 {
389 int swaps = 0;
390 for (int j = 0; j < dimen()-1; j++)
391 {
392 if (static_cast<double>(kpd.rows[j])/kpd.cols[j]
393 < static_cast<double>(kpd.rows[j+1])/kpd.cols[j+1])
394 {
395 // swap dimensions and matrices at j and j+1
396 int s = kpd.rows[j+1];
397 kpd.rows[j+1] = kpd.rows[j];
398 kpd.rows[j] = s;
399 s = kpd.cols[j+1];
400 kpd.cols[j+1] = kpd.cols[j];
401 kpd.cols[j] = s;
402 const TwoDMatrix *m = matlist[j+1];
403 matlist[j+1] = matlist[j];
404 matlist[j] = m;
405
406 // project the swap to the permutation ‘oper’
407 s = oper.getMap()[j+1];
408 oper.getMap()[j+1] = oper.getMap()[j];
409 oper.getMap()[j] = s;
410 swaps++;
411 }
412 }
413 if (swaps == 0)
414 return;
415 }
416 }
417