dvector.cpp

Go to the documentation of this file.

/*
 * $Id$
 *
 * Author: David Fournier
 * Copyright (c) 2008-2012 Regents of the University of California
 */
#ifndef OPT_LIB
  #ifndef __SUNPRO_CC
    #include <algorithm>
    using std::fill_n;
  #endif
#endif

#include "fvar.hpp"

#ifdef DIAG
  int heapcheck(void){return 0;}
#else
  int heapcheck(void);
#endif

#ifdef _MSC_VER
#include <memory.h>
#endif
#ifdef DEBUG
  #include <cassert>
  #include <climits>
#endif

double avg(double x,double y)
{
  return 0.5*(x+y);
}

dvector& dvector::shift(int min)
{
  v += indexmin()-min;
  index_max=index_max-index_min+min;
  index_min=min;
  shape->shift(min);
  return *this;
}
dvector::~dvector()
{
  if (shape)
  {
    if (shape->ncopies)
    {
      (shape->ncopies)--;
    }
    else
    {
#ifdef DEBUG
  #ifdef DIAG
      myheapcheck(" Entering ~dvector");
  #endif
      if (v == NULL)
      {
        cerr << " Trying to delete NULL pointer in ~dvector\n";
        ad_exit(21);
      }
#endif
      deallocate();
    }
  }
}
void dvector::deallocate()
{
  if (shape)
  {
    v = (double*)(shape->trueptr);
    if (v)
    {
      delete [] v;
      v = NULL;
    }

    delete shape;
    shape = NULL;
  }
}
void dvector::safe_deallocate()
{
  if (shape)
  {
    if (shape->ncopies)
    {
      cerr << "trying to deallocate a dvector with copies" << endl;
      ad_exit(1);
    }

    deallocate();
  }
}

dvector::dvector(const dvector& other)
{
  shallow_copy(other);
}
void dvector::shallow_copy(const dvector& other)
{
#if defined(DIAG)
  cout << "starting out in dvector contructor\n";
#endif
  shape = other.shape;
  if (shape)
  {
    (shape->ncopies)++;
  }
#ifdef DEBUG
  else
  {
    cerr << "Making a copy of an unallocated dvector"<<endl;
  }
#endif
  v = other.v;
  index_min = other.index_min;
  index_max = other.index_max;
}
dvector::dvector(const predvector& pdv)
{
  int mmin = pdv.p->indexmin();
  int mmax = pdv.p->indexmax();
  if (pdv.lb < mmin || pdv.ub > mmax)
  {
    cerr << "index out of bounds in dvector subvector operator" << endl;
    ad_exit(1);
  }

  shape = pdv.p->shape;
  if (shape)
  {
    (shape->ncopies)++;
  }
  else
  {
    cerr << "Waring: Taking a subvector of an unallocated dvector.\n";
  }

  v = pdv.p->v;
  index_min = pdv.lb;
  index_max = pdv.ub;
}

dvector& dvector::operator=(const double x)
{
#ifdef DIAG
  myheapcheck("Entering dvector=");
#endif

  int min = index_min;
  int max = index_max;
  double* pvi = v + min;
  for (int i = min; i <= max; ++i)
  {
    *pvi = x;
    ++pvi;
  }
#ifdef DIAG
  myheapcheck("Leaving dvector=");
#endif
   return *this;
}

dvector& dvector::operator=(const dvector& t)
{
  if (!(*this))
  {
    allocatec(t);
  }
#if defined (AD_FAST_ASSIGN)
  else if (!(t.shape->ncopies))
  {
    deallocate();
    allocatec(t);
  }
#endif
  else
  {
    int min = indexmin();

#ifndef OPT_LIB
    int max = indexmax();
    if (min != t.indexmin() || max != t.indexmax())
    {
       cerr << "Index bounds do not match in "
       "dvector& operator = (const dvector&)\n";
       ad_exit(1);
    }
#endif

    if (v != t.v)
    {
      constexpr size_t sizeofdouble = sizeof(double);
      double* pv = v + min;
      double* pt = t.get_v() + min;
      memcpy(pv, pt, size() * sizeofdouble);
    }
  }
  return *this;
}

independent_variables& independent_variables::operator=(const dvector& t)
{
#ifdef DIAG
  myheapcheck("Entering dvector =");
#endif

  int min = indexmin();
#ifndef OPT_LIB
  int max = indexmax();
  if (min != t.indexmin() || max != t.indexmax())
  {
    cerr << "Index bounds do not match in "
            "independent_variables& independent_variables::operator=(const dvector& t)\n";
    ad_exit(1);
  }
#endif
     //double tmp;
   // disconnect dvector  pointer  from old array
  if (v != t.address())
  {
    constexpr size_t sizeofdouble = sizeof(double);
    double* pv = v + min;
    double* pt = t.get_v() + min;
    memcpy(pv, pt, size() * sizeofdouble);
  }
#ifdef DIAG
  myheapcheck("Leaving dvector =");
#endif
  return *this;
}

dvector::dvector(unsigned int sz, double* x)
{
#ifdef DEBUG
  assert(sz > 0 && sz - 1 <= INT_MAX);
#endif
  allocate(0, (int)(sz - 1));
  for (unsigned int i = 0; i < sz; i++)
  {
    //cout << "Doing the assignment in constructor\n";
    v[i] = x[i];
  }
}

/*
 dvector::operator double* ( )
 {
   return(v);
 }
*/

dvector::dvector(const dvar_vector_position& dvp, const kkludge_object& kk)
 {
   allocate(dvp.indexmin(),dvp.indexmax());
 }

 dvector::dvector( int ncl,  int nch)
 {
   if (ncl>nch)
     allocate();
   else
     allocate(ncl,nch);
 }

 dvector::dvector(void)
 {
   allocate();
 }

 void dvector::safe_allocate(int ncl,int nch)
 {
   if (allocated())
   {
     cerr << "trying to allocate an already allocated dvector " << endl;
     ad_exit(1);
   }
   else
   {
     allocate(ncl,nch);
   }
 }

void dvector::allocate(int ncl, int nch)
{
  if (ncl > nch)
  {
    return allocate();
  }
  //Originally +2
  //int size = nch - ncl + 2;

  unsigned int size = static_cast<unsigned int>(nch - ncl + 1);
  v = new double[size];

  if (v == NULL)
  {
    cerr << " Error: Unable to allocate v in"
         << " dvector::allocate(int, int).\n";
    ad_exit(1);
  }

#if defined(THREAD_SAFE)
  shape = new ts_vector_shapex(ncl, nch, v);
#else
  shape = new vector_shapex(ncl, nch, v);
#endif
  if (shape == NULL)
  {
    cerr << " Error: Unable to allocate shape in"
         << " dvector::allocate(int, int).\n";
    ad_exit(1);
  }

  //reset v(index_min) to v[0]
  v -= ncl;

  index_min = ncl;
  index_max = nch;

  initialize();
}
void dvector::allocate(const dvector& other)
{
  allocate(other.indexmin(), other.indexmax());
}

void dvector::allocate(const dvar_vector& dv)
{
  allocate(dv.indexmin(),dv.indexmax());
}

void dvector::allocatec(const dvector& t)
{
  if (!(*this))
  {
    if (t.shape)
    {
      shape=t.shape;
      (shape->ncopies)++;
    }
    v = t.v;
    index_min=t.index_min;
    index_max=t.index_max;
  }
  else
  {
    cerr << "Trying to alocate to an already allocated dvector" << endl;
  }
}

void dvector::allocate()
{
  shape = NULL;
  v = NULL;
  index_min = 1;
  index_max = 0;
}

double operator*(const dvector& t1, const dvector& t2)
{
  int min = t1.indexmin();
  int max = t1.indexmax();
#ifndef OPT_LIB
  if (min != t2.indexmin() || max != t2.indexmax())
  {
    cerr << "Index bounds do not match in "
            "dvector operator * (const dvector&, const dvector&)\n";
    ad_exit(1);
  }
#endif

  double tmp{0};

#ifdef USE_ASSEMBLER
  int min=t1.indexmin();
  int n=t1.indexmax()-min+1;
  dp_dotproduct(&tmp,&(t1(min)),&(t2(min)),n);
#else
  double* pt1 = t1.get_v() + min;
  double* pt2 = t2.get_v() + min;
  for (int i = min; i <= max; ++i)
  {
    tmp += *pt1 * *pt2;
    ++pt1;
    ++pt2;
  }
#endif

  return tmp;
}

dvector operator+(const dvector& t1, const dvector& t2)
{
  int min = t1.indexmin();
  int max = t1.indexmax();

#ifndef OPT_LIB
  if (min != t2.indexmin() || max != t2.indexmax())
  {
    cerr << "Index bounds do not match in "
       "dvector operator+(const dvector&, const dvector&)\n";
    ad_exit(1);
  }
#endif
  dvector tmp(min, max);

  double* pt1 = t1.get_v() + min;
  double* pt2 = t2.get_v() + min;
  double* ptmp = tmp.get_v() + min;
  for (int i = min; i <= max; ++i)
  {
    *ptmp += *pt1 + *pt2;

    ++pt1;
    ++pt2;
    ++ptmp;
  }

  return tmp;
}

dvector operator-(const dvector& t1, const dvector& t2)
{
  int min = t1.indexmin();
  int max = t1.indexmax();

#ifndef OPT_LIB
  if (min != t2.indexmin() || max != t2.indexmax())
  {
    cerr << "Index bounds do not match in "
       "dvector operator+(const dvector&, const dvector&)\n";
    ad_exit(1);
  }
#endif
  dvector tmp(min, max);

  double* pt1 = t1.get_v() + min;
  double* pt2 = t2.get_v() + min;
  double* ptmp = tmp.get_v() + min;
  for (int i = min; i <= max; ++i)
  {
    *ptmp += *pt1 - *pt2;

    ++pt1;
    ++pt2;
    ++ptmp;
  }

  return tmp;
}

dvector operator*(const double x, const dvector& t1)
{
  int min = t1.indexmin();
  int max = t1.indexmax();

  dvector tmp(min, max);

  double* pt1 = t1.get_v() + min;
  double* ptmp = tmp.get_v() + min;
  for (int i = min; i <= max; ++i)
  {
    *ptmp = x * *pt1;

    ++pt1;
    ++ptmp;
  }

  return tmp;
}
/*
#ifdef __TURBOC__
   void myheapcheck(char * msg)
   {
     if( ::heapcheck() == _HEAPCORRUPT )
     {
       cerr << msg << "Heap is corrupted.\n";
     }
     else
     {
       cerr << msg << "Heap is OK.\n";
     }
   }
#else
*/
   void myheapcheck(char * msg){}

/*
#endif
*/

int max(int a,int b)
{
  if (a>b)
    return a;
  else
    return b;
}

double cube(const double m)
{
  return std::pow(m, 3);
}

double fourth(const double m)
{
  double m2=m*m;
  return m2*m2;
}