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 "pyramid_prod2.hh"
22 #include "rfs_tensor.hh"
23
24 /* Here we only call sp.createPackedColumns(c, cols, unit_flag) which
25 fills ‘cols’ and ‘unit_flag’ for the given column ‘c’. Then we set
26 ‘end_seq’ according to ‘unit_flag’ and columns lengths. */
27
IrregTensorHeader(const StackProduct<FGSTensor> & sp,const IntSequence & c)28 IrregTensorHeader::IrregTensorHeader(const StackProduct<FGSTensor> &sp,
29 const IntSequence &c)
30 : nv(sp.getAllSize()),
31 unit_flag(sp.dimen()),
32 cols(sp.createPackedColumns(c, unit_flag)),
33 end_seq(sp.dimen())
34 {
35 for (int i = 0; i < sp.dimen(); i++)
36 {
37 end_seq[i] = cols[i]->length();
38 if (unit_flag[i] != -1)
39 end_seq[i] = unit_flag[i]+1;
40 }
41 }
42
43 /* Here we have to increment the given integer sequence. We do it by
44 the following code, whose pattern is valid for all tensor. The only
45 difference is how we increment item of coordinates. */
46
47 void
increment(IntSequence & v) const48 IrregTensorHeader::increment(IntSequence &v) const
49 {
50 TL_RAISE_IF(v.size() != dimen(),
51 "Wrong size of coordinates in IrregTensorHeader::increment");
52
53 if (v.size() == 0)
54 return;
55 int i = v.size()-1;
56
57 // increment i-th item in coordinate ‘v’
58 /* Here we increment item of coordinates. Whenever we reached end of
59 column coming from matrices, and ‘unit_flag’ is not -1, we have to
60 jump to that ‘unit_flag’. */
61 v[i]++;
62 if (unit_flag[i] != -1 && v[i] == cols[i]->length()-1)
63 v[i] = unit_flag[i];
64
65 while (i > 0 && v[i] == end_seq[i])
66 {
67 v[i] = 0;
68 i--;
69 // increment i-th item in coordinate ‘v’
70 /* Same code as above */
71 v[i]++;
72 if (unit_flag[i] != -1 && v[i] == cols[i]->length()-1)
73 v[i] = unit_flag[i];
74 }
75 }
76
77 /* It is a product of all column lengths. */
78
79 int
calcMaxOffset() const80 IrregTensorHeader::calcMaxOffset() const
81 {
82 int res = 1;
83 for (int i = 0; i < dimen(); i++)
84 res *= cols[i]->length();
85 return res;
86 }
87
88 /* Everything is done in IrregTensorHeader, only we have to Kronecker
89 multiply all columns of the header. */
90
IrregTensor(const IrregTensorHeader & h)91 IrregTensor::IrregTensor(const IrregTensorHeader &h)
92 : Tensor(indor::along_row, IntSequence(h.dimen(), 0), h.end_seq,
93 h.calcMaxOffset(), 1, h.dimen()),
94 header(h)
95 {
96 if (header.dimen() == 1)
97 {
98 getData() = *(header.cols[0]);
99 return;
100 }
101
102 auto last = std::make_unique<Vector>(*(header.cols[header.dimen()-1]));
103 for (int i = header.dimen()-2; i > 0; i--)
104 {
105 auto newlast = std::make_unique<Vector>(last->length()*header.cols[i]->length());
106 KronProd::kronMult(ConstVector(*(header.cols[i])),
107 ConstVector(*last), *newlast);
108 last = std::move(newlast);
109 }
110 KronProd::kronMult(ConstVector(*(header.cols[0])),
111 ConstVector(*last), getData());
112 }
113
114 void
addTo(FRSingleTensor & out) const115 IrregTensor::addTo(FRSingleTensor &out) const
116 {
117 for (index it = begin(); it != end(); ++it)
118 {
119 IntSequence tmp(it.getCoor());
120 tmp.sort();
121 Tensor::index ind(out, tmp);
122 out.get(*ind, 0) += get(*it, 0);
123 }
124 }
125