/*
 * $Id: dvector.cpp 494 2012-06-13 20:41:16Z johnoel $
 *
 * Author: David Fournier
 * Copyright (c) 2008-2012 Regents of the University of California 
 */
/**
 * \file
 * Description not yet available.
 */
#include "fvar.hpp"

long int farptr_tolong(void *);

#ifdef DOSX286
  int heapcheck(void){return 0;}
#else
  /// Does nothing.
  int heapcheck(void);
#endif

#ifdef __MSVC32__
#include <memory.h>
#endif

/** Average of two numbers; constant objects.
\ingroup misc
\param x A double
\param y A double
\return \f$0.5*(x+y)\f$
*/
double avg(double x,double y)
{
  return 0.5*(x+y);
}

/** Shift valid range of subscripts.
\ingroup matop
\param min Ingeger value subtracted from the subscripts
\return A %devector containing shifted vector.
*/
dvector& dvector::shift(int min)
{
  v += indexmin()-min;
  index_max=index_max-index_min+min;
  index_min=min;
  shape->shift(min);
  return *this;
}

 /** 
 Default destructor. Invoked by the compiler. Only frees allocated memory
 if all shallow copies in scope have been removed.
 */
 dvector::~dvector()
 {
   if (shape)
   {
     if (shape->ncopies)
     {
       (shape->ncopies)--;
     }
     else
     {
  
       #ifdef DIAG
         myheapcheck(" Entering ~dvector");
       #endif
  
       if ( v == NULL)
       {
         cerr << " Trying to delete NULL pointer in ~dvector\n";
         ad_exit(21);
       }
       deallocate();
     }
   }
 }

static int testflag=0;
static int ycounter=0;
  /** Deallocates memory for a dvector object.
  Called by the ~dvector(); produces an error if the double * member
  is NULL.
  */
  void dvector::deallocate()
  {
   if (testflag)
   {
     ycounter++;
     cout << " C " << ycounter << endl;
     if (ycounter==41)
       cout << "HERE" << endl;
     test_the_pointer();
   }
    if(shape)
    {
      v = (double*) (shape->trueptr);
      if (v)
      {
       if (testflag)
       {
         test_the_pointer();
       }
        delete [] v;
       if (testflag)
       {
         test_the_pointer();
       }
        v = NULL;
      }
      else
      {
        cerr << "error in dvector shape not zero but v=0" << endl;
      }
      delete shape;
      if (testflag)
      {
        test_the_pointer();
      }
    }
    else
    {
      //cerr << "Warning -- trying to deallocate an unitialized dvector"
       //    << endl;
    }
    shape=NULL;
   if (testflag)
   {
     cout << " D " << ycounter << endl;
     test_the_pointer();
   }
  }

  /**
  Deallocate memory safely. Checks if shallow copies are in scope.
  */
  void dvector::safe_deallocate()
  {
    if (shape)
    {
      if (shape->ncopies)
      {
        cerr << "trying to deallocate a dvector with copies" << endl;
        ad_exit(1);
      }
    }
    else
    {
      deallocate();
    }
  }

  /**
    Copy constructor. This constructor make a "shallow" copy. 
    Values of existing pointers are copied, and no memory is allocated.
    Invoked by the compiler in some circumstances such as function 
    call containing an argument of class %dvector passed by value. 

    The following code segments creates the \c dvector object \c x and fills it with
    uniform random numbers (using the constant 79 as seed; \ref fill_randu). The new \c dvector 
    object \c y will reference the same memory locations as \c x, and any changes in \c y
    will also appear in \c x. The new \c dvector object \c z does not refer to 
    the same memory locations as \c x, and any changes in \c z will be unique to \c z.

    \code
    dvector x(1,10);
    y.fill_randu(79L);
    dvector y = x; // shallow copy
    dvector z;
    z = x;         // "deep" copy
    \endcode
  */
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;
 }

 /**
 Explicit shallow copy.
 \param t %dvecotor to be copied
 */
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;
 }

 /**
 Creates a %dvector from an instance of class %predvector.
 Creates a shallow copy.
 \param pdv an instance of class %predvector.
 */
dvector::dvector(const predvector& pdv)
 {
   #ifdef DIAG
    // cout << "starting out in dvector contructor\n";
   #endif
   shape=pdv.p->shape;
   if (shape)
   {
     (shape->ncopies)++;
   }
   else
   {
     cerr << "Taking a subvector  of an unallocated dvector"<<endl;
   }
   v = pdv.p->v;
   index_min=pdv.lb;
   index_max=pdv.ub;
 }

 /**
  Assignment operator for double argument.
  \ingroup matop
  Assigns a the same value to all members of the vector object.
  \param x A double constant.
  \return Reference to a %dvector object.
 */
 dvector& dvector::operator = ( CGNU_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);
 }

 /**
  \ingroup matop
  Assignment operator for %dvector argument.
  Assigns the values of the argument to a target %dvector in the LHS of the
  assignment operator. The range of valid subscripts in the argument and the
  target must be identical.
  \param t A %dvector constant
  \return Reference to a %dvector object.
 */
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();
         memcpy(&elem(min),&t.elem(min),(max-min+1)*sizeof(double));
     }
   }
   return (*this);
 }

 /**
  \ingroup matop
  Assignment operator for %dvector argument.
  Assigns the values of the argument to a %independent_variables target in the LHS of the
  assignment operator. The range of valid subscripts in the argument and the
  target must be identical.
  \param t A %dvector constant
  \return Reference to a %independent_variables object.
 */
 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 dvector& operator = (const dvector&)\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);
 }

 /**
 Construct a %dvector object from a C style array of doubles.
 \param sz Number of valid memory locations allocated in the array
 The range of valid subscripts for the %dvector object will be [0,sz-1].
 \param x Pointer to double pointing to the first element in the array.
 */
 dvector::dvector( unsigned int sz, double * x )
 {
   allocate(0,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);
 }
*/

 /**
 Construct a %dvector object from a dvar_vector_position object.
 Used in writing adjoint functions callled by the autodif library,
 \param dvp Reference to a %dvar_vector_position object (usually read from stack)
 \param kk Kludge to avoid ambiguous function references
 */
dvector::dvector(const dvar_vector_position& dvp, const kkludge_object& kk)
 {
   allocate(dvp.indexmin(),dvp.indexmax());
 }

 /**
 Construct a %dvector with specified range of valid subscript.
 \param ncl Integer specifying lowest valid subscript.
 \param nch Integer specifying highest valid subscript.
 */
 dvector::dvector( int ncl,  int nch)
 {
   if (ncl>nch)
     allocate();
   else
     allocate(ncl,nch);
 }

 /**
 Construct a %dvector without allocating memory. 
 Useful in creating classes containing dvectors.
 */
 dvector::dvector(void)
 {
   allocate();
 }

 /**
 Safely allocate memory for a %dvector.
 Exits with an error message if memory for this instance has already been allocated.
 \param ncl Integer specifying lowest valid subscript.
 \param nch Integer specifying highest valid subscript.
 */
 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);
   }
 }

 /**
 Allocate memory for a %dvector.
 Exits with an error message if subscript range makes no sense.
 \param ncl Integer specifying lowest valid subscript.
 \param nch Integer specifying highest valid subscript.
 */
 void dvector::allocate(int ncl,int nch)
 {
   if (ncl>nch)
     allocate();
   else
   {
     int itemp=nch-ncl;
     if (itemp<0)
     {
       cerr << "Error in dvector constructor max index must be >= minindex\n"
            << "minindex = " << ncl << " maxindex = " << nch <<endl;
       ad_exit(1);
     }
     if ( (v = new double [itemp+2]) ==0)
     {
       cerr << " Error trying to allocate memory for dvector\n";
       ad_exit(21);
     }
#if defined(THREAD_SAFE)
   if ( (shape=new ts_vector_shapex(ncl,nch,v)) == NULL)
#else
   if ( (shape=new vector_shapex(ncl,nch,v)) == NULL)
#endif
     {
       cerr << "Error trying to allocate memory for dvector\n";
       ad_exit(1);
     }
  
     //int align= ((int) v)%8 ;
     //if (align)
     //{
     //  int diff=(8-align)%8;
     //  v=(double*)( ((char*)v)+diff);
     //}     

     index_min=ncl;
     index_max=nch;
     v -= indexmin();
     #ifdef SAFE_INITIALIZE
       for ( int i=indexmin(); i<=indexmax(); i++)
       {
         v[i]=0.;
       }
     #endif
   }
 }

 /**
 Allocate memory for a %dvector the same size as it's argument.
 \param dv Reference to a %dvector.
 */
void dvector::allocate(const dvector& dv)
{
  allocate(dv.indexmin(),dv.indexmax());
}

 /**
 Allocate memory for a %dvector the same size as it's argument.
 \param dv Reference to a dvar_vector.
 */
void dvector::allocate(const dvar_vector& dv)
{
  allocate(dv.indexmin(),dv.indexmax());
}

/**
Make shallow copy of %dvector shape.
Copies the shape of its argument.
\param t Reference to a %dvector.
*/
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;
  }
}

 /**
 Allocate %dvector without allocating memory.
 All pointers set to NULL and subscript range is invalid.
 */
 void dvector::allocate(void)
 {
   shape=NULL;
   v = NULL;
  index_min=1;
  index_max=0;
 }

  /** Compute the dot product of two vectors. The minimum and maxium
  legal subscripts of the arguments must agree; otherwize an error message
   is printed and execution terminates.
  \ingroup matop
  \param t1 A %dvector, \f$a\f$.
  \param t2 A %dvector, \f$b\f$.
  \return A double, \f$z = a\cdot b = \sum_i a_i\cdot b_i\f$  containing 
  the value of the dot product of the two arguments.
  */
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);
  }

  /**
  \ingroup matop
  Add two dvectors.
  Exits with error if bounds of the two arguments differ.
  \param t1 %dvector reference, \f$a\f$.
  \param t2 %dvector reference, \f$b\f$.
  \return A %dvector, \f$z_i = a_i + b_i\f$  containing 
  the value of the sum of the two arguments.
  */
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);
  }

  /**
  \ingroup matop
  Subtract two dvectors.
  Exits with error if bounds of the two arguments differ.
  \param t1 %dvector reference, \f$a\f$.
  \param t2 %dvector reference, \f$b\f$.
  \return A %dvector, \f$z_i = a_i - b_i\f$  containing 
  the value of the difference of the two arguments.
  */
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);
  }

  /**
  \ingroup matop
  Multiply a %dvector by a constant.
  \param x Double constant, \f$x\f$.
  \param t1 %dvector reference, \f$y\f$.
  \return A %dvector \f$z_i = x*y_i\f$.
  */
dvector operator*(CGNU_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 
*/
   /** Does nothing. This function is only defined for older Borland compilers.
   The user could provide an implementation that might be useful in certain circumstances.
   \param msg Pointer to character array.
   */
   void myheapcheck(char * msg){}

/*
#endif
*/

/** Largest of two integer; constant objects
\ingroup misc
\param a An integer
\param b An integer
\return A integer containing  \f$ z = \max(a,b)\f$
*/
int max(int a,int b)
{
  if (a>b)
    return a;
  else
    return b;
}

/** Cube of a number; constant objects.
\ingroup misc
\param m Number to be cubed
\return \f$m^3\f$
*/
double cube( CGNU_DOUBLE m)
{
  return m*m*m;
}

/** Fourth power of a number; constant objects.
\ingroup misc
\param m Number to be taken to the fourth power.
\return \f$m^4\f$
*/
double fourth( CGNU_DOUBLE m)
{
  double m2=m*m;
  return m2*m2;
}

