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 // Sparse tensor.
22 
23 /* Here we declare a sparse full and general symmetry tensors with the
24    multidimensional index along columns. We implement them as a std::multimap
25    associating to each sequence of coordinates IntSequence a set of pairs (row,
26    number). This is very convenient but not optimal in terms of memory
27    consumption. So the implementation can be changed.
28 
29    The current std::multimap implementation allows insertions. Another
30    advantage of this approach is that we do not need to calculate column
31    numbers from the IntSequence, since the column is accessed directly via the
32    key which is IntSequence.
33 
34    The only operation we need to do with the full symmetry sparse tensor
35    is a left multiplication of a row oriented single column tensor. The
36    result of such an operation is a column of the same size as the sparse
37    tensor. Other important operations are slicing operations. We need to
38    do sparse and dense slices of full symmetry sparse tensors. In fact,
39    the only constructor of general symmetry sparse tensor is slicing from
40    the full symmetry sparse. */
41 
42 #ifndef SPARSE_TENSOR_H
43 #define SPARSE_TENSOR_H
44 
45 #include "symmetry.hh"
46 #include "tensor.hh"
47 #include "gs_tensor.hh"
48 #include "Vector.hh"
49 
50 #include <map>
51 
52 struct ltseq
53 {
54   bool
operator ()ltseq55   operator()(const IntSequence &s1, const IntSequence &s2) const
56   {
57     return s1 < s2;
58   }
59 };
60 
61 /* This is a super class of both full symmetry and general symmetry sparse
62    tensors. It contains a std::multimap and implements insertions. It tracks
63    maximum and minimum row, for which there is an item. */
64 
65 class SparseTensor
66 {
67 public:
68   using Map = std::multimap<IntSequence, std::pair<int, double>, ltseq>;
69 protected:
70   Map m;
71   int dim;
72   int nr;
73   int nc;
74   int first_nz_row;
75   int last_nz_row;
76 public:
SparseTensor(int d,int nnr,int nnc)77   SparseTensor(int d, int nnr, int nnc)
78     : dim(d), nr(nnr), nc(nnc), first_nz_row(nr), last_nz_row(-1)
79   {
80   }
81   void insert(IntSequence s, int r, double c);
82   const Map &
getMap() const83   getMap() const
84   {
85     return m;
86   }
87   int
dimen() const88   dimen() const
89   {
90     return dim;
91   }
92   int
nrows() const93   nrows() const
94   {
95     return nr;
96   }
97   int
ncols() const98   ncols() const
99   {
100     return nc;
101   }
102   double
getFillFactor() const103   getFillFactor() const
104   {
105     return static_cast<double>(m.size())/nrows()/ncols();
106   }
107   double getFoldIndexFillFactor() const;
108   double getUnfoldIndexFillFactor() const;
109   int
getNumNonZero() const110   getNumNonZero() const
111   {
112     return m.size();
113   }
114   int
getFirstNonZeroRow() const115   getFirstNonZeroRow() const
116   {
117     return first_nz_row;
118   }
119   int
getLastNonZeroRow() const120   getLastNonZeroRow() const
121   {
122     return last_nz_row;
123   }
124   virtual const Symmetry &getSym() const = 0;
125   void print() const;
126   bool isFinite() const;
127 };
128 
129 /* This is a full symmetry sparse tensor. It implements multColumnAndAdd() and,
130    in addition to SparseTensor, it has ‘nv’ (number of variables) and symmetry
131    (basically it is a dimension). */
132 
133 class FSSparseTensor : public SparseTensor
134 {
135 private:
136   int nv;
137   Symmetry sym;
138 public:
139   FSSparseTensor(int d, int nvar, int r);
140   void insert(IntSequence s, int r, double c);
141   void multColumnAndAdd(const Tensor &t, Vector &v) const;
142   const Symmetry &
getSym() const143   getSym() const override
144   {
145     return sym;
146   }
147   int
nvar() const148   nvar() const
149   {
150     return nv;
151   }
152   void print() const;
153 };
154 
155 /* This is a general symmetry sparse tensor. It has TensorDimens and can be
156    constructed as a slice of the full symmetry sparse tensor. The slicing
157    constructor takes the same form as the slicing FGSTensor constructor from
158    full symmetry sparse tensor. */
159 
160 class GSSparseTensor : public SparseTensor
161 {
162 private:
163   TensorDimens tdims;
164 public:
165   GSSparseTensor(const FSSparseTensor &t, const IntSequence &ss,
166                  const IntSequence &coor, TensorDimens td);
167   void insert(IntSequence s, int r, double c);
168   const Symmetry &
getSym() const169   getSym() const override
170   {
171     return tdims.getSym();
172   }
173   const TensorDimens &
getDims() const174   getDims() const
175   {
176     return tdims;
177   }
178   void print() const;
179 
180 };
181 
182 #endif
183