gslpp_matrix_complex.cpp

Go to the documentation of this file.

/* 
 * Copyright (C) 2012 HEPfit Collaboration
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 */

#include <math.h>
#ifndef __GSL_BLAS_H__
#include <gsl/gsl_blas.h>
#endif
#ifndef __GSL_LINALG_H__
#include <gsl/gsl_linalg.h>
#endif
#ifndef __GSL_EIGEN_H__
#include <gsl/gsl_eigen.h>
#endif
#ifndef GSLPP_MATRIX_COMPLEX_H
#include "gslpp_matrix_complex.h"
#endif

namespace gslpp
{
  matrix<complex>::matrix(const size_t& size_i, const size_t& size_j, const complex& z)
  {
    _matrix = gsl_matrix_complex_alloc(size_i, size_j);
    gsl_matrix_complex_set_all(_matrix, z.as_gsl_type());
  }

  matrix<complex>::matrix(const size_t& size_i, const complex& z)
  {
    _matrix = gsl_matrix_complex_alloc(size_i, size_i);
    gsl_matrix_complex_set_all(_matrix, z.as_gsl_type());
  }

  matrix<complex>::matrix(const size_t& size_i, const size_t& size_j, const double& a)
  {
    complex z(a);
    _matrix = gsl_matrix_complex_alloc(size_i, size_j);
    gsl_matrix_complex_set_all(_matrix, z.as_gsl_type());
  }

  matrix<complex>::matrix(const size_t& size_i, const double& a)
  {
    complex z(a);
    _matrix = gsl_matrix_complex_alloc(size_i, size_i);
    gsl_matrix_complex_set_all(_matrix, z.as_gsl_type());
  }
  matrix<complex>::matrix(const matrix<complex>& m)
  {
    _matrix = gsl_matrix_complex_alloc(m.size_i(), m.size_j());
    gsl_matrix_complex_memcpy(_matrix, m.as_gsl_type_ptr());
  }

  matrix<complex>::matrix(const matrix<double>& m)
  {
    size_t i, j, size_i, size_j;
    size_i=m.size_i();
    size_j=m.size_j();
    _matrix = gsl_matrix_complex_alloc(size_i, size_j);
    for(i=0; i<size_i; ++i)
      for(j=0; j<size_j; ++j)
        gsl_matrix_complex_set(_matrix, i, j, m(i,j)+0.*complex::i());
  }

  matrix<complex>::matrix(const vector<complex>& v)
  {
    size_t i, size = v.size();
    complex z(0,0,false);
    _matrix = gsl_matrix_complex_alloc(size, size);
    gsl_matrix_complex_set_all(_matrix, z.as_gsl_type());
    for(i=0; i<size; ++i)
        gsl_matrix_complex_set(_matrix, i, i, v(i));
  }

  matrix<complex>::matrix(const vector<double>& v)
  {
    size_t i, size = v.size();
    complex z(0,0,false);
    _matrix = gsl_matrix_complex_alloc(size, size);
    gsl_matrix_complex_set_all(_matrix, z.as_gsl_type());
    for(i=0; i<size; ++i)
        gsl_matrix_complex_set(_matrix, i, i, v(i)+0.*complex::i());
  }

  matrix<complex>::matrix(const gsl_matrix_complex& m)
  {
    _matrix = gsl_matrix_complex_alloc(m.size1, m.size2);
    gsl_matrix_complex_memcpy(_matrix, &m);
  }

  matrix<complex>::matrix(const gsl_matrix_complex* m)
  {
    _matrix = gsl_matrix_complex_alloc(m->size1, m->size2);
    gsl_matrix_complex_memcpy(_matrix, m);
  }

  matrix<complex>::~matrix()
  {
    gsl_matrix_complex_free(_matrix);
  }

  const complex matrix<complex>::operator()(const size_t& i, const size_t& j) const
  {
    gsl_complex *x = gsl_matrix_complex_ptr(_matrix, i, j);
    return complex(x);
  }

//  complex& matrix<complex>::operator()(const size_t& i, const size_t& j)
//  {
//    gsl_complex *x = gsl_matrix_complex_ptr(_matrix, i, j);
//    return complex(x);
//  }

  matrix<complex>& matrix<complex>::operator=(const matrix<complex>& m)
  {
      if(size_i()==m.size_i() && size_j()==m.size_j())
          gsl_matrix_complex_memcpy(_matrix, m.as_gsl_type_ptr());
      else
      {
          std::cout << "\n ** Wrong size assign in matrix<complex> operator =" << std::endl;
          exit(EXIT_FAILURE);
      }
    return *this;
  }

  void matrix<complex>::assign(const size_t& i, const size_t& j, const complex& z)
  {
    gsl_complex *x = gsl_matrix_complex_ptr(_matrix, i, j);
    *x = z.as_gsl_type();
  }

  void matrix<complex>::assign(const size_t& i, const size_t& j,const double& a)
  {
    gsl_complex *x = gsl_matrix_complex_ptr(_matrix, i, j);
    *x = complex(a).as_gsl_type();
  }
  
  void matrix<complex>::assignre(const size_t& i, const size_t& j,const double& a)
  {
    gsl_complex *x = gsl_matrix_complex_ptr(_matrix, i, j);
    GSL_SET_REAL(x,a);
  }
  
  void matrix<complex>::assignim(const size_t& i, const size_t& j,const double& a)
  {
    gsl_complex *x = gsl_matrix_complex_ptr(_matrix, i, j);
    GSL_SET_IMAG(x,a);
  }
  
  void matrix<complex>::assign(const size_t& i, const size_t& j, const matrix<complex>& z)
  {
      size_t ki,kj;
      gsl_complex *x;
      if(i+z.size_i() <= size_i() && j+z.size_j() <= size_j())
          for(ki=i;ki<i+z.size_i();ki++)
              for(kj=j;kj<j+z.size_j();kj++)
              {
                  x = gsl_matrix_complex_ptr(_matrix, ki, kj);
                  *x = z(ki-i,kj-j).as_gsl_type();
              }
      else
      {
          std::cout << "\n ** Wrong size assign in matrix<complex> assign submatrix" << std::endl;
          exit(EXIT_FAILURE);
      }
  }

  void matrix<complex>::assign(const size_t& i, const size_t& j, const matrix<double>& a)
  {
      matrix<complex> z(a);
      assign(i,j,z);
  }

  size_t matrix<complex>::size_i() const
  {
    return _matrix->size1;
  }

  size_t matrix<complex>::size_j() const
  {
    return _matrix->size2;
  }

  matrix<complex> matrix<complex>::Id(size_t size)
  {
    return matrix<double>::Id(size) * (1.+0.*complex::i());
  }

  matrix<complex> matrix<complex>::transpose() const
  {
    matrix<complex> m1(size_j(), size_i(), 0.);
    if (gsl_matrix_complex_transpose_memcpy(m1.as_gsl_type_ptr(), _matrix))
    {
      std::cout << "\n ** Error in matrix<complex> transpose" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }

  matrix<complex> matrix<complex>::hconjugate() const
  {
    matrix<complex> m1(*this);
    gsl_complex z1, z2;

    GSL_SET_COMPLEX(&z1,1.,0.);
    GSL_SET_COMPLEX(&z2,0.,0.);

    if(gsl_blas_zgemm(CblasConjTrans,CblasNoTrans, z1, _matrix,
       Id(size_j()).as_gsl_type_ptr(), z2, m1.as_gsl_type_ptr()) != 0)
    {
      std::cout << "\n ** Error in matrix<complex> conjugate" << std::endl;
      exit(EXIT_FAILURE);
    }

    return m1;
  }

  matrix<complex> matrix<complex>::inverse() const
  {
    matrix<complex> m1(_matrix);
    matrix<complex> m2(_matrix);
    int signum;
    gsl_permutation *p;

    if (size_j() != size_i())
    {
      std::cout << "\n ** Size mismatch in matrix<complex> inverse" << std::endl;
      exit(EXIT_FAILURE);
    }

    if ((p = gsl_permutation_alloc(size_i())) == NULL)
    {
      std::cout << "\n ** Error in matrix<complex> inverse" << std::endl;
      exit(EXIT_FAILURE);
    }

    if(gsl_linalg_complex_LU_decomp (m1.as_gsl_type_ptr(), p, &signum))
    {
      std::cout << "\n ** Error in matrix<complex> inverse" << std::endl;
      gsl_permutation_free(p);
      exit(EXIT_FAILURE);
    }

    if(gsl_linalg_complex_LU_invert(m1.as_gsl_type_ptr(), p,
       m2.as_gsl_type_ptr()))
    {
      std::cout << "\n ** Error in matrix<complex> inverse" << std::endl;
      gsl_permutation_free(p);
      exit(EXIT_FAILURE);
    }
    gsl_permutation_free(p);
    return m2;
  }
  
    matrix<double> matrix<complex>::real() const 
    {
        matrix<double> res(size_i(), size_j());
        for(size_t i = 0; i < size_i(); i++)  {
            gsl_vector_complex_view tmp = gsl_matrix_complex_row(_matrix, i);
            gsl_vector_view tmpd = gsl_vector_complex_real(&tmp.vector);
            gsl_matrix_set_row(res.as_gsl_type_ptr(),i,&tmpd.vector);
        }
       return res;
    }

    matrix<double> matrix<complex>::imag() const 
    {
        matrix<double> res(size_i(), size_j());
        for(size_t i = 0; i < size_i(); i++)  {
            gsl_vector_complex_view tmp = gsl_matrix_complex_row(_matrix, i);
            gsl_vector_view tmpd = gsl_vector_complex_imag(&tmp.vector);
            gsl_matrix_set_row(res.as_gsl_type_ptr(),i,&tmpd.vector);
        }
       return res;
    }
  
  void matrix<complex>::eigensystem(matrix<complex> &U, vector<double> &S) const
  {
      matrix<complex> m(*this);

      gsl_eigen_hermv_workspace *ws;

      ws = gsl_eigen_hermv_alloc(size_i());

      gsl_eigen_hermv(m.as_gsl_type_ptr(), S.as_gsl_type_ptr(), U.as_gsl_type_ptr(), ws);

      gsl_eigen_hermv_sort(S.as_gsl_type_ptr(), U.as_gsl_type_ptr(),
             GSL_EIGEN_SORT_VAL_ASC);

      gsl_eigen_hermv_free(ws);
  }

  void matrix<complex>::singularvalue(matrix<complex> &U, matrix<complex> &V,
          vector<double> &S) const
  {
     size_t i;
     matrix<complex> m(*this);
     vector<complex> v1(size_i(),0.);

     m = (*this)*hconjugate();
     m.eigensystem(U, S);
     m = hconjugate()*(*this);
     m.eigensystem(V, S);
     m = U.hconjugate()*(*this)*V;
     for(i=0; i<m.size_i(); i++)
     {
         v1.assign(i,complex(1.,-m(i,i).arg(),true));
         S(i)=m(i,i).abs();
     }
     V=V*matrix<complex>(v1);
  }

  gsl_matrix_complex* matrix<complex>::as_gsl_type_ptr() const
  {
    return _matrix;
  }

  gsl_matrix_complex& matrix<complex>::as_gsl_type()
  {
    return const_cast<gsl_matrix_complex&>(*_matrix);
  }

  const gsl_matrix_complex& matrix<complex>::as_gsl_type() const
  {
    return const_cast<gsl_matrix_complex&>(*_matrix);
  }

//  bool matrix<complex>::is_equal(const matrix<complex>& m1, const matrix<complex>& m2){
//      return gsl_matrix_complex_equal(m1.as_gsl_type_ptr(),m2.as_gsl_type_ptr());
//  }
  
  matrix<complex> matrix<complex>::operator-() const
  {
    matrix<complex> m1(_matrix);
    gsl_complex z1;
    GSL_SET_COMPLEX(&z1,-1.,0.);
    if (gsl_matrix_complex_scale(m1.as_gsl_type_ptr(), z1))
    {
      std::cout << "\n ** Error in matrix<complex> unary -" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }

  matrix<complex> matrix<complex>::operator+(const matrix<complex>& m) const
  {
    matrix<complex> m1(_matrix);
    if (gsl_matrix_complex_add(m1.as_gsl_type_ptr(), m.as_gsl_type_ptr()))
    {
      std::cout << "\n ** Error in matrix<complex> +" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }
  matrix<complex> matrix<complex>::operator-(const matrix<complex>& m) const 
  {
    matrix<complex> m1(_matrix);
    matrix<complex> m2 = -m;
    if (gsl_matrix_complex_add(m1.as_gsl_type_ptr(), m2.as_gsl_type_ptr()))
    {
      std::cout << "\n ** Error in matrix<complex> +" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }

  matrix<complex> matrix<complex>::operator*(const matrix<complex>& m) const 
  {
    matrix<complex> m1(_matrix);
    gsl_complex z1, z2;

    if (size_j() != m.size_i())
    {
      std::cout << "\n ** Size mismatch in matrix<complex> *" << std::endl;
      exit(EXIT_FAILURE);
    }
    GSL_SET_COMPLEX(&z1,1.,0.);
    GSL_SET_COMPLEX(&z2,0.,0.);
    if(gsl_blas_zgemm(CblasNoTrans,CblasNoTrans, z1, _matrix,
       m.as_gsl_type_ptr(), z2, m1.as_gsl_type_ptr()))
    {
      std::cout << "\n ** Error in matrix<complex> *" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }

  matrix<complex> matrix<complex>::operator+(const matrix<double>& m) const 
  {
    matrix<complex> m1(_matrix);
    matrix<complex> m2(m);
    if (gsl_matrix_complex_add(m1.as_gsl_type_ptr(), m2.as_gsl_type_ptr()))
    {
      std::cout << "\n ** Error in matrix<complex> +" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }
  matrix<complex> matrix<complex>::operator-(const matrix<double>& m) const
  {
    matrix<complex> m1(_matrix);
    matrix<complex> m2(-m);
    if (gsl_matrix_complex_add(m1.as_gsl_type_ptr(), m2.as_gsl_type_ptr()))
    {
      std::cout << "\n ** Error in matrix<complex> +" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }

  matrix<complex> matrix<complex>::operator*(const matrix<double>& m) const
  {
    matrix<complex> m1(m);
    return *this * m1;
  }

  vector<complex> matrix<complex>::operator*(const vector<complex>& v) const
  {
    gsl_complex z1, z2;
    vector<complex> v1(size_i(),0.);

    if (size_j() != v.size())
    {
      std::cout << "\n ** Size mismatch in matrix<complex> * (vector complex)"
          << std::endl;
      exit(EXIT_FAILURE);
    }
    GSL_SET_COMPLEX(&z1, 1., 0.);
    GSL_SET_COMPLEX(&z2, 0., 0.);
    if(gsl_blas_zgemv(CblasNoTrans, z1, _matrix, v.as_gsl_type_ptr(),
       z2, v1.as_gsl_type_ptr()))
    {
      std::cout << "\n ** Error in matrix<complex> * (vector complex)"
          << std::endl;
      exit(EXIT_FAILURE);
    }
    return v1;
  }

  vector<complex> matrix<complex>::operator*(const vector<double>& v) const
  {
    vector<complex> v1(v);
    return *this * v1;
  }

  matrix<complex>& matrix<complex>::operator+=(const matrix<complex>& m)
  {
    *this = *this + m;
    return *this;
  }
  matrix<complex>& matrix<complex>::operator-=(const matrix<complex>& m)
  {
    *this = *this - m;
    return *this;
  }

  matrix<complex>& matrix<complex>::operator*=(const matrix<complex>& m)
  {
    if (!((size_i() == size_j()) && (size_i() == m.size_i()) &&
           (size_j() == m.size_j())))
    {
      std::cout << "\n ** Size mismatch in matrix<complex> *= (matrix)"
          << std::endl;
      exit(EXIT_FAILURE);
    }
    *this = *this * m;
    return *this;
  }

  matrix<complex> matrix<complex>::operator+(const complex& z) const
  {
    matrix<complex> m1(_matrix);
    if (gsl_matrix_complex_add_constant(_matrix, z.as_gsl_type()))
    {
      std::cout << "\n ** Error in matrix<complex> + (complex)" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }

  matrix<complex> matrix<complex>::operator-(const complex& z) const
  {
    matrix<complex> m1(_matrix);
    if (gsl_matrix_complex_add_constant(_matrix, (-z).as_gsl_type()))
    {
      std::cout << "\n ** Error in matrix<complex> - (complex)" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }

  matrix<complex> matrix<complex>::operator*(const complex& z) const
  {
    matrix<complex> m1(_matrix);
    if (gsl_matrix_complex_scale(m1.as_gsl_type_ptr(), z.as_gsl_type()))
    {
      std::cout << "\n ** Error in matrix<complex> * (complex)" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }

  matrix<complex> matrix<complex>::operator/(const complex& z) const
  {
    matrix<complex> m1(_matrix);
    if (gsl_matrix_complex_scale(m1.as_gsl_type_ptr(), z.inverse().as_gsl_type()))
    {
      std::cout << "\n ** Error in matrix<complex> / (complex)" << std::endl;
      exit(EXIT_FAILURE);
    }
    return m1;
  }
  
  matrix<complex>& matrix<complex>::operator+=(const complex& z)
  {
    *this = *this + z;
    return *this;
  }
  matrix<complex>& matrix<complex>::operator-=(const complex& z)
  {
    *this = *this - z;
    return *this;
  }
  matrix<complex>& matrix<complex>::operator*=(const complex& z)
  {
    *this = *this * z;
    return *this;
  }
  matrix<complex>& matrix<complex>::operator/=(const complex& z)
  {
    *this = *this / z;
    return *this;
  }

  matrix<complex> matrix<complex>::operator+(const double& a) const
  {
    complex z(a);
    return *this + z;
  }

  matrix<complex> matrix<complex>::operator-(const double& a) const
  {
    complex z(a);
    return *this - z;
  }

  matrix<complex> matrix<complex>::operator*(const double& a) const
  {
    complex z(a);
    return *this * z;
  }

  matrix<complex> matrix<complex>::operator/(const double& a) const
  {
    complex z(a);
    return *this / z;
  }
  matrix<complex>& matrix<complex>::operator+=(const double& a)
  {
    *this = *this + a;
    return *this;
  }
  matrix<complex>& matrix<complex>::operator-=(const double& a)
  {
    *this = *this - a;
    return *this;
  }
  matrix<complex>& matrix<complex>::operator*=(const double& a)
  {
    *this = *this * a;
    return *this;
  }
  matrix<complex>& matrix<complex>::operator/=(const double& a)
  {
    *this = *this / a;
    return *this;
  }
  std::ostream& operator<<(std::ostream& output, const matrix<complex>& m)
  {
    size_t i,j;
    for (i=0; i<m.size_i(); i++)
    {
      output << std::endl;
      output << "\t(";
      for (j=0; j<m.size_j()-1; j++)
        output << m(i,j) << ",";
      output << m(i,j) << ")";
    }
    return output;
  }

  matrix<complex> operator+(const matrix<double> m1, const matrix<complex> m2) {
      return m2 + m1;
  }

  matrix<complex> operator-(const matrix<double> m1, const matrix<complex> m2) {
      return -m2 + m1;
  }

  matrix<complex> operator+(const complex& z, const matrix<complex> m)
  {
    return m + z;
  }

  matrix<complex> operator-(const complex& z, const matrix<complex> m)
  {
    return -m + z;
  }

  matrix<complex> operator*(const complex& z, const matrix<complex> m)
  {
    return m * z;
  }

  vector<complex> operator*(const vector<complex>& v, const matrix<complex> m) 
  {
    return m.transpose() * v;
  }

  vector<complex> operator*(const vector<double>& v, const matrix<complex> m)
  {
    return m.transpose() * v;
  }

  matrix<complex> operator+(const double& a, const matrix<complex> m)
  {
    return m + a;
  }

  matrix<complex> operator-(const double& a, const matrix<complex> m)
  {
    return -m + a;
  }

  matrix<complex> operator*(const double& a, const matrix<complex> m)
  {
    return m * a;
  }
}