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 DOSX286
  int heapcheck(void){return 0;}
#else

  int heapcheck(void);
#endif

#ifdef _MSC_VER
#include <memory.h>
#endif
#ifndef OPT_LIB
  #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 SAFE_ALL
  #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& t)
{
#ifdef DIAG
  cout << "starting out in dvector contructor\n";
#endif
  shape=t.shape;
  if (shape)
  {
    (shape->ncopies)++;
  }
  else
  {
    cerr << "Making a copy of an unallocated dvector"<<endl;
  }
  v = t.v;
  index_min = t.index_min;
  index_max = t.index_max;
}

void dvector::shallow_copy(const dvector& t)
 {
   #ifdef DIAG
    // cout << "starting out in dvector contructor\n";
   #endif
   shape=t.shape;
   if (shape)
   {
     (shape->ncopies)++;
   }
   else
   {
     cerr << "Making a copy of an unallocated dvector"<<endl;
   }
   v = t.v;
   index_min=t.index_min;
   index_max=t.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

   {
     for (int i=indexmin();i<=indexmax();i++)
     {
       elem(i)=x;
     }
   }
   #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
  {
    if (indexmin() != t.indexmin() ||  indexmax() != t.indexmax() )
    {
       cerr << "Index bounds do not match in "
       "dvector& operator = (const dvector&)\n";
       ad_exit(1);
    }

    if (v != t.v)
    {
      int min=indexmin();
      int max=indexmax();
#ifndef OPT_LIB
      assert(max >= min);
#endif
      size_t size = (size_t)(max - min + 1);
      memcpy(&elem(min), &t.elem(min), size * sizeof(double));
    }
  }
  return *this;
}

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

   if (indexmin() != t.indexmin() ||  indexmax() != t.indexmax() )
   {
     cerr << "Index bounds do not match in "
  "independent_variables& independent_variables::operator=(const dvector& t)\n";
     ad_exit(1);
   }
     //double tmp;
   // disconnect dvector  pointer  from old array
   if (v != t.address())
   {
     for (int i=indexmin();i<=indexmax();i++)
     {
       (*this)[i]=t[i];
     }
   }
   #ifdef DIAG
     myheapcheck("Leaving dvector =");
   #endif
   return (*this);
 }

dvector::dvector(unsigned int sz, double* x)
{
#ifndef OPT_LIB
  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)
{
  //\todo Should check if v and shape are already allocated.

  if (nch < ncl)
  {
    allocate();
  }
  else
  {
    //Originally +2
    //int size = nch - ncl + 2;

    int size = nch - ncl + 1;

    v = new double[size];
    if (v == NULL)
    {
      cerr << " Error trying to allocate memory for dvector\n";
      ad_exit(21);
    }
#ifndef OPT_LIB
    assert(size >= 0);
    memset(v, 0, sizeof(double) * (unsigned int)size);
#endif

#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 trying to allocate memory for dvector\n";
      ad_exit(1);
    }

    index_min = ncl;
    index_max = nch;

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

void dvector::allocate(const dvector& dv)
{
  allocate(dv.indexmin(),dv.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)
  {
     if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
     {
       cerr << "Index bounds do not match in "
       "dvector operator * (const dvector&, const dvector&)\n";
       ad_exit(1);
     }
     double tmp;
     tmp=0;
   #ifdef OPT_LIB
     const double * pt1=&t1[t1.indexmin()];
     const double * pt1m=&t1[t1.indexmax()];
     const double * pt2=&t2[t2.indexmin()];
     do
     {
       tmp+=*pt1++ * *pt2++;
     }
     while (pt1<=pt1m);

   #else
     #ifndef USE_ASSEMBLER
       int min=t1.indexmin();
       int max=t1.indexmax();
       for (int i=min; i<=max; i++)
       {
         tmp+=t1[i]*t2[i];
       }
     #else
       int min=t1.indexmin();
       int n=t1.indexmax()-min+1;
       dp_dotproduct(&tmp,&(t1(min)),&(t2(min)),n);
     #endif
   #endif

     return(tmp);
  }

dvector operator+(const dvector& t1, const dvector& t2)
  {
     if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
     {
       cerr << "Index bounds do not match in "
       "dvector operator+(const dvector&, const dvector&)\n";
       ad_exit(1);
     }
     dvector tmp(t1.indexmin(),t1.indexmax());
   #ifdef OPT_LIB
     const double * pt1=&t1[t1.indexmin()];
     const double * pt1m=&t1[t1.indexmax()];
     const double * pt2=&t2[t2.indexmin()];
     double * ptmp=&tmp[t1.indexmin()];
     do
     {
       *ptmp++=*pt1++ + *pt2++;
     }
     while (pt1<=pt1m);

   #else

     for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
     {
       tmp[i]=t1[i]+t2[i];
     }
   #endif
     return(tmp);
  }

dvector operator-(const dvector& t1, const dvector& t2)
  {
     if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
     {
       cerr << "Index bounds do not match in "
       "dvector operator - (const dvector&, const dvector&)\n";
       ad_exit(1);
     }
     dvector tmp(t1.indexmin(),t1.indexmax());
   #ifdef OPT_LIB
     const double * pt1=&t1[t1.indexmin()];
     const double * pt1m=&t1[t1.indexmax()];
     const double * pt2=&t2[t2.indexmin()];
     double * ptmp=&tmp[t1.indexmin()];
     do
     {
       *ptmp++=*pt1++ - *pt2++;
     }
     while (pt1<=pt1m);

   #else

     for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
     {
       tmp[i]=t1[i]-t2[i];
     }
   #endif
     return(tmp);
  }

dvector operator*(const double x, const dvector& t1)
  {
     dvector tmp(t1.indexmin(),t1.indexmax());
   #ifdef OPT_LIB
     const double * pt1=&t1[t1.indexmin()];
     const double * pt1m=&t1[t1.indexmax()];
     double * ptmp=&tmp[t1.indexmin()];
     do
     {
       *ptmp++=x * *pt1++;
     }
     while (pt1<=pt1m);

   #else

     for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
     {
       tmp[i]=x*t1[i];
     }
   #endif
     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 m*m*m;
}

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