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_prod.hh"
22 #include "permutation.hh"
23 #include "tl_exception.hh"
24
25 /* Here we construct the USubTensor object. We allocate space via the parent
26 URTensor. Number of columns is a length of the list of indices ‘lst’, number
27 of variables and dimensions are of the tensor $h$, this is given by ‘hdims’.
28
29 We go through all equivalences with number of classes equal to dimension of
30 B. For each equivalence we make a permutation ‘per’. Then we fetch all the
31 necessary tensors g with symmetries implied by symmetry of B and the
32 equivalence. Then we go through the list of indices, permute them by the
33 permutation and add the Kronecker product of the selected columns. This is
34 done by addKronColumn(). */
35
USubTensor(const TensorDimens & bdims,const TensorDimens & hdims,const FGSContainer & cont,const std::vector<IntSequence> & lst)36 USubTensor::USubTensor(const TensorDimens &bdims,
37 const TensorDimens &hdims,
38 const FGSContainer &cont,
39 const std::vector<IntSequence> &lst)
40 : URTensor(lst.size(), hdims.getNVX(0), hdims.dimen())
41 {
42 TL_RAISE_IF(!hdims.getNVX().isConstant(),
43 "Tensor has not full symmetry in USubTensor()");
44 const EquivalenceSet &eset = TLStatic::getEquiv(bdims.dimen());
45 zeros();
46 for (const auto &it : eset)
47 {
48 if (it.numClasses() == hdims.dimen())
49 {
50 Permutation per(it);
51 std::vector<const FGSTensor *> ts
52 = cont.fetchTensors(bdims.getSym(), it);
53 for (int i = 0; i < static_cast<int>(lst.size()); i++)
54 {
55 IntSequence perindex(lst[i].size());
56 per.apply(lst[i], perindex);
57 addKronColumn(i, ts, perindex);
58 }
59 }
60 }
61 }
62
63 /* This makes a Kronecker product of appropriate columns from tensors in ‘fs’
64 and adds such data to i-th column of this matrix. The appropriate columns
65 are defined by ‘pindex’ sequence. A column of a tensor has index created
66 from a corresponding part of ‘pindex’. The sizes of these parts are given by
67 dimensions of the tensors in ‘ts’.
68
69 Here we break the given index ‘pindex’ according to the dimensions of the
70 tensors in ‘ts’, and for each subsequence of the ‘pindex’ we find an index
71 of the folded tensor, which involves calling getOffset() for folded tensor,
72 which might be costly. We gather all columns to a vector ‘tmpcols’ which are
73 Kronecker multiplied in constructor of URSingleTensor. Finally we add data
74 of URSingleTensor to the i-th column. */
75
76 void
addKronColumn(int i,const std::vector<const FGSTensor * > & ts,const IntSequence & pindex)77 USubTensor::addKronColumn(int i, const std::vector<const FGSTensor *> &ts,
78 const IntSequence &pindex)
79 {
80 std::vector<ConstVector> tmpcols;
81 int lastdim = 0;
82 for (auto t : ts)
83 {
84 IntSequence ind(pindex, lastdim, lastdim+t->dimen());
85 lastdim += t->dimen();
86 index in(*t, ind);
87 tmpcols.push_back(t->getCol(*in));
88 }
89
90 URSingleTensor kronmult(tmpcols);
91 Vector coli{getCol(i)};
92 coli.add(1.0, kronmult.getData());
93 }
94