Skip to content
Snippets Groups Projects
Select Git revision
  • 38824dc1e5c3afd38b25bad8f04f9f527770f6c8
  • master default
  • fcst_shock_decomp
  • 4.5
  • clang+openmp
  • local_state_space_iteration_k
  • dynamic-striated
  • occbin
  • exo_steady_state
  • filter_initial_state
  • declare_vars_in_model_block
  • exceptions
  • rmExtraExo
  • julia
  • error_msg_undeclared_model_vars
  • static_aux_vars
  • slice
  • aux_func
  • penalty
  • 4.4
  • separateM_
  • 4.5.7
  • 4.5.6
  • 4.5.5
  • 4.5.4
  • 4.5.3
  • 4.5.2
  • 4.5.1
  • 4.5.0
  • 4.4.3
  • 4.4.2
  • 4.4.1
  • 4.4.0
  • 4.4-beta1
  • 4.3.3
  • 4.3.2
  • 4.3.1
  • 4.3.0
  • 4.2.5
  • 4.2.4
  • 4.2.3
41 results

ExprNode.hh

Blame
  • Forked from Dynare / dynare
    Source project has a limited visibility.
    Vector.cc 8.03 KiB
    // Copyright (C) 2004-2011, Ondra Kamenik
    
    #include "Vector.hh"
    #include "GeneralMatrix.hh"
    #include "SylvException.hh"
    
    #include <dynblas.h>
    
    #include <cstdio>
    #include <cstring>
    #include <cstdlib>
    #include <cmath>
    #include <algorithm>
    #include <limits>
    
    using namespace std;
    
    ZeroPad zero_pad;
    
    Vector::Vector(const Vector &v)
      : len(v.length()),  data(new double[len]), destroy(true)
    {
      copy(v.base(), v.skip());
    }
    
    Vector::Vector(const ConstVector &v)
      : len(v.length()),  data(new double[len]), destroy(true)
    {
      copy(v.base(), v.skip());
    }
    
    const Vector &
    Vector::operator=(const Vector &v)
    {
      if (this == &v)
        return *this;
    
      if (v.length() != length())
        {
          throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
        }
      if (s == v.s
          && (data <= v.data && v.data < data+len*s
              || v.data <= data && data < v.data+v.len*v.s)
          && (data-v.data) % s == 0)
        {
          printf("this destroy=%d, v destroy=%d, data-v.data=%lu, len=%d\n", destroy, v.destroy, (unsigned long) (data-v.data), len);
          throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors.");
        }
      copy(v.base(), v.skip());
      return *this;
    }
    
    const Vector &
    Vector::operator=(const ConstVector &v)
    {
      if (v.length() != length())
        {
          throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
        }
      if (v.skip() == 1 && skip() == 1 && (
                                           (base() < v.base() + v.length() && base() >= v.base())
                                           || (base() + length() < v.base() + v.length()
                                               && base() + length() > v.base())))
        {
          throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors.");
        }
      copy(v.base(), v.skip());
      return *this;
    }
    
    void
    Vector::copy(const double *d, int inc)
    {
      blas_int n = length();
      blas_int incy = skip();
      blas_int inc2 = inc;
      dcopy(&n, d, &inc2, base(), &incy);
    }
    
    Vector::Vector(Vector &v, int off, int l)
      : len(l), s(v.skip()), data(v.base()+off*v.skip()) 
    {
      if (off < 0 || off + length() > v.length())
        throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
    }
    
    Vector::Vector(const Vector &v, int off, int l)
      : len(l),  data(new double[len]), destroy(true)
    {
      if (off < 0 || off + length() > v.length())
        throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
      copy(v.base()+off*v.skip(), v.skip());
    }
    
    Vector::Vector(GeneralMatrix &m, int col)
      : len(m.numRows()),  data(&(m.get(0, col))) 
    {
    }
    
    Vector::Vector(int row, GeneralMatrix &m)
      : len(m.numCols()), s(m.getLD()), data(&(m.get(row, 0))) 
    {
    }
    
    bool
    Vector::operator==(const Vector &y) const
    {
      return ConstVector(*this) == y;
    }
    
    bool
    Vector::operator!=(const Vector &y) const
    {
      return ConstVector(*this) != y;
    }
    
    bool
    Vector::operator<(const Vector &y) const
    {
      return ConstVector(*this) < y;
    }
    
    bool
    Vector::operator<=(const Vector &y) const
    {
      return ConstVector(*this) <= y;
    }
    
    bool
    Vector::operator>(const Vector &y) const
    {
      return ConstVector(*this) > y;
    }
    
    bool
    Vector::operator>=(const Vector &y) const
    {
      return ConstVector(*this) >= y;
    }
    
    void
    Vector::zeros()
    {
      if (skip() == 1)
        {
          double *p = base();
          for (int i = 0; i < length()/ZeroPad::length;
               i++, p += ZeroPad::length)
            memcpy(p, zero_pad.getBase(), sizeof(double)*ZeroPad::length);
          for (; p < base()+length(); p++)
            *p = 0.0;
        }
      else
        {
          for (int i = 0; i < length(); i++)
            operator[](i) = 0.0;
        }
    }
    
    void
    Vector::nans()
    {
      for (int i = 0; i < length(); i++)
        operator[](i) = std::numeric_limits<double>::quiet_NaN();
    }
    
    void
    Vector::infs()
    {
      for (int i = 0; i < length(); i++)
        operator[](i) = std::numeric_limits<double>::infinity();
    }
    
    Vector::~Vector()
    {
      if (destroy)
        {
          delete [] data;
        }
    }
    
    void
    Vector::rotatePair(double alpha, double beta1, double beta2, int i)
    {
      double tmp = alpha*operator[](i) - beta1*operator[](i+1);
      operator[](i+1) = alpha*operator[](i+1) - beta2*operator[](i);
      operator[](i) = tmp;
    }
    
    void
    Vector::add(double r, const Vector &v)
    {
      add(r, ConstVector(v));
    }
    
    void
    Vector::add(double r, const ConstVector &v)
    {
      blas_int n = length();
      blas_int incx = v.skip();
      blas_int incy = skip();
      daxpy(&n, &r, v.base(), &incx, base(), &incy);
    }
    
    void
    Vector::add(const double *z, const Vector &v)
    {
      add(z, ConstVector(v));
    }
    
    void
    Vector::add(const double *z, const ConstVector &v)
    {
      blas_int n = length()/2;
      blas_int incx = v.skip();
      blas_int incy = skip();
      zaxpy(&n, z, v.base(), &incx, base(), &incy);
    }
    
    void
    Vector::mult(double r)
    {
      blas_int n = length();
      blas_int incx = skip();
      dscal(&n, &r, base(), &incx);
    }
    
    void
    Vector::mult2(double alpha, double beta1, double beta2,
                  Vector &x1, Vector &x2,
                  const Vector &b1, const Vector &b2)
    {
      x1.zeros();
      x2.zeros();
      mult2a(alpha, beta1, beta2, x1, x2, b1, b2);
    }
    
    void
    Vector::mult2a(double alpha, double beta1, double beta2,
                   Vector &x1, Vector &x2,
                   const Vector &b1, const Vector &b2)
    {
      x1.add(alpha, b1);
      x1.add(-beta1, b2);
      x2.add(alpha, b2);
      x2.add(-beta2, b1);
    }
    
    double
    Vector::getNorm() const
    {
      ConstVector v(*this);
      return v.getNorm();
    }
    
    double
    Vector::getMax() const
    {
      ConstVector v(*this);
      return v.getMax();
    }
    
    double
    Vector::getNorm1() const
    {
      ConstVector v(*this);
      return v.getNorm1();
    }
    
    double
    Vector::dot(const Vector &y) const
    {
      return ConstVector(*this).dot(ConstVector(y));
    }
    
    bool
    Vector::isFinite() const
    {
      return (ConstVector(*this)).isFinite();
    }
    
    void
    Vector::print() const
    {
      for (int i = 0; i < length(); i++)
        {
          printf("%d\t%8.4g\n", i, operator[](i));
        }
    }
    
    ConstVector::ConstVector(const Vector &v, int off, int l)
      : BaseConstVector(l, v.skip(), v.base() + v.skip()*off)
    {
      if (off < 0 || off + length() > v.length())
        {
          throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
        }
    }
    
    ConstVector::ConstVector(const ConstVector &v, int off, int l)
      : BaseConstVector(l, v.skip(), v.base() + v.skip()*off)
    {
      if (off < 0 || off + length() > v.length())
        {
          throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
        }
    }
    
    ConstVector::ConstVector(const double *d, int skip, int l)
      : BaseConstVector(l, skip, d)
    {
    }
    
    ConstVector::ConstVector(const ConstGeneralMatrix &m, int col)
      : BaseConstVector(m.numRows(), 1, &(m.get(0, col)))
    {
    }
    
    ConstVector::ConstVector(int row, const ConstGeneralMatrix &m)
      : BaseConstVector(m.numCols(), m.getLD(), &(m.get(row, 0)))
    {
    }
    
    bool
    ConstVector::operator==(const ConstVector &y) const
    {
      if (length() != y.length())
        return false;
      if (length() == 0)
        return true;
      int i = 0;
      while (i < length() && operator[](i) == y[i])
        i++;
      return i == length();
    }
    
    bool
    ConstVector::operator<(const ConstVector &y) const
    {
      int i = std::min(length(), y.length());
      int ii = 0;
      while (ii < i && operator[](ii) == y[ii])
        ii++;
      if (ii < i)
        return operator[](ii) < y[ii];
      else
        return length() < y.length();
    }
    
    double
    ConstVector::getNorm() const
    {
      double s = 0;
      for (int i = 0; i < length(); i++)
        {
          s += operator[](i)*operator[](i);
        }
      return sqrt(s);
    }
    
    double
    ConstVector::getMax() const
    {
      double r = 0;
      for (int i = 0; i < length(); i++)
        {
          if (abs(operator[](i)) > r)
            r = abs(operator[](i));
        }
      return r;
    }
    
    double
    ConstVector::getNorm1() const
    {
      double norm = 0.0;
      for (int i = 0; i < length(); i++)
        {
          norm += abs(operator[](i));
        }
      return norm;
    }
    
    double
    ConstVector::dot(const ConstVector &y) const
    {
      if (length() != y.length())
        throw SYLV_MES_EXCEPTION("Vector has different length in ConstVector::dot.");
      blas_int n = length();
      blas_int incx = skip();
      blas_int incy = y.skip();
      return ddot(&n, base(), &incx, y.base(), &incy);
    }
    
    bool
    ConstVector::isFinite() const
    {
      int i = 0;
      while (i < length() && isfinite(operator[](i)))
        i++;
      return i == length();
    }
    
    void
    ConstVector::print() const
    {
      for (int i = 0; i < length(); i++)
        {
          printf("%d\t%8.4g\n", i, operator[](i));
        }
    }
    
    ZeroPad::ZeroPad()
    {
      pad.fill(0.0);
    }