Defining a Permutation Matrix

In this example we show how you can define your own general matrix type and support the matrix-vector product and matrix-matrix product for it. The matrix type we define represents a permutation.

Our Permutation Matrix

The permutation is derived from GeneralMatrix and we define functions blas::mv and blas::mm for it.

flens/examples/permutation.h

#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

Example

That's it. Now we just use it.

flens/examples/permutation.cc

#include <cxxstd/cmath.h>

#include <cxxstd/iostream.h>

#include <flens/examples/permutation.h>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    typedef DenseVector<Array<double> >     RealDenseVector;

    typedef GeMatrix<FullStorage<double> >  RealGeMatrix;

    const int   n = 5;

    Permutation P(n);

    P.vector = 5, 3, 2, 1, 4;

    RealDenseVector                 x(5), y;

    x = 1, 2, 3, 4, 5;

    y = P*x;

    cout << "x = " << x << endl;

    cout << "y = P*x = " << y << endl;

    x = transpose(P)*y;

    cout << "x = P^T*y = " << x << endl;

    RealGeMatrix  I(n, n), B(n,n);

    I = 0;

    I.diag(0) = 1;

    B = P*I;

    cout << "B = P*I = " << B << endl;

}

Let's compile and run

$shell> cd flens/examples                                                      
$shell> g++ -Wall -std=c++11 -I../.. permutation.cc                            
$shell> ./a.out                                                                
x = 
            1              2              3              4              5 
y = P*x = 
            5              3              2              1              4 
x = P^T*y = 
            1              2              3              4              5 
B = P*I = 
            0             0             0             0             1 
            0             0             1             0             0 
            0             1             0             0             0 
            1             0             0             0             0 
            0             0             0             1             0