[ublas] vector function

Subject: [ublas] vector function
From: Heiko Bauke (heiko.bauke_at_[hidden])
Date: 2009-03-13 03:40:00

Hi,

I have a question similar to the question by Denis Taniguchi.

I am working on a project where special (sparse) matrix vector
multiplications are required. The matrices are fixed 4x4 matrices (the
Dirac matrices). At the moment matrix vector multiplications are
implemented as functions, e.g.

template<typename T>
blas::vector<T> mult_alpha_x(const blas::vector<T> &v) {
  blas::vector<T> v_new(4);
  v_new(0)=v(3);
  v_new(1)=v(2);
  v_new(2)=v(1);
  v_new(3)=v(0);
  return v_new;
}

I would like to take advantage of expression templates and implement
these functions such that they take a vector expression as argument
and return a vector expression as result. How can I implement such 
functions?

        Heiko

#include <cstdlib>
#include <iostream>
#include <complex>
#include <boost/numeric/ublas/vector.hpp>

namespace blas=boost::numeric::ublas;

typedef std::complex<double> complex;

template<typename T>
blas::vector<T> mult_alpha_x(const blas::vector<T> &v) {
  blas::vector<T> v_new(4);
  v_new(0)=v(3);
  v_new(1)=v(2);
  v_new(2)=v(1);
  v_new(3)=v(0);
  return v_new;
}

template<typename T>
blas::vector<T> mult_alpha_y(const blas::vector<T> &v) {
  blas::vector<T> v_new(4);
  v_new(0)=complex(0, -1)*v(3);
  v_new(1)=complex(0,  1)*v(2);
  v_new(2)=complex(0, -1)*v(1);
  v_new(3)=complex(0,  1)*v(0);
  return v_new;
}

template<typename T>
blas::vector<T> mult_alpha_z(const blas::vector<T> &v) {
  blas::vector<T> v_new(4);
  v_new(0)=v(2);
  v_new(1)=-v(3);
  v_new(2)=v(0);
  v_new(3)=-v(1);
  return v_new;
}

int main() {
  blas::vector<complex> v(4);
  v(0)=1;  v(1)=2;
  v(2)=3;  v(3)=4;
  mult_alpha_x(v);
  mult_alpha_y(v);
  mult_alpha_z(v);
  return EXIT_SUCCESS;
}

-- 
-- Jäten ist Zensur an der Natur.
-- (Oskar Kokoschka)
-- Cluster Computing @ http://www.clustercomputing.de
--       Heiko Bauke @ http://www.mpi-hd.mpg.de/personalhomes/bauke