permutation.h (flens/examples/permutation.h)

      1
      2
      3
      4
<doc 5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
     56
     57
     58
     59
     60
     61
     62
     63
     64
     65
     66
     67
     68
     69
     70
     71
     72
     73
     74
     75
     76
     77
     78
     79

#include <flens/flens.h>

namespace flens {

struct Permutation
    : public GeneralMatrix<Permutation>
{
    typedef int   ElementType;
    typedef int   IndexType;

    explicit
    Permutation(int n)
        : vector(n)
    {
    }

    DenseVector<Array<int> >  vector;
};

namespace blas {

template <typename ALPHA, typename VX, typename BETA, typename VY>
void
mv(Transpose trans, const ALPHA &alpha,
   const Permutation &P, const DenseVector<VX> &x,
   const BETA &beta, DenseVector<VY> &y)
{
    if (y.length()==0) {
        y.resize(x.length(), x.firstIndex());
    }
    if (trans==NoTrans) {
        int ix0 = x.firstIndex();
        int iy0 = y.firstIndex();
        for (int i=1, ix=ix0, iy=iy0; i<=x.length(); ++i, ++ix, ++iy) {
            y(i) = beta*y(i) + alpha*x(P.vector(i));
        }
    } else {
        int ix0 = x.firstIndex();
        int iy0 = y.firstIndex();
        for (int i=1, ix=ix0, iy=iy0; i<=x.length(); ++i, ++ix, ++iy) {
            y(P.vector(i)) = beta*y(i) + alpha*x(i);
        }
    }
}

template <typename ALPHA, typename MA, typename BETA, typename MB>
void
mm(Transpose transP, Transpose transA, const ALPHA &alpha,
   const Permutation &P, const GeMatrix<MA> &A, const BETA &beta,
   GeMatrix<MB> &B)
{
    typedef typename MB::IndexType   IndexType;

    if (B.numRows()==0 || B.numCols()==0) {
        B.resize(A.numRows(), A.numCols(), A.firstRow(), A.firstCol());
    }

    const Underscore<IndexType> _;

    if (transA==NoTrans) {
        for (IndexType j=1; j<=A.numCols(); ++j) {
            const auto x = A(_,j);
            auto       y = B(_,j);
            mv(transP, alpha, P, x, beta, y);
        }
    } else {
        for (IndexType i=1; i<=A.numRows(); ++i) {
            const auto x = A(i,_);
            auto       y = B(_,i);
            mv(transP, alpha, P, x, beta, y);
        }
    }
}

} // namespace blas

} // namespace flens