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 // General symmetry tensor. 22 23 /* Here we define tensors for general symmetry. All tensors from here are 24 identifying the multidimensional index with columns. Thus all 25 symmetries regard to columns. The general symmetry here is not the most 26 general. It captures all symmetries of indices which are given by 27 continuous partitioning of indices. Two items are symmetric if they 28 belong to the same group. The continuity implies that if two items 29 belong to one group, then all items between them belong to that 30 group. This continuous partitioning of indices is described by 31 Symmetry class. 32 33 The dimension of the tensors here are described (besides the symmetry) 34 also by number of variables for each group. This is dealt in the class 35 for tensor dimensions defined also here. */ 36 37 #ifndef GS_TENSOR_H 38 #define GS_TENSOR_H 39 40 #include "tensor.hh" 41 #include "fs_tensor.hh" 42 #include "symmetry.hh" 43 #include "rfs_tensor.hh" 44 45 class FGSTensor; 46 class UGSTensor; 47 class FSSparseTensor; 48 49 /* This class encapsulates symmetry information for the general 50 symmetry tensor. It maintains a vector of variable numbers ‘nvs’, and 51 symmetry ‘sym’. For example, let the symmetry be y²u³, and 52 variable numbers be 10 for y, and 5 for u. Then the ‘nvs’ is 53 (10,5), and ‘sym’ is (2,3). Also it maintains ‘nvmax’ unfolded ‘nvs’ with 54 respect to the symmetry, this is (10,10,5,5,5). 55 56 The class is able to calculate number of offsets (columns or rows depending 57 what matrix coordinate we describe) in unfolded and folded tensors 58 with the given symmetry. */ 59 60 class TensorDimens 61 { 62 protected: 63 IntSequence nvs; 64 Symmetry sym; 65 IntSequence nvmax; 66 public: TensorDimens(Symmetry s,IntSequence nvars)67 TensorDimens(Symmetry s, IntSequence nvars) 68 : nvs(std::move(nvars)), sym(std::move(s)), nvmax(nvs.unfold(sym)) 69 { 70 } 71 // Full-symmetry special case TensorDimens(int nvar,int dimen)72 TensorDimens(int nvar, int dimen) 73 : nvs{nvar}, sym{dimen}, nvmax(dimen, nvar) 74 { 75 } 76 // Constructs the tensor dimensions for slicing (see the implementation for details) 77 TensorDimens(const IntSequence &ss, const IntSequence &coor); 78 79 bool operator ==(const TensorDimens & td) const80 operator==(const TensorDimens &td) const 81 { 82 return nvs == td.nvs && sym == td.sym; 83 } 84 bool operator !=(const TensorDimens & td) const85 operator!=(const TensorDimens &td) const 86 { 87 return !operator==(td); 88 } 89 90 int dimen() const91 dimen() const 92 { 93 return sym.dimen(); 94 } 95 int getNVX(int i) const96 getNVX(int i) const 97 { 98 return nvmax[i]; 99 } 100 const IntSequence & getNVS() const101 getNVS() const 102 { 103 return nvs; 104 } 105 const IntSequence & getNVX() const106 getNVX() const 107 { 108 return nvmax; 109 } 110 const Symmetry & getSym() const111 getSym() const 112 { 113 return sym; 114 } 115 116 int calcUnfoldMaxOffset() const; 117 int calcFoldMaxOffset() const; 118 int calcFoldOffset(const IntSequence &v) const; 119 void decrement(IntSequence &v) const; 120 }; 121 122 /* Here is a class for folded general symmetry tensor. It only contains 123 tensor dimensions, it defines types for indices, implement virtual 124 methods of super class FTensor. */ 125 126 class GSSparseTensor; 127 class FGSTensor : public FTensor 128 { 129 friend class UGSTensor; 130 131 const TensorDimens tdims; 132 public: FGSTensor(int r,TensorDimens td)133 FGSTensor(int r, TensorDimens td) 134 : FTensor(indor::along_col, td.getNVX(), r, 135 td.calcFoldMaxOffset(), td.dimen()), tdims(std::move(td)) 136 { 137 } 138 FGSTensor(const FGSTensor &) = default; 139 FGSTensor(FGSTensor &&) = default; 140 FGSTensor(int first_row,int num,FGSTensor & t)141 FGSTensor(int first_row, int num, FGSTensor &t) 142 : FTensor(first_row, num, t), tdims(t.tdims) 143 { 144 } 145 146 // Constructs a slice from a fully symmetric sparse tensor 147 FGSTensor(const FSSparseTensor &t, const IntSequence &ss, 148 const IntSequence &coor, TensorDimens td); 149 150 // Constructs a slice from a fully symmetric dense tensor 151 FGSTensor(const FFSTensor &t, const IntSequence &ss, 152 const IntSequence &coor, TensorDimens td); 153 154 // Converting constructors 155 explicit FGSTensor(const UGSTensor &ut); 156 explicit FGSTensor(const GSSparseTensor &sp); FGSTensor(FFSTensor & t)157 explicit FGSTensor(FFSTensor &t) 158 : FTensor(0, t.nrows(), t), tdims(t.nvar(), t.dimen()) 159 { 160 } 161 162 ~FGSTensor() override = default; 163 164 void increment(IntSequence &v) const override; 165 void decrement(IntSequence & v) const166 decrement(IntSequence &v) const override 167 { 168 tdims.decrement(v); 169 } 170 std::unique_ptr<UTensor> unfold() const override; 171 const TensorDimens & getDims() const172 getDims() const 173 { 174 return tdims; 175 } 176 const Symmetry & getSym() const177 getSym() const 178 { 179 return getDims().getSym(); 180 } 181 182 /* Performs a contraction of one variable in the tensor. This is, for 183 instance: 184 185 [r_xⁱzᵏ]_α₁…αᵢγ₁…γₖ = [t_xⁱyʲzᵏ]_α₁…αᵢβ₁…βⱼγ₁…γₖ·[c]^β₁…βⱼ 186 */ 187 void contractAndAdd(int i, FGSTensor &out, 188 const FRSingleTensor &col) const; 189 190 int getOffset(const IntSequence & v) const191 getOffset(const IntSequence &v) const override 192 { 193 return tdims.calcFoldOffset(v); 194 } 195 }; 196 197 /* Besides similar things that has FGSTensor, we have here also 198 method unfoldData(), and helper method getFirstIndexOf() 199 which corresponds to sorting coordinates in fully symmetric case (here 200 the action is more complicated, so we put it to the method). */ 201 202 class UGSTensor : public UTensor 203 { 204 friend class FGSTensor; 205 206 const TensorDimens tdims; 207 public: UGSTensor(int r,TensorDimens td)208 UGSTensor(int r, TensorDimens td) 209 : UTensor(indor::along_col, td.getNVX(), r, 210 td.calcUnfoldMaxOffset(), td.dimen()), tdims(std::move(td)) 211 { 212 } 213 214 UGSTensor(const UGSTensor &) = default; 215 UGSTensor(UGSTensor &&) = default; 216 UGSTensor(int first_row,int num,UGSTensor & t)217 UGSTensor(int first_row, int num, UGSTensor &t) 218 : UTensor(first_row, num, t), tdims(t.tdims) 219 { 220 } 221 222 // Constructs a slice from fully symmetric sparse tensor 223 UGSTensor(const FSSparseTensor &t, const IntSequence &ss, 224 const IntSequence &coor, TensorDimens td); 225 226 // Constructs a slice from fully symmetric dense unfolded tensor 227 UGSTensor(const UFSTensor &t, const IntSequence &ss, 228 const IntSequence &coor, TensorDimens td); 229 230 // Converting constructors 231 explicit UGSTensor(const FGSTensor &ft); UGSTensor(UFSTensor & t)232 explicit UGSTensor(UFSTensor &t) 233 : UTensor(0, t.nrows(), t), tdims(t.nvar(), t.dimen()) 234 { 235 } 236 237 ~UGSTensor() override = default; 238 239 void increment(IntSequence &v) const override; 240 void decrement(IntSequence &v) const override; 241 std::unique_ptr<FTensor> fold() const override; 242 const TensorDimens & getDims() const243 getDims() const 244 { 245 return tdims; 246 } 247 const Symmetry & getSym() const248 getSym() const 249 { 250 return getDims().getSym(); 251 } 252 253 void contractAndAdd(int i, UGSTensor &out, 254 const URSingleTensor &col) const; 255 int getOffset(const IntSequence &v) const override; 256 private: 257 void unfoldData(); 258 public: 259 index getFirstIndexOf(const index &in) const; 260 }; 261 262 #endif 263