Vector.cpp

Go to the documentation of this file.

/* 
 * @BEGIN LICENSE
 *
 * Copyright (C) 2014-2015  Ward Poelmans
 *
 * This file is part of v2DM-DOCI.
 * 
 * v2DM-DOCI is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Foobar is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
 *
 * @END LICENSE
 */

#include <iostream>
#include <fstream>
#include <algorithm>
#include <time.h>
#include <cmath>
#include <cstring>
#include <assert.h>

using std::endl;
using std::ostream;
using std::ofstream;
using std::ifstream;

#include "Vector.h"
#include "lapack.h"

using namespace doci2DM;

Vector::Vector(int n)
{
   this->n = n;

   vector.reset(new double[n]);
}

Vector::Vector(Matrix &matrix)
{
   //allocate
   this->n = matrix.gn();

   vector.reset(new double[n]);

   diagonalize(matrix);
}

Vector::Vector(const Vector &vec_copy)
{
   this->n = vec_copy.n;

   vector.reset(new double[n]);

   std::memcpy(vector.get(), vec_copy.vector.get(), n*sizeof(double));
}

Vector::Vector(Vector &&vec_copy)
{
   this->n = vec_copy.n;

   vector = std::move(vec_copy.vector);
}


Vector::~Vector()
{
}

Vector &Vector::operator=(const Vector &vec_copy)
{
   assert(vec_copy.n == n);
   std::memcpy(vector.get(), vec_copy.vector.get(), n*sizeof(double));

   return *this;
}

Vector &Vector::operator=(Vector &&vec_copy)
{
   assert(vec_copy.n == n);

   vector = std::move(vec_copy.vector);

   return *this;
}

Vector &Vector::operator=(double a)
{
   for(int i = 0;i < n;++i)
      vector[i] = a;

   return *this;
}

Vector &Vector::operator+=(const Vector &vector_pl)
{
   daxpy(1,vector_pl);

   return *this;
}

Vector &Vector::operator-=(const Vector &vector_pl)
{
   daxpy(-1,vector_pl);

   return *this;
}

Vector &Vector::daxpy(double alpha,const Vector &vector_pl)
{
   int inc = 1;

   daxpy_(&n,&alpha,vector_pl.vector.get(),&inc,vector.get(),&inc);

   return *this;
}

Vector &Vector::operator/=(double c)
{
   dscal(1.0/c);

   return *this;
}

Vector &Vector::operator*=(double c)
{
   dscal(c);

   return *this;
}

double &Vector::operator[](int i)
{
   assert(i<n);

   return vector[i];
}

double Vector::operator[](int i) const
{
   assert(i<n);

   return vector[i];
}

void Vector::diagonalize(Matrix &matrix)
{
   *this = matrix.diagonalize();
}

double *Vector::gVector()
{
   return vector.get();
}

const double *Vector::gVector() const
{
   return vector.get();
}

int Vector::gn() const
{
   return n;
}

double Vector::sum() const
{
   double ward = 0;

   for(int i = 0;i < n;++i)
      ward += vector[i];

   return ward;
} 

double Vector::trace() const
{
   return sum();
}

double Vector::log_product() const
{
   double ward = 0;

   for(int i = 0;i < n;++i)
      ward += log(vector[i]);

   return ward;
}

double Vector::ddot(const Vector &vector_i) const
{
   int inc = 1;

   return ddot_(&n,vector.get(),&inc,vector_i.vector.get(),&inc);
}

void Vector::dscal(double alpha)
{
   int inc = 1;

   dscal_(&n,&alpha,vector.get(),&inc);
}

void Vector::fill_Random()
{
   fill_Random(time(NULL));
}

void Vector::fill_Random(int seed)
{
   srand(seed);

   for(int i = 0;i < n;++i)
      vector[i] = (double) rand()/RAND_MAX;
}

namespace doci2DM 
{
   std::ostream &operator<<(std::ostream &output,const doci2DM::Vector &vector_p)
   {
      for(int i = 0;i < vector_p.gn();++i)
         output << i << "\t" << vector_p[i] << endl;

      return output;
   }
}

double Vector::min() const
{
   double min = vector[0];

   for(int i=1;i<n;i++)
      if(vector[i]<min)
         min = vector[i];

   return min;
}

double Vector::max() const
{
   double max = vector[n-1];

   for(int i=n-2;i>=0;i--)
      if(vector[i]>max)
         max = vector[i];

   return max;
}

void Vector::invert()
{
   for(int i = 0;i < n;++i)
      vector[i] = 1.0/vector[i];
}

void Vector::sqrt(int option)
{
   if(option == 1)
      for(int i = 0;i < n;++i)
         vector[i] = std::sqrt(vector[i]);
   else 
      for(int i = 0;i < n;++i)
         vector[i] = 1.0/std::sqrt(vector[i]);
}

void Vector::symmetrize()
{
}

void Vector::L_map(const Vector &a,const Vector &b)
{
   assert(a.n == n && b.n == n);

   for(int i = 0;i < n;++i)
      vector[i] = a[i] * b[i] * a[i];
}

Vector &Vector::mprod(const Vector &x,const Vector &y)
{
   assert(x.n == n && y.n == n);

   for(int i = 0;i < n;++i)
      vector[i] = x[i] * y[i];

   return *this;
}

void Vector::sort()
{
   std::sort(vector.get(), &vector.get()[n] );
}

void Vector::sep_pm(Vector &pos, Vector &neg)
{
   assert(pos.n == n && neg.n == n);
   pos = 0;
   neg = 0;

   for(int i=0;i<n;i++)
      if(vector[i] < 0)
         neg[i] = vector[i];
      else
         pos[i] = vector[i];
}

/* vim: set ts=3 sw=3 expandtab :*/