MVector.cpp

Go to the documentation of this file.

/*
 *  Copyright (c) 2015, Chris Rogers
 *  All rights reserved.
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *  1. Redistributions of source code must retain the above copyright notice,
 *     this list of conditions and the following disclaimer.
 *  2. Redistributions in binary form must reproduce the above copyright notice,
 *     this list of conditions and the following disclaimer in the documentation
 *     and/or other materials provided with the distribution.
 *  3. Neither the name of STFC nor the names of its contributors may be used to
 *     endorse or promote products derived from this software without specific
 *     prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 */

#include "Fields/Interpolation/MVector.h"
#include "Fields/Interpolation/MMatrix.h"

namespace interpolation {

std::ostream& operator<<(std::ostream& out, m_complex  c)
{
  out << re(c) << " r " << im(c) << " i";
  return out;
}

std::istream& operator>>(std::istream& in,  m_complex& c)
{
  std::string dummy;
  in >> re(c) >> dummy >> im(c) >> dummy;
  return in;
}



template MVector<double>   ::MVector(size_t i);
template MVector<m_complex>::MVector(size_t i);

template MVector<double>   ::MVector(const MVector<double>&);
template MVector<m_complex>::MVector(const MVector<m_complex>&);
template MVector<double>   ::MVector(size_t i, double    value);
template MVector<m_complex>::MVector(size_t i, m_complex value);

template void MVector<double>::build_vector( const double* data_begin, const double* data_end );
template void MVector<m_complex>::build_vector( const m_complex* data_begin, const m_complex* data_end );

template std::ostream& operator<<(std::ostream& out, MVector<double>    v);
template std::ostream& operator<<(std::ostream& out, MVector<m_complex> v);
template std::istream& operator>>(std::istream& out, MVector<double>&    v);
template std::istream& operator>>(std::istream& out, MVector<m_complex>& v);
template MVector<double>    MVector<double>   ::sub(size_t n1, size_t n2) const;
template MVector<m_complex> MVector<m_complex>::sub(size_t n1, size_t n2) const;


template <typename Tmplt>
MVector<Tmplt>::MVector( size_t i ) : _vector(NULL)
{
  build_vector(i);
}


template <typename Tmplt>
MVector<Tmplt>::MVector( const MVector<Tmplt>& mv) : _vector(NULL)
{ *this = mv; }


template <typename Tmplt>
MVector<Tmplt>::MVector( size_t size, Tmplt  value ) : _vector(NULL)
{
  build_vector(size);
  for(size_t i=0; i<size; i++) operator()(i+1) = value;
}


template <>
void MVector<double>::build_vector   ( size_t size )
{
  if(_vector != NULL) gsl_vector_free((gsl_vector*)_vector);
  _vector = gsl_vector_alloc(size);
}

template <>
void MVector<m_complex>::build_vector( size_t size )
{
  if(_vector != NULL) gsl_vector_complex_free((gsl_vector_complex*)_vector);
  _vector = gsl_vector_complex_alloc(size);
}

template <class Tmplt>
void MVector<Tmplt>::build_vector   ( const Tmplt* data_begin, const Tmplt* data_end )
{
  build_vector(data_end - data_begin);
  for(size_t i=0; i<num_row(); i++) operator()(i+1) = data_begin[i];
}


template <class Tmplt>
size_t MVector<Tmplt>::num_row() const
{ if(_vector != NULL) return ((gsl_vector*)_vector)->size; else return 0;}
template size_t MVector<double>   ::num_row() const;
template size_t MVector<m_complex>::num_row() const;

template <>
const double& MVector<double>::operator()(const size_t i) const
{return *gsl_vector_ptr( (gsl_vector*)_vector, i-1);}
template <>
const m_complex& MVector<m_complex>::operator()(const size_t i) const
{return *gsl_vector_complex_ptr((gsl_vector_complex*)_vector, i-1);}
template <>
double& MVector<double>::operator()(const size_t i)
{return *gsl_vector_ptr( (gsl_vector*)_vector, i-1);}
template <>
m_complex& MVector<m_complex>::operator()(const size_t i)
{return *gsl_vector_complex_ptr((gsl_vector_complex*)_vector, i-1);}


template <class Tmplt>
MMatrix<Tmplt> MVector<Tmplt>::T() const
{
  MMatrix<Tmplt> mat = MMatrix<Tmplt>(1,num_row());
  for(size_t i=1; i<=num_row(); i++)
    mat(1, i) = this->operator()(i);
  return mat;
}
template MMatrix<double>    MVector<double>   ::T() const;
template MMatrix<m_complex> MVector<m_complex>::T() const;

template <class Tmplt> bool operator==(const MVector<Tmplt>& c1, const MVector<Tmplt>& c2)
{
  if(c1.num_row() != c2.num_row())
    return false;
  for(size_t i=0; i<c1.num_row(); i++)
    if(c1(i+1) != c2(i+1)) return false;
  return true;
}
template bool operator==(const MVector<double>&    c1, const MVector<double>& c2);
template bool operator==(const MVector<m_complex>& c1, const MVector<m_complex>& c2);

template <>
MVector<double>& MVector<double>::operator= (const MVector<double>& mv)
{
  if (&mv == this) return *this;
  delete_vector();
  if(!mv._vector) { _vector = NULL; return *this; }
  _vector = gsl_vector_alloc( mv.num_row() );
  gsl_vector_memcpy((gsl_vector*)_vector, (const gsl_vector*)mv._vector);
  return *this;
}

template <>
MVector<m_complex>& MVector<m_complex>::operator= (const MVector<m_complex>& mv)
{
  if (&mv == this) return *this;
  delete_vector();
  if(!mv._vector) { _vector = NULL; return *this; }
  _vector = gsl_vector_complex_alloc( mv.num_row() );
  gsl_vector_complex_memcpy((gsl_vector_complex*)_vector, (const gsl_vector_complex*)mv._vector);
  return *this;
}

template <class Tmplt> std::ostream& operator<<(std::ostream& out, MVector<Tmplt> v)
{
  out << v.num_row() << "\n";
  for(size_t i=0; i<v.num_row(); i++) out << "  " << v(i+1) << "\n";
  return out;
}

template <class Tmplt> std::istream& operator>>(std::istream& in, MVector<Tmplt>& v)
{
  size_t n;
  in >> n;
  v = MVector<Tmplt>(n);
  for(size_t i=1; i<=v.num_row(); i++) in >> v(i);
  return in;
}

const gsl_vector*         MVector_to_gsl(const MVector<double>&      vd)
{return vd.get_vector(vd);}
const gsl_vector_complex* MVector_to_gsl(const MVector<gsl_complex>& vc)
{return vc.get_vector(vc);}

template <class Tmplt>
MVector<Tmplt> MVector<Tmplt>::sub(size_t n1, size_t n2) const
{
  MVector<Tmplt> temp(n2-n1+1);
  for(size_t i=n1; i<=n2; i++) temp(i-n1+1) = operator()(i);
  return temp;
}

MVector<m_complex> complex(MVector<double> real)
{
  MVector<m_complex> c(real.num_row());
  for(size_t i=1; i<=real.num_row(); i++) c(i) = m_complex_build(real(i));
  return c;
}

MVector<m_complex> complex(MVector<double> real, MVector<double> im)
{
  MVector<m_complex> c(real.num_row());
  for(size_t i=1; i<=real.num_row(); i++) c(i) = m_complex_build(real(i), im(i));
  return c;
}

MVector<double>    re     (MVector<m_complex> c)
{
  MVector<double> d(c.num_row());
  for(size_t i=1; i<=c.num_row(); i++) d(i) = re( c(i) );
  return d;
}

MVector<double>    im     (MVector<m_complex> c)
{
  MVector<double> d(c.num_row());
  for(size_t i=1; i<=c.num_row(); i++) d(i) = im( c(i) );
  return d;
}

}