stack_container.hh 24.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/*
 * Copyright © 2004 Ondra Kamenik
 * Copyright © 2019 Dynare Team
 *
 * This file is part of Dynare.
 *
 * Dynare is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Dynare is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
18
 * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
19
 */
20
21
22
23

// Stack of containers.

/* Here we develop abstractions for stacked containers of tensors. For
24
25
   instance, in perturbation methods for DSGE models, we need the function:

26
27
28
29
                   ⎛ G(y*,u,u′,σ) ⎞
    z(y*,u,u′,σ) = ⎢  g(y*,u,σ)   ⎥
                   ⎢      y*      ⎥
                   ⎝      u       ⎠
30
31
32

   and we need to calculate one step of Faà Di Bruno formula:

33

34
    [B_sᵏ]_α₁…αₗ = [f_zˡ]_β₁…βₗ    ∑     ∏  [z_{s^|cₘ|}]_cₘ(α)^βₘ
35
                                c∈ℳₗ,ₖ ᵐ⁼¹
36
37

   where we have containers for derivatives of G and g.
38
39

   The main purpose of this file is to define abstractions for stack of
40
41
42
   containers and possibly raw variables, and code multAndAdd() method
   calculating (one step of) the Faà Di Bruno formula for folded and
   unfolded tensors. Note also, that tensors [f_zˡ] are sparse.
43
44
45
46
47

   The abstractions are built as follows. At the top, there is an
   interface describing stack of columns. It contains pure virtual
   methods needed for manipulating the container stack. For technical
   reasons it is a template. Both versions (folded, and unfolded) provide
48
   all interface necessary for implementation of multAndAdd(). The second
49
   way of inheritance is first general implementation of the interface
50
51
52
   StackContainer, and then specific (ZContainer for our specific z).
   The only method which is virtual also after StackContainer is
   getType(), which is implemented in the specialization and determines
53
54
55
   behaviour of the stack. The complete classes are obtained by
   inheriting from the both branches, as it is drawn below:

56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

                        StackContainerInterface<FGSTensor>
                          ↙                           ↘
            StackContainer<FGSTensor>             FoldedStackContainer

              ZContainer<FGSTensor>                 ↙

                                   FoldedZContainer


                        StackContainerInterface<UGSTensor>
                          ↙                           ↘
            StackContainer<UGSTensor>           UnfoldedStackContainer

              ZContainer<UGSTensor>                 ↙

                                  UnfoldedZContainer

74

75
   We have also two supporting classes StackProduct and KronProdStack
76
77
78
79
80
81
82
83
84
85
86
87
88
89
   and a number of worker classes used as threads. */

#ifndef STACK_CONTAINER_H
#define STACK_CONTAINER_H

#include "int_sequence.hh"
#include "equivalence.hh"
#include "tl_static.hh"
#include "t_container.hh"
#include "kron_prod.hh"
#include "permutation.hh"
#include "sthread.hh"

/* Here is the general interface to stack container. The subclasses
90
91
92
   maintain IntSequence of stack sizes, i.e. size of G, g, y, and
   u. Then a convenience IntSequence of stack offsets. Then vector of
   pointers to containers, in our example G, and g.
93

94
   A non-virtual subclass must implement getType() which determines
95
96
97
98
   dependency of stack items on symmetries. There are three possible types
   for a symmetry. Either the stack item derivative wrt. the symmetry is
   a matrix, or a unit matrix, or zero.

99
100
   Method isZero() returns true if the derivative of a given stack item
   wrt. to given symmetry is zero as defined by getType() or the
101
102
103
104
   derivative is not present in the container. In this way, we can
   implement the formula conditional some of the tensors are zero, which
   is not true (they are only missing).

105
   Method createPackedColumn() returns a vector of stack derivatives with
106
   respect to the given symmetry and of the given column, where all zeros
107
108
   from zero types, or unit matrices are deleted. See kron_prod.hh for an
   explanation. */
109

110
template<class _Ttype>
111
112
113
class StackContainerInterface
{
public:
114
  using _Ctype = TensorContainer<_Ttype>;
115
  enum class itype { matrix, unit, zero };
116
public:
117
  StackContainerInterface() = default;
118
119
120
121
122
  virtual ~StackContainerInterface() = default;
  virtual const IntSequence &getStackSizes() const = 0;
  virtual IntSequence &getStackSizes() = 0;
  virtual const IntSequence &getStackOffsets() const = 0;
  virtual IntSequence &getStackOffsets() = 0;
123
  virtual int numConts() const = 0;
124
  virtual const _Ctype &getCont(int i) const = 0;
125
126
127
  virtual itype getType(int i, const Symmetry &s) const = 0;
  virtual int numStacks() const = 0;
  virtual bool isZero(int i, const Symmetry &s) const = 0;
128
  virtual const _Ttype &getMatrix(int i, const Symmetry &s) const = 0;
129
130
  virtual int getLengthOfMatrixStacks(const Symmetry &s) const = 0;
  virtual int getUnitPos(const Symmetry &s) const = 0;
131
132
133
  virtual std::unique_ptr<Vector> createPackedColumn(const Symmetry &s,
                                                     const IntSequence &coor,
                                                     int &iu) const = 0;
134
135
136
137
138
139
140
141
  int
  getAllSize() const
  {
    return getStackOffsets()[numStacks()-1]
      + getStackSizes()[numStacks()-1];
  }
};

142
143
/* Here is StackContainer, which implements almost all interface
   StackContainerInterface but one method getType() which is left for
144
145
146
   implementation to specializations.

   It does not own its tensors. */
147

148
template<class _Ttype>
149
150
151
class StackContainer : virtual public StackContainerInterface<_Ttype>
{
public:
152
153
154
  using _Stype = StackContainerInterface<_Ttype>;
  using _Ctype = typename StackContainerInterface<_Ttype>::_Ctype;
  using itype = typename StackContainerInterface<_Ttype>::itype;
155
156
157
158
protected:
  int num_conts;
  IntSequence stack_sizes;
  IntSequence stack_offsets;
159
  std::vector<const _Ctype *> conts;
160
161
public:
  StackContainer(int ns, int nc)
162
163
    : stack_sizes(ns, 0), stack_offsets(ns, 0),
      conts(nc)
164
165
166
  {
  }
  const IntSequence &
167
  getStackSizes() const override
168
169
170
171
  {
    return stack_sizes;
  }
  IntSequence &
172
  getStackSizes() override
173
174
175
176
  {
    return stack_sizes;
  }
  const IntSequence &
177
  getStackOffsets() const override
178
179
180
181
  {
    return stack_offsets;
  }
  IntSequence &
182
  getStackOffsets() override
183
184
185
186
  {
    return stack_offsets;
  }
  int
187
  numConts() const override
188
  {
189
    return conts.size();
190
  }
191
  const _Ctype &
192
  getCont(int i) const override
193
  {
194
    return *(conts[i]);
195
  }
196
  itype getType(int i, const Symmetry &s) const override = 0;
197
  int
198
  numStacks() const override
199
200
201
202
  {
    return stack_sizes.size();
  }
  bool
203
  isZero(int i, const Symmetry &s) const override
204
205
206
  {
    TL_RAISE_IF(i < 0 || i >= numStacks(),
                "Wrong index to stack in StackContainer::isZero.");
207
208
    return (getType(i, s) == itype::zero
            || (getType(i, s) == itype::matrix && !conts[i]->check(s)));
209
210
  }

211
  const _Ttype &
212
  getMatrix(int i, const Symmetry &s) const override
213
  {
214
    TL_RAISE_IF(isZero(i, s) || getType(i, s) == itype::unit,
215
216
217
218
219
                "Matrix is not returned in StackContainer::getMatrix");
    return conts[i]->get(s);
  }

  int
220
  getLengthOfMatrixStacks(const Symmetry &s) const override
221
222
223
  {
    int res = 0;
    int i = 0;
224
    while (i < numStacks() && getType(i, s) == itype::matrix)
225
226
227
228
229
      res += stack_sizes[i++];
    return res;
  }

  int
230
  getUnitPos(const Symmetry &s) const override
231
232
233
234
  {
    if (s.dimen() != 1)
      return -1;
    int i = numStacks()-1;
235
    while (i >= 0 && getType(i, s) != itype::unit)
236
237
238
239
      i--;
    return i;
  }

240
  std::unique_ptr<Vector>
241
  createPackedColumn(const Symmetry &s,
242
                     const IntSequence &coor, int &iu) const override
243
244
245
246
247
248
249
250
251
252
253
254
255
  {
    TL_RAISE_IF(s.dimen() != coor.size(),
                "Incompatible coordinates for symmetry in StackContainer::createPackedColumn");

    int len = getLengthOfMatrixStacks(s);
    iu = -1;
    int i = 0;
    if (-1 != (i = getUnitPos(s)))
      {
        iu = stack_offsets[i] + coor[0];
        len++;
      }

256
    auto res = std::make_unique<Vector>(len);
257
    i = 0;
258
    while (i < numStacks() && getType(i, s) == itype::matrix)
259
      {
260
261
        const _Ttype &t = getMatrix(i, s);
        Tensor::index ind(t, coor);
262
        Vector subres(*res, stack_offsets[i], stack_sizes[i]);
263
        subres = ConstGeneralMatrix(t).getCol(*ind);
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
        i++;
      }
    if (iu != -1)
      (*res)[len-1] = 1;

    return res;
  }

protected:
  void
  calculateOffsets()
  {
    stack_offsets[0] = 0;
    for (int i = 1; i < stack_offsets.size(); i++)
      stack_offsets[i] = stack_offsets[i-1] + stack_sizes[i-1];
  }
};

class WorkerFoldMAADense;
class WorkerFoldMAASparse1;
class WorkerFoldMAASparse2;
class WorkerFoldMAASparse4;
class FoldedStackContainer : virtual public StackContainerInterface<FGSTensor>
{
  friend class WorkerFoldMAADense;
  friend class WorkerFoldMAASparse1;
  friend class WorkerFoldMAASparse2;
  friend class WorkerFoldMAASparse4;
public:
293
  static constexpr double fill_threshold = 0.00005;
294
295
296
297
  void
  multAndAdd(int dim, const TensorContainer<FSSparseTensor> &c,
             FGSTensor &out) const
  {
298
    if (c.check(Symmetry{dim}))
299
      multAndAdd(c.get(Symmetry{dim}), out);
300
301
302
303
304
305
306
307
308
  }
  void multAndAdd(const FSSparseTensor &t, FGSTensor &out) const;
  void multAndAdd(int dim, const FGSContainer &c, FGSTensor &out) const;
protected:
  void multAndAddSparse1(const FSSparseTensor &t, FGSTensor &out) const;
  void multAndAddSparse2(const FSSparseTensor &t, FGSTensor &out) const;
  void multAndAddSparse3(const FSSparseTensor &t, FGSTensor &out) const;
  void multAndAddSparse4(const FSSparseTensor &t, FGSTensor &out) const;
  void multAndAddStacks(const IntSequence &fi, const FGSTensor &g,
309
                        FGSTensor &out, std::mutex &mut) const;
310
  void multAndAddStacks(const IntSequence &fi, const GSSparseTensor &g,
311
                        FGSTensor &out, std::mutex &mut) const;
312
313
314
315
316
317
318
319
320
321
322
};

class WorkerUnfoldMAADense;
class WorkerUnfoldMAASparse1;
class WorkerUnfoldMAASparse2;
class UnfoldedStackContainer : virtual public StackContainerInterface<UGSTensor>
{
  friend class WorkerUnfoldMAADense;
  friend class WorkerUnfoldMAASparse1;
  friend class WorkerUnfoldMAASparse2;
public:
323
  static constexpr double fill_threshold = 0.00005;
324
325
326
327
  void
  multAndAdd(int dim, const TensorContainer<FSSparseTensor> &c,
             UGSTensor &out) const
  {
328
    if (c.check(Symmetry{dim}))
329
      multAndAdd(c.get(Symmetry{dim}), out);
330
331
332
333
334
335
336
  }
  void multAndAdd(const FSSparseTensor &t, UGSTensor &out) const;
  void multAndAdd(int dim, const UGSContainer &c, UGSTensor &out) const;
protected:
  void multAndAddSparse1(const FSSparseTensor &t, UGSTensor &out) const;
  void multAndAddSparse2(const FSSparseTensor &t, UGSTensor &out) const;
  void multAndAddStacks(const IntSequence &fi, const UGSTensor &g,
337
                        UGSTensor &out, std::mutex &mut) const;
338
339
};

340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
/* Here is the specialization of the StackContainer. We implement
   here the x needed in DSGE context for welfare assessment. We implement getType() and define a constructor feeding the data and sizes.

   It depends on four variables U(y,u,u',σ), the variable u' being introduced to enable additions with 4-variable tensors*/

template<class _Ttype>
class XContainer : public StackContainer<_Ttype>
{
public:
  using _Tparent = StackContainer<_Ttype>;
  using _Stype = StackContainerInterface<_Ttype>;
  using _Ctype = typename _Tparent::_Ctype;
  using itype = typename _Tparent::itype;
  XContainer(const _Ctype *g, int ng)
    : _Tparent(1, 1)
  {
    _Tparent::stack_sizes = { ng };
    _Tparent::conts[0] = g;
    _Tparent::calculateOffsets();
  }

  /* Here we say, what happens if we derive z. recall the top of the
     file, how z looks, and code is clear. */

  itype
  getType(int i, const Symmetry &s) const override
  {
    if (i==0)
      if (s[2] > 0)
        return itype::zero;
      else
        return itype::matrix;

    TL_RAISE("Wrong stack index in XContainer::getType");
  }

};

class FoldedXContainer : public XContainer<FGSTensor>,
                         public FoldedStackContainer
{
public:
  using _Ctype = TensorContainer<FGSTensor>;
  FoldedXContainer(const _Ctype *g, int ng)
    : XContainer<FGSTensor>(g, ng)
  {
  }
};

class UnfoldedXContainer : public XContainer<UGSTensor>,
                           public UnfoldedStackContainer
{
public:
  using _Ctype = TensorContainer<UGSTensor>;
  UnfoldedXContainer(const _Ctype *g, int ng)
    : XContainer<UGSTensor>(g, ng)
  {
  }
};
399
400
/* Here is the specialization of the StackContainer. We implement
   here the z needed in DSGE context. We implement getType() and define
401
402
   a constructor feeding the data and sizes.

403
404
405
406
407
   Note that it has two containers, the first is dependent on four variables
   G(y*,u,u′,σ), and the second dependent on three variables g(y*,u,σ). So that
   we would be able to stack them, we make the second container g be dependent
   on four variables, the third being u′ a dummy and always returning zero if
   dimension of u′ is positive. */
408

409
template<class _Ttype>
410
411
412
class ZContainer : public StackContainer<_Ttype>
{
public:
413
414
415
416
  using _Tparent = StackContainer<_Ttype>;
  using _Stype = StackContainerInterface<_Ttype>;
  using _Ctype = typename _Tparent::_Ctype;
  using itype = typename _Tparent::itype;
417
418
419
420
  ZContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng,
             int ny, int nu)
    : _Tparent(4, 2)
  {
421
    _Tparent::stack_sizes = { ngss, ng, ny, nu };
422
423
424
425
426
    _Tparent::conts[0] = gss;
    _Tparent::conts[1] = g;
    _Tparent::calculateOffsets();
  }

427
428
  /* Here we say, what happens if we derive z. recall the top of the
     file, how z looks, and code is clear. */
429
430

  itype
431
  getType(int i, const Symmetry &s) const override
432
433
  {
    if (i == 0)
434
      return itype::matrix;
435
436
    if (i == 1)
      if (s[2] > 0)
437
        return itype::zero;
438
      else
439
        return itype::matrix;
440
    if (i == 2)
441
      if (s == Symmetry{1, 0, 0, 0})
442
        return itype::unit;
443
      else
444
        return itype::zero;
445
    if (i == 3)
446
      if (s == Symmetry{0, 1, 0, 0})
447
        return itype::unit;
448
      else
449
        return itype::zero;
450
451
452
453
454
455
456
457
458
459

    TL_RAISE("Wrong stack index in ZContainer::getType");
  }

};

class FoldedZContainer : public ZContainer<FGSTensor>,
                         public FoldedStackContainer
{
public:
460
  using _Ctype = TensorContainer<FGSTensor>;
461
462
463
464
465
466
467
468
469
470
471
  FoldedZContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng,
                   int ny, int nu)
    : ZContainer<FGSTensor>(gss, ngss, g, ng, ny, nu)
  {
  }
};

class UnfoldedZContainer : public ZContainer<UGSTensor>,
                           public UnfoldedStackContainer
{
public:
472
  using _Ctype = TensorContainer<UGSTensor>;
473
474
475
476
477
478
479
480
  UnfoldedZContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng,
                     int ny, int nu)
    : ZContainer<UGSTensor>(gss, ngss, g, ng, ny, nu)
  {
  }
};

/* Here we have another specialization of container used in context of
481
   DSGE. We define a container for
482

483
484
485
486
487
    G(y,u,u′,σ)=g**(g*(y,u,σ),u′,σ)

   For some reason, the symmetry of g** has length 4 although it
   is really dependent on three variables (To now the reason, consult
   the ZContainer class declaration). So, it has four stacks, the
488
   third one is dummy, and always returns zero. The first stack
489
   corresponds to a container of g*. */
490

491
template<class _Ttype>
492
493
494
class GContainer : public StackContainer<_Ttype>
{
public:
495
496
497
498
  using _Tparent = StackContainer<_Ttype>;
  using _Stype = StackContainerInterface<_Ttype>;
  using _Ctype = typename StackContainer<_Ttype>::_Ctype;
  using itype = typename StackContainer<_Ttype>::itype;
499
500
501
  GContainer(const _Ctype *gs, int ngs, int nu)
    : StackContainer<_Ttype>(4, 1)
  {
502
    _Tparent::stack_sizes = { ngs, nu, nu, 1 };
503
504
505
506
    _Tparent::conts[0] = gs;
    _Tparent::calculateOffsets();
  }

507
508
  /* Here we define the dependencies in g**(g*(y,u,σ),u′,σ). Also note, that
     first derivative of g* w.r.t. σ is always zero, so we also add this
509
510
511
     information. */

  itype
512
  getType(int i, const Symmetry &s) const override
513
514
  {
    if (i == 0)
515
      if (s[2] > 0 || s == Symmetry{0, 0, 0, 1})
516
        return itype::zero;
517
      else
518
        return itype::matrix;
519
    if (i == 1)
520
      if (s == Symmetry{0, 0, 1, 0})
521
        return itype::unit;
522
      else
523
        return itype::zero;
524
    if (i == 2)
525
      return itype::zero;
526
    if (i == 3)
527
      if (s == Symmetry{0, 0, 0, 1})
528
        return itype::unit;
529
      else
530
        return itype::zero;
531
532
533
534
535
536
537
538
539
540

    TL_RAISE("Wrong stack index in GContainer::getType");
  }

};

class FoldedGContainer : public GContainer<FGSTensor>,
                         public FoldedStackContainer
{
public:
541
  using _Ctype = TensorContainer<FGSTensor>;
542
543
544
545
546
547
548
549
550
551
  FoldedGContainer(const _Ctype *gs, int ngs, int nu)
    : GContainer<FGSTensor>(gs, ngs, nu)
  {
  }
};

class UnfoldedGContainer : public GContainer<UGSTensor>,
                           public UnfoldedStackContainer
{
public:
552
  using _Ctype = TensorContainer<UGSTensor>;
553
554
555
556
557
558
  UnfoldedGContainer(const _Ctype *gs, int ngs, int nu)
    : GContainer<UGSTensor>(gs, ngs, nu)
  {
  }
};

559
560
/* Here we have a support class for product of StackContainers. It
   only adds a dimension to StackContainer. It selects the symmetries
561
562
563
564
   according to equivalence classes passed to the constructor. The
   equivalence can have permuted classes by some given
   permutation. Nothing else is interesting. */

565
template<class _Ttype>
566
567
568
class StackProduct
{
public:
569
570
571
  using _Stype = StackContainerInterface<_Ttype>;
  using _Ctype = typename _Stype::_Ctype;
  using itype = typename _Stype::itype;
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
protected:
  const _Stype &stack_cont;
  InducedSymmetries syms;
  Permutation per;
public:
  StackProduct(const _Stype &sc, const Equivalence &e,
               const Symmetry &os)
    : stack_cont(sc), syms(e, os), per(e)
  {
  }
  StackProduct(const _Stype &sc, const Equivalence &e,
               const Permutation &p, const Symmetry &os)
    : stack_cont(sc), syms(e, p, os), per(e, p)
  {
  }
  int
  dimen() const
  {
    return syms.size();
  }
  int
  getAllSize() const
  {
    return stack_cont.getAllSize();
  }
  const Symmetry &
  getProdSym(int ip) const
  {
    return syms[ip];
  }
  bool
  isZero(const IntSequence &istacks) const
  {
    TL_RAISE_IF(istacks.size() != dimen(),
                "Wrong istacks coordinates for StackProduct::isZero");

    bool res = false;
    int i = 0;
    while (i < dimen() && !(res = stack_cont.isZero(istacks[i], syms[i])))
      i++;
    return res;
  }

  itype
  getType(int is, int ip) const
  {
    TL_RAISE_IF(is < 0 || is >= stack_cont.numStacks(),
                "Wrong index to stack in StackProduct::getType");
    TL_RAISE_IF(ip < 0 || ip >= dimen(),
                "Wrong index to stack container in StackProduct::getType");
    return stack_cont.getType(is, syms[ip]);
  }

625
  const _Ttype &
626
627
628
629
630
  getMatrix(int is, int ip) const
  {
    return stack_cont.getMatrix(is, syms[ip]);
  }

631
632
  std::vector<std::unique_ptr<Vector>>
  createPackedColumns(const IntSequence &coor, IntSequence &iu) const
633
634
  {
    TL_RAISE_IF(iu.size() != dimen(),
635
                "Wrong storage length for unit flags in StackProduct::createPackedColumns");
636
    TL_RAISE_IF(coor.size() != per.size(),
637
                "Wrong size of index coor in StackProduct::createPackedColumns");
638
    IntSequence perindex(coor.size());
639
    std::vector<std::unique_ptr<Vector>> vs;
640
641
642
643
644
    per.apply(coor, perindex);
    int off = 0;
    for (int i = 0; i < dimen(); i++)
      {
        IntSequence percoor(perindex, off, syms[i].dimen() + off);
645
        vs.push_back(stack_cont.createPackedColumn(syms[i], percoor, iu[i]));
646
647
        off += syms[i].dimen();
      }
648
    return vs;
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
  }

  int
  getSize(int is) const
  {
    return stack_cont.getStackSizes()[is];
  }

  int
  numMatrices(const IntSequence &istacks) const
  {
    TL_RAISE_IF(istacks.size() != dimen(),
                "Wrong size of stack coordinates in StackContainer::numMatrices");
    int ret = 0;
    int ip = 0;
    while (ip < dimen() && getType(istacks[ip], ip) == _Stype::matrix)
      {
        ret++;
        ip++;
      }
    return ret;
  }
};

673
674
/* Here we only inherit from Kronecker product KronProdAllOptim, only to
   allow for a constructor constructing from StackProduct. */
675

676
template<class _Ttype>
677
678
679
class KronProdStack : public KronProdAllOptim
{
public:
680
681
  using _Ptype = StackProduct<_Ttype>;
  using _Stype = StackContainerInterface<_Ttype>;
682

683
  /* Here we construct KronProdAllOptim from StackContainer and given
684
685
686
     selections of stack items from stack containers in the product. We
     only decide whether to insert matrix, or unit matrix.

687
688
     At this point, we do not call KronProdAllOptim::optimizeOrder(), so
     the KronProdStack behaves like KronProdAll (i.e. no optimization
689
690
691
692
693
694
695
696
697
698
     is done). */

  KronProdStack(const _Ptype &sp, const IntSequence &istack)
    : KronProdAllOptim(sp.dimen())
  {
    TL_RAISE_IF(sp.dimen() != istack.size(),
                "Wrong stack product dimension for KronProdStack constructor");

    for (int i = 0; i < sp.dimen(); i++)
      {
699
        TL_RAISE_IF(sp.getType(istack[i], i) == _Stype::itype::zero,
700
                    "Attempt to construct KronProdStack from zero matrix");
701
        if (sp.getType(istack[i], i) == _Stype::itype::unit)
702
          setUnit(i, sp.getSize(istack[i]));
703
        if (sp.getType(istack[i], i) == _Stype::itype::matrix)
704
          {
705
706
            const TwoDMatrix &m = sp.getMatrix(istack[i], i);
            TL_RAISE_IF(m.nrows() != sp.getSize(istack[i]),
707
                        "Wrong size of returned matrix in KronProdStack constructor");
708
            setMat(i, m);
709
710
711
712
713
          }
      }
  }
};

714
class WorkerFoldMAADense : public sthread::detach_thread
715
716
717
718
719
720
721
{
  const FoldedStackContainer &cont;
  Symmetry sym;
  const FGSContainer &dense_cont;
  FGSTensor &out;
public:
  WorkerFoldMAADense(const FoldedStackContainer &container,
722
                     Symmetry s,
723
724
                     const FGSContainer &dcontainer,
                     FGSTensor &outten);
725
  void operator()(std::mutex &mut) override;
726
727
};

728
class WorkerFoldMAASparse1 : public sthread::detach_thread
729
730
731
732
733
734
735
736
{
  const FoldedStackContainer &cont;
  const FSSparseTensor &t;
  FGSTensor &out;
  IntSequence coor;
public:
  WorkerFoldMAASparse1(const FoldedStackContainer &container,
                       const FSSparseTensor &ten,
737
                       FGSTensor &outten, IntSequence c);
738
  void operator()(std::mutex &mut) override;
739
740
};

741
class WorkerFoldMAASparse2 : public sthread::detach_thread
742
743
744
745
746
747
748
749
{
  const FoldedStackContainer &cont;
  const FSSparseTensor &t;
  FGSTensor &out;
  IntSequence coor;
public:
  WorkerFoldMAASparse2(const FoldedStackContainer &container,
                       const FSSparseTensor &ten,
750
                       FGSTensor &outten, IntSequence c);
751
  void operator()(std::mutex &mut) override;
752
753
};

754
class WorkerFoldMAASparse4 : public sthread::detach_thread
755
756
757
758
759
760
761
762
{
  const FoldedStackContainer &cont;
  const FSSparseTensor &t;
  FGSTensor &out;
  IntSequence coor;
public:
  WorkerFoldMAASparse4(const FoldedStackContainer &container,
                       const FSSparseTensor &ten,
763
                       FGSTensor &outten, IntSequence c);
764
  void operator()(std::mutex &mut) override;
765
766
};

767
class WorkerUnfoldMAADense : public sthread::detach_thread
768
769
770
771
772
773
774
{
  const UnfoldedStackContainer &cont;
  Symmetry sym;
  const UGSContainer &dense_cont;
  UGSTensor &out;
public:
  WorkerUnfoldMAADense(const UnfoldedStackContainer &container,
775
                       Symmetry s,
776
777
                       const UGSContainer &dcontainer,
                       UGSTensor &outten);
778
  void operator()(std::mutex &mut) override;
779
780
};

781
class WorkerUnfoldMAASparse1 : public sthread::detach_thread
782
783
784
785
786
787
788
789
{
  const UnfoldedStackContainer &cont;
  const FSSparseTensor &t;
  UGSTensor &out;
  IntSequence coor;
public:
  WorkerUnfoldMAASparse1(const UnfoldedStackContainer &container,
                         const FSSparseTensor &ten,
790
                         UGSTensor &outten, IntSequence c);
791
  void operator()(std::mutex &mut) override;
792
793
};

794
class WorkerUnfoldMAASparse2 : public sthread::detach_thread
795
796
797
798
799
800
801
802
{
  const UnfoldedStackContainer &cont;
  const FSSparseTensor &t;
  UGSTensor &out;
  IntSequence coor;
public:
  WorkerUnfoldMAASparse2(const UnfoldedStackContainer &container,
                         const FSSparseTensor &ten,
803
                         UGSTensor &outten, IntSequence c);
804
  void operator()(std::mutex &mut) override;
805
806
807
};

#endif