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 // Multiplying tensor columns.
22 
23 /* In here, we implement the Faà Di Bruno for folded
24    tensors. Recall, that one step of the Faà Di Bruno is a formula:
25 
26                                  ₗ
27     [B_sᵏ]_α₁…αₗ = [h_yˡ]_γ₁…γₗ  ∏  [g_{s^|cₘ|}]_cₘ(α)^γₘ
28                                 ᵐ⁼¹
29 
30    In contrast to unfolded implementation of UGSContainer::multAndAdd()
31    with help of KronProdAll and UPSTensor, we take a completely
32    different strategy. We cannot afford full instantiation of
33 
34              ₗ
35        ∑     ∏  [g_{s^|cₘ|}]_cₘ(α)^γₘ
36     c∈ℳₗ,ₖ ᵐ⁼¹
37 
38    and therefore we do it per partes. We select some number of columns,
39    for instance 10, calculate 10 continuous iterators of tensor B. Then we
40    form unfolded tensor
41 
42                   ⎡         ₗ                       ⎤
43     [G]_S^γ₁…γₗ = ⎢   ∑     ∏  [g_{s^|cₘ|}]_cₘ(α)^γₘ⎥
44                   ⎣c∈ℳₗ,ₖ ᵐ⁼¹                       ⎦_S
45 
46    where S is the selected set of 10 indices. This is done as Kronecker
47    product of vectors corresponding to selected columns. Note that, in
48    general, there is no symmetry in G, its type is special class for
49    this purpose.
50 
51    If g is folded, then we have to form the folded version of G. There is
52    no symmetry in G data, so we sum all unfolded indices corresponding
53    to folded index together. This is perfectly OK, since we multiply
54    these groups of (equivalent) items with the same number in fully
55    symmetric g.
56 
57    After this, we perform ordinary matrix multiplication to obtain a
58    selected set of columns of B.
59 
60    In here, we define a class for forming and representing [G]_S^γ₁…γₗ.
61    Basically, this tensor is row-oriented (multidimensional index is along
62    rows), and it is fully symmetric. So we inherit from URTensor. If we need
63    its folded version, we simply use a suitable conversion. The new abstraction
64    will have only a new constructor allowing a construction from the given set
65    of indices S, and given set of tensors g. The rest of the process is
66    implemented in FGSContainer::multAndAdd() unfolded code or
67    FGSContainer::multAndAdd() folded code. */
68 
69 #ifndef PYRAMID_PROD_H
70 #define PYRAMID_PROD_H
71 
72 #include "int_sequence.hh"
73 #include "rfs_tensor.hh"
74 #include "gs_tensor.hh"
75 #include "t_container.hh"
76 
77 #include <vector>
78 
79 /* Here we define the new tensor for representing [G]_S^γ₁…γₗ. It allows a
80    construction from container of folded general symmetry tensors ‘cont’, and
81    set of indices ‘ts’. Also we have to supply dimensions of resulting tensor
82    B, and dimensions of tensor h. */
83 
84 class USubTensor : public URTensor
85 {
86 public:
87   USubTensor(const TensorDimens &bdims, const TensorDimens &hdims,
88              const FGSContainer &cont, const std::vector<IntSequence> &lst);
89   void addKronColumn(int i, const std::vector<const FGSTensor *> &ts,
90                      const IntSequence &pindex);
91 };
92 
93 #endif
94