ADMB Documentation  11.2.2826
 All Classes Files Functions Variables Typedefs Friends Defines
dvector.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dvector.cpp 2633 2014-11-12 22:04:36Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00013 long int farptr_tolong(void *);
00014 
00015 #ifdef DOSX286
00016   int heapcheck(void){return 0;}
00017 #else
00018 
00019   int heapcheck(void);
00020 #endif
00021 
00022 #ifdef _MSC_VER
00023 #include <memory.h>
00024 #endif
00025 #ifndef OPT_LIB
00026   #include <cassert>
00027   #include <climits>
00028 #endif
00029 
00036 double avg(double x,double y)
00037 {
00038   return 0.5*(x+y);
00039 }
00040 
00046 dvector& dvector::shift(int min)
00047 {
00048   v += indexmin()-min;
00049   index_max=index_max-index_min+min;
00050   index_min=min;
00051   shape->shift(min);
00052   return *this;
00053 }
00054 
00059 dvector::~dvector()
00060 {
00061   if (shape)
00062   {
00063     if (shape->ncopies)
00064     {
00065       (shape->ncopies)--;
00066     }
00067     else
00068     {
00069       #ifdef DIAG
00070       myheapcheck(" Entering ~dvector");
00071       #endif
00072       if (v == NULL)
00073       {
00074         cerr << " Trying to delete NULL pointer in ~dvector\n";
00075         ad_exit(21);
00076       }
00077       deallocate();
00078     }
00079   }
00080 }
00081 
00087 void dvector::deallocate()
00088 {
00089   if (shape)
00090   {
00091     v = (double*)(shape->trueptr);
00092     if (v)
00093     {
00094       delete [] v;
00095       v = NULL;
00096     }
00097 
00098     delete shape;
00099     shape = NULL;
00100   }
00101 }
00102 
00106   void dvector::safe_deallocate()
00107   {
00108     if (shape)
00109     {
00110       if (shape->ncopies)
00111       {
00112         cerr << "trying to deallocate a dvector with copies" << endl;
00113         ad_exit(1);
00114       }
00115     }
00116     else
00117     {
00118       deallocate();
00119     }
00120   }
00121 
00146 dvector::dvector(const dvector& t)
00147  {
00148    #ifdef DIAG
00149     // cout << "starting out in dvector contructor\n";
00150    #endif
00151    shape=t.shape;
00152    if (shape)
00153    {
00154      (shape->ncopies)++;
00155    }
00156    else
00157    {
00158      cerr << "Making a copy of an unallocated dvector"<<endl;
00159    }
00160    v = t.v;
00161    index_min=t.index_min;
00162    index_max=t.index_max;
00163  }
00164 
00169 void dvector::shallow_copy(const dvector& t)
00170  {
00171    #ifdef DIAG
00172     // cout << "starting out in dvector contructor\n";
00173    #endif
00174    shape=t.shape;
00175    if (shape)
00176    {
00177      (shape->ncopies)++;
00178    }
00179    else
00180    {
00181      cerr << "Making a copy of an unallocated dvector"<<endl;
00182    }
00183    v = t.v;
00184    index_min=t.index_min;
00185    index_max=t.index_max;
00186  }
00187 
00194 dvector::dvector(const predvector& pdv)
00195 {
00196   int mmin = pdv.p->indexmin();
00197   int mmax = pdv.p->indexmax();
00198   if (pdv.lb < mmin || pdv.ub > mmax)
00199   {
00200     cerr << "index out of bounds in dvector subvector operator" << endl;
00201     ad_exit(1);
00202   }
00203 
00204   shape = pdv.p->shape;
00205   if (shape)
00206   {
00207     (shape->ncopies)++;
00208   }
00209   else
00210   {
00211     cerr << "Waring: Taking a subvector of an unallocated dvector.\n";
00212   }
00213 
00214   v = pdv.p->v;
00215   index_min = pdv.lb;
00216   index_max = pdv.ub;
00217 }
00218 
00226 dvector& dvector::operator=(const double x)
00227  {
00228    #ifdef DIAG
00229      myheapcheck("Entering dvector =");
00230    #endif
00231 
00232    {
00233      for (int i=indexmin();i<=indexmax();i++)
00234      {
00235        elem(i)=x;
00236      }
00237    }
00238    #ifdef DIAG
00239      myheapcheck("Leaving dvector =");
00240    #endif
00241    return (*this);
00242  }
00243 
00253 dvector& dvector::operator=(const dvector& t)
00254  {
00255    if (!(*this))
00256    {
00257      allocatec(t);
00258    }
00259 #  if defined (AD_FAST_ASSIGN)
00260    else if (!(t.shape->ncopies))
00261    {
00262      deallocate();
00263      allocatec(t);
00264    }
00265 #  endif
00266    else
00267    {
00268      if (indexmin() != t.indexmin() ||  indexmax() != t.indexmax() )
00269      {
00270        cerr << "Index bounds do not match in "
00271        "dvector& operator = (const dvector&)\n";
00272        ad_exit(1);
00273      }
00274 
00275      if (v != t.v)
00276      {
00277          int min=indexmin();
00278          int max=indexmax();
00279 #ifndef OPT_LIB
00280          assert(max >= min);
00281 #endif
00282          size_t size = (size_t)(max - min + 1);
00283          memcpy(&elem(min), &t.elem(min), size * sizeof(double));
00284      }
00285    }
00286    return (*this);
00287  }
00288 
00299  independent_variables& independent_variables::operator=(const dvector& t)
00300  {
00301    #ifdef DIAG
00302      myheapcheck("Entering dvector =");
00303    #endif
00304 
00305    if (indexmin() != t.indexmin() ||  indexmax() != t.indexmax() )
00306    {
00307      cerr << "Index bounds do not match in "
00308   "independent_variables& independent_variables::operator=(const dvector& t)\n";
00309      ad_exit(1);
00310    }
00311      //double tmp;
00312    // disconnect dvector  pointer  from old array
00313    if (v != t.address())
00314    {
00315      for (int i=indexmin();i<=indexmax();i++)
00316      {
00317        (*this)[i]=t[i];
00318      }
00319    }
00320    #ifdef DIAG
00321      myheapcheck("Leaving dvector =");
00322    #endif
00323    return (*this);
00324  }
00325 
00332 dvector::dvector(unsigned int sz, double* x)
00333 {
00334 #ifndef OPT_LIB
00335   assert(sz > 0 && sz - 1 <= INT_MAX);
00336 #endif
00337   allocate(0, (int)(sz - 1));
00338   for (unsigned int i = 0; i < sz; i++)
00339   {
00340     //cout << "Doing the assignment in constructor\n";
00341     v[i] = x[i];
00342   }
00343 }
00344 
00345 /*
00346  dvector::operator double* ( )
00347  {
00348    return(v);
00349  }
00350 */
00351 
00359 dvector::dvector(const dvar_vector_position& dvp, const kkludge_object& kk)
00360  {
00361    allocate(dvp.indexmin(),dvp.indexmax());
00362  }
00363 
00369  dvector::dvector( int ncl,  int nch)
00370  {
00371    if (ncl>nch)
00372      allocate();
00373    else
00374      allocate(ncl,nch);
00375  }
00376 
00381  dvector::dvector(void)
00382  {
00383    allocate();
00384  }
00385 
00393  void dvector::safe_allocate(int ncl,int nch)
00394  {
00395    if (allocated())
00396    {
00397      cerr << "trying to allocate an already allocated dvector " << endl;
00398      ad_exit(1);
00399    }
00400    else
00401    {
00402      allocate(ncl,nch);
00403    }
00404  }
00405 
00412  void dvector::allocate(int ncl,int nch)
00413  {
00414    if (ncl>nch)
00415      allocate();
00416    else
00417    {
00418      int itemp=nch-ncl;
00419      if (itemp<0)
00420      {
00421        cerr << "Error in dvector constructor max index must be >= minindex\n"
00422             << "minindex = " << ncl << " maxindex = " << nch <<endl;
00423        ad_exit(1);
00424      }
00425      if ( (v = new double [itemp+2]) ==0)
00426      {
00427        cerr << " Error trying to allocate memory for dvector\n";
00428        ad_exit(21);
00429      }
00430 #if defined(THREAD_SAFE)
00431    if ( (shape=new ts_vector_shapex(ncl,nch,v)) == NULL)
00432 #else
00433    if ( (shape=new vector_shapex(ncl,nch,v)) == NULL)
00434 #endif
00435      {
00436        cerr << "Error trying to allocate memory for dvector\n";
00437        ad_exit(1);
00438      }
00439 
00440      //int align= ((int) v)%8 ;
00441      //if (align)
00442      //{
00443      //  int diff=(8-align)%8;
00444      //  v=(double*)( ((char*)v)+diff);
00445      //}
00446 
00447      index_min=ncl;
00448      index_max=nch;
00449      v -= indexmin();
00450      #ifdef SAFE_INITIALIZE
00451        for ( int i=indexmin(); i<=indexmax(); i++)
00452        {
00453          v[i]=0.;
00454        }
00455      #endif
00456    }
00457  }
00458 
00463 void dvector::allocate(const dvector& dv)
00464 {
00465   allocate(dv.indexmin(),dv.indexmax());
00466 }
00467 
00472 void dvector::allocate(const dvar_vector& dv)
00473 {
00474   allocate(dv.indexmin(),dv.indexmax());
00475 }
00476 
00482 void dvector::allocatec(const dvector& t)
00483 {
00484   if (!(*this))
00485   {
00486     if (t.shape)
00487     {
00488       shape=t.shape;
00489       (shape->ncopies)++;
00490     }
00491     v = t.v;
00492     index_min=t.index_min;
00493     index_max=t.index_max;
00494   }
00495   else
00496   {
00497     cerr << "Trying to alocate to an already allocated dvector" << endl;
00498   }
00499 }
00500 
00505  void dvector::allocate(void)
00506  {
00507    shape=NULL;
00508    v = NULL;
00509   index_min=1;
00510   index_max=0;
00511  }
00512 
00522 double operator*(const dvector& t1, const dvector& t2)
00523   {
00524      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00525      {
00526        cerr << "Index bounds do not match in "
00527        "dvector operator * (const dvector&, const dvector&)\n";
00528        ad_exit(1);
00529      }
00530      double tmp;
00531      tmp=0;
00532    #ifdef OPT_LIB
00533      const double * pt1=&t1[t1.indexmin()];
00534      const double * pt1m=&t1[t1.indexmax()];
00535      const double * pt2=&t2[t2.indexmin()];
00536      do
00537      {
00538        tmp+=*pt1++ * *pt2++;
00539      }
00540      while (pt1<=pt1m);
00541 
00542    #else
00543      #ifndef USE_ASSEMBLER
00544        int min=t1.indexmin();
00545        int max=t1.indexmax();
00546        for (int i=min; i<=max; i++)
00547        {
00548          tmp+=t1[i]*t2[i];
00549        }
00550      #else
00551        int min=t1.indexmin();
00552        int n=t1.indexmax()-min+1;
00553        dp_dotproduct(&tmp,&(t1(min)),&(t2(min)),n);
00554      #endif
00555    #endif
00556 
00557      return(tmp);
00558   }
00559 
00569 dvector operator+(const dvector& t1, const dvector& t2)
00570   {
00571      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00572      {
00573        cerr << "Index bounds do not match in "
00574        "dvector operator+(const dvector&, const dvector&)\n";
00575        ad_exit(1);
00576      }
00577      dvector tmp(t1.indexmin(),t1.indexmax());
00578    #ifdef OPT_LIB
00579      const double * pt1=&t1[t1.indexmin()];
00580      const double * pt1m=&t1[t1.indexmax()];
00581      const double * pt2=&t2[t2.indexmin()];
00582      double * ptmp=&tmp[t1.indexmin()];
00583      do
00584      {
00585        *ptmp++=*pt1++ + *pt2++;
00586      }
00587      while (pt1<=pt1m);
00588 
00589    #else
00590 
00591      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00592      {
00593        tmp[i]=t1[i]+t2[i];
00594      }
00595    #endif
00596      return(tmp);
00597   }
00598 
00608 dvector operator-(const dvector& t1, const dvector& t2)
00609   {
00610      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00611      {
00612        cerr << "Index bounds do not match in "
00613        "dvector operator - (const dvector&, const dvector&)\n";
00614        ad_exit(1);
00615      }
00616      dvector tmp(t1.indexmin(),t1.indexmax());
00617    #ifdef OPT_LIB
00618      const double * pt1=&t1[t1.indexmin()];
00619      const double * pt1m=&t1[t1.indexmax()];
00620      const double * pt2=&t2[t2.indexmin()];
00621      double * ptmp=&tmp[t1.indexmin()];
00622      do
00623      {
00624        *ptmp++=*pt1++ - *pt2++;
00625      }
00626      while (pt1<=pt1m);
00627 
00628    #else
00629 
00630      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00631      {
00632        tmp[i]=t1[i]-t2[i];
00633      }
00634    #endif
00635      return(tmp);
00636   }
00637 
00645 dvector operator*(const double x, const dvector& t1)
00646   {
00647      dvector tmp(t1.indexmin(),t1.indexmax());
00648    #ifdef OPT_LIB
00649      const double * pt1=&t1[t1.indexmin()];
00650      const double * pt1m=&t1[t1.indexmax()];
00651      double * ptmp=&tmp[t1.indexmin()];
00652      do
00653      {
00654        *ptmp++=x * *pt1++;
00655      }
00656      while (pt1<=pt1m);
00657 
00658    #else
00659 
00660      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00661      {
00662        tmp[i]=x*t1[i];
00663      }
00664    #endif
00665      return(tmp);
00666   }
00667 /*
00668 #ifdef __TURBOC__
00669    void myheapcheck(char * msg)
00670    {
00671      if( ::heapcheck() == _HEAPCORRUPT )
00672      {
00673        cerr << msg << "Heap is corrupted.\n";
00674      }
00675      else
00676      {
00677        cerr << msg << "Heap is OK.\n";
00678      }
00679    }
00680 #else
00681 */
00687    void myheapcheck(char * msg){}
00688 
00689 /*
00690 #endif
00691 */
00692 
00699 int max(int a,int b)
00700 {
00701   if (a>b)
00702     return a;
00703   else
00704     return b;
00705 }
00706 
00712 double cube(const double m)
00713 {
00714   return m*m*m;
00715 }
00716 
00722 double fourth(const double m)
00723 {
00724   double m2=m*m;
00725   return m2*m2;
00726 }
00727