Home     Back     Up Browse Files

Content

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.

#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.

#include <cmath>
#include <iostream>
#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
Home     Back     Up