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 "t_container.hh"
22 #include "kron_prod.hh"
23 #include "ps_tensor.hh"
24 #include "pyramid_prod.hh"
25
26 // UGSContainer conversion from FGSContainer
UGSContainer(const FGSContainer & c)27 UGSContainer::UGSContainer(const FGSContainer &c)
28 : TensorContainer<UGSTensor>(c.num())
29 {
30 for (const auto &it : c)
31 insert(std::make_unique<UGSTensor>(*(it.second)));
32 }
33
34 /* We set ‘l’ to dimension of ‘t’, this is a tensor which multiplies tensors
35 from the container from the left. Also we set ‘k’ to a dimension of the
36 resulting tensor. We go through all equivalences on ‘k’ element set and
37 pickup only those which have ‘l’ classes.
38
39 In each loop, we fetch all necessary tensors for the product to the vector
40 ‘ts’. Then we form Kronecker product KronProdAll and feed it with tensors
41 from ‘ts’. Then we form unfolded permuted symmetry tensor UPSTensor as
42 matrix product of ‘t’ and Kronecker product ‘kp’. Then we add the permuted
43 data to ‘out’. This is done by UPSTensor method addTo(). */
44
45 void
multAndAdd(const UGSTensor & t,UGSTensor & out) const46 UGSContainer::multAndAdd(const UGSTensor &t, UGSTensor &out) const
47 {
48 int l = t.dimen();
49 int k = out.dimen();
50 const EquivalenceSet &eset = TLStatic::getEquiv(k);
51
52 for (const auto &it : eset)
53 if (it.numClasses() == l)
54 {
55 std::vector<const UGSTensor *> ts = fetchTensors(out.getSym(), it);
56 KronProdAllOptim kp(l);
57 for (int i = 0; i < l; i++)
58 kp.setMat(i, *(ts[i]));
59 kp.optimizeOrder();
60 UPSTensor ups(out.getDims(), it, t, kp);
61 ups.addTo(out);
62 }
63 }
64
65 // FGSContainer conversion from UGSContainer
FGSContainer(const UGSContainer & c)66 FGSContainer::FGSContainer(const UGSContainer &c)
67 : TensorContainer<FGSTensor>(c.num())
68 {
69 for (const auto &it : c)
70 insert(std::make_unique<FGSTensor>(*(it.second)));
71 }
72
73 // FGSContainer::multAndAdd() folded code
74 /* Here we perform one step of the Faà Di Bruno operation. We call the
75 multAndAdd() for unfolded tensor. */
76 void
multAndAdd(const FGSTensor & t,FGSTensor & out) const77 FGSContainer::multAndAdd(const FGSTensor &t, FGSTensor &out) const
78 {
79 UGSTensor ut(t);
80 multAndAdd(ut, out);
81 }
82
83 // FGSContainer::multAndAdd() unfolded code
84 /* This is the same as UGSContainer::multAndAdd() but we do not construct
85 UPSTensor from the Kronecker product, but FPSTensor. */
86 void
multAndAdd(const UGSTensor & t,FGSTensor & out) const87 FGSContainer::multAndAdd(const UGSTensor &t, FGSTensor &out) const
88 {
89 int l = t.dimen();
90 int k = out.dimen();
91 const EquivalenceSet &eset = TLStatic::getEquiv(k);
92
93 for (const auto &it : eset)
94 if (it.numClasses() == l)
95 {
96 std::vector<const FGSTensor *> ts
97 = fetchTensors(out.getSym(), it);
98 KronProdAllOptim kp(l);
99 for (int i = 0; i < l; i++)
100 kp.setMat(i, *(ts[i]));
101 kp.optimizeOrder();
102 FPSTensor fps(out.getDims(), it, t, kp);
103 fps.addTo(out);
104 }
105 }
106
107 /* This fills a given vector with integer sequences corresponding to first
108 ‘num’ indices from interval ‘start’ (including) to ‘end’ (excluding). If
109 there are not ‘num’ of such indices, the shorter vector is returned. */
110 Tensor::index
getIndices(int num,std::vector<IntSequence> & out,const Tensor::index & start,const Tensor::index & end)111 FGSContainer::getIndices(int num, std::vector<IntSequence> &out,
112 const Tensor::index &start,
113 const Tensor::index &end)
114 {
115 out.clear();
116 int i = 0;
117 Tensor::index run = start;
118 while (i < num && run != end)
119 {
120 out.push_back(run.getCoor());
121 i++;
122 ++run;
123 }
124 return run;
125 }
126