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 stacked tensor columns. 22 23 /* We need to calculate the following tensor product: 24 ⱼ ₗ 25 [f_sʲ]_α₁…αⱼ = ∑ [f_zˡ]_β₁…βₗ ∑ ∏ [z_{s^|cₘ|}]_cₘ(α)^βₘ 26 ˡ⁼¹ c∈ℳₗ,ₖ ᵐ⁼¹ 27 where s=(y,u,u′,σ), and z is a composition of four variables, say (v,w,y,u). 28 Note that z ends with y and u, and the only non-zero derivative of the 29 trailing part of z involving y or u is the first derivative and is the unit 30 matrix y_y=(1) or uᵤ=(1). Also, we suppose that the dependence of v, and w 31 on s is such that whenever derivative of w is nonzero, then also of v. This 32 means that there for any derivative and any index there is a continuous part 33 of derivatives of v and optionally of w followed by column of zeros 34 containing at most one 1. 35 36 This structure can be modelled and exploited with some costs at 37 programming. For example, let us consider the following product: 38 39 [B_y²u³]_α₁α₂β₁β₂β₃ = … + [f_z³]_γ₁γ₂γ₃ [z_yu]^γ₁_α₁β₁ [z_y]^γ₂_α₂ [zᵤᵤ]^γ₃_β₂β₃ + … 40 41 The term corresponds to equivalence { {0,2}, {1}, {3,4} }. For the fixed 42 index α₁α₂β₁β₂β₃ we have to make a Kronecker product of the columns 43 44 [z_yu]_α₁β₁ ⊗ [z_y]_α₂ ⊗ [zᵤᵤ]_β₂β₃ 45 46 which can be written as: 47 48 ⎛[v_yu]_α₁β₁⎞ ⎛[v_y]_α₂⎞ ⎛[vᵤᵤ]_β₂β₃⎞ 49 ⎢[w_yu]_α₁β₁⎥ ⎢[w_y]_α₂⎥ ⎢[wᵤᵤ]_β₂β₃⎥ 50 ⎢ 0 ⎥⊗⎢ 1_α₂ ⎥⊗⎢ 0 ⎥ 51 ⎝ 0 ⎠ ⎝ 0 ⎠ ⎝ 0 ⎠ 52 where 1_α₂ is a column of zeros having the only 1 at α₂ index. 53 54 This file develops the abstraction for this Kronecker product column 55 without multiplication of the zeros at the top. Basically, it will be 56 a column which is a Kronecker product of the columns without the 57 zeros: 58 59 ⎛[v_y]_α₂⎞ 60 ⎛[v_yu]_α₁β₁⎞⊗ ⎢[w_y]_α₂⎥⊗⎛[vᵤᵤ]_β₂β₃⎞ 61 ⎝[w_yu]_α₁β₁⎠ ⎝ 1 ⎠ ⎝[wᵤᵤ]_β₂β₃⎠ 62 63 The class will have a tensor infrastructure introducing ‘index’ which 64 iterates over all items in the column with γ₁γ₂γ₃ 65 as coordinates in [f_z³]. The data of such a tensor is 66 not suitable for any matrix operation and will have to be accessed 67 only through the ‘index’. Note that this does not matter, since 68 [f_zˡ] are sparse. */ 69 70 #ifndef PYRAMID_PROD2_H 71 #define PYRAMID_PROD2_H 72 73 #include "permutation.hh" 74 #include "tensor.hh" 75 #include "tl_exception.hh" 76 #include "rfs_tensor.hh" 77 #include "stack_container.hh" 78 79 #include "Vector.hh" 80 81 #include <vector> 82 83 /* First we declare a helper class for the tensor. Its purpose is to 84 gather the columns which are going to be Kronecker multiplied. The 85 input of this helper class is StackProduct<FGSTensor> and coordinate 86 ‘c’ of the column. 87 88 It maintains ‘unit_flag’ array which says for what columns we must 89 stack 1 below v and w. In this case, the value of ‘unit_flag’ is 90 an index of the 1, otherwise the value of ‘unit_flag’ is -1. 91 92 Also we have storage for the stacked columns ‘cols’. The object is 93 responsible for memory management associated to this storage. That is 94 why we do not allow any copy constructor, since we need to be sure 95 that no accidental copies take place. We declare the copy constructor 96 as private and not implement it. */ 97 98 class IrregTensor; 99 class IrregTensorHeader 100 { 101 friend class IrregTensor; 102 int nv; 103 IntSequence unit_flag; 104 std::vector<std::unique_ptr<Vector>> cols; 105 IntSequence end_seq; 106 public: 107 IrregTensorHeader(const StackProduct<FGSTensor> &sp, const IntSequence &c); 108 int dimen() const109 dimen() const 110 { 111 return unit_flag.size(); 112 } 113 void increment(IntSequence &v) const; 114 int calcMaxOffset() const; 115 }; 116 117 /* Here we declare the irregular tensor. There is no special logic here. We 118 inherit from Tensor and we must implement three methods, increment(), 119 decrement() and getOffset(). The last two are not implemented now, since 120 they are not needed, and they raise an exception. The first just calls 121 increment() of the header. Also we declare a method addTo() which adds this 122 unfolded irregular single column tensor to folded (regular) single column 123 tensor. 124 125 The header IrregTensorHeader lives with an object by a reference. This is 126 dangerous. However, we will use this class only in a simple loop and both 127 IrregTensor and IrregTensorHeader will be destructed at the end of a block. 128 Since the super class Tensor must be initialized before any member, we could 129 do either a save copy of IrregTensorHeader, or relatively dangerous the 130 reference member. For the reason above we chose the latter. */ 131 132 class IrregTensor : public Tensor 133 { 134 const IrregTensorHeader &header; 135 public: 136 IrregTensor(const IrregTensorHeader &h); 137 void addTo(FRSingleTensor &out) const; 138 void increment(IntSequence & v) const139 increment(IntSequence &v) const override 140 { 141 header.increment(v); 142 } 143 void decrement(IntSequence & v) const144 decrement(IntSequence &v) const override 145 { 146 TL_RAISE("Not implemented error in IrregTensor::decrement"); 147 } 148 int getOffset(const IntSequence & v) const149 getOffset(const IntSequence &v) const override 150 { 151 TL_RAISE("Not implemented error in IrregTensor::getOffset"); 152 } 153 }; 154 155 #endif 156