ADMB Documentation  11.1.2428
 All Classes Files Functions Variables Typedefs Friends Defines
dvector.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dvector.cpp 2280 2014-09-03 21:53:43Z 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 
00032 double avg(double x,double y)
00033 {
00034   return 0.5*(x+y);
00035 }
00036 
00042 dvector& dvector::shift(int min)
00043 {
00044   v += indexmin()-min;
00045   index_max=index_max-index_min+min;
00046   index_min=min;
00047   shape->shift(min);
00048   return *this;
00049 }
00050 
00055 dvector::~dvector()
00056 {
00057   if (shape)
00058   {
00059     if (shape->ncopies)
00060     {
00061       (shape->ncopies)--;
00062     }
00063     else
00064     {
00065       #ifdef DIAG
00066       myheapcheck(" Entering ~dvector");
00067       #endif
00068       if (v == NULL)
00069       {
00070         cerr << " Trying to delete NULL pointer in ~dvector\n";
00071         ad_exit(21);
00072       }
00073       deallocate();
00074     }
00075   }
00076 }
00077 
00083 void dvector::deallocate()
00084 {
00085   if (shape)
00086   {
00087     v = (double*)(shape->trueptr);
00088     if (v)
00089     {
00090       delete [] v;
00091       v = NULL;
00092     }
00093 
00094     delete shape;
00095     shape = NULL;
00096   }
00097 }
00098 
00102   void dvector::safe_deallocate()
00103   {
00104     if (shape)
00105     {
00106       if (shape->ncopies)
00107       {
00108         cerr << "trying to deallocate a dvector with copies" << endl;
00109         ad_exit(1);
00110       }
00111     }
00112     else
00113     {
00114       deallocate();
00115     }
00116   }
00117 
00142 dvector::dvector(const dvector& t)
00143  {
00144    #ifdef DIAG
00145     // cout << "starting out in dvector contructor\n";
00146    #endif
00147    shape=t.shape;
00148    if (shape)
00149    {
00150      (shape->ncopies)++;
00151    }
00152    else
00153    {
00154      cerr << "Making a copy of an unallocated dvector"<<endl;
00155    }
00156    v = t.v;
00157    index_min=t.index_min;
00158    index_max=t.index_max;
00159  }
00160 
00165 void dvector::shallow_copy(const dvector& t)
00166  {
00167    #ifdef DIAG
00168     // cout << "starting out in dvector contructor\n";
00169    #endif
00170    shape=t.shape;
00171    if (shape)
00172    {
00173      (shape->ncopies)++;
00174    }
00175    else
00176    {
00177      cerr << "Making a copy of an unallocated dvector"<<endl;
00178    }
00179    v = t.v;
00180    index_min=t.index_min;
00181    index_max=t.index_max;
00182  }
00183 
00190 dvector::dvector(const predvector& pdv)
00191 {
00192   int mmin = pdv.p->indexmin();
00193   int mmax = pdv.p->indexmax();
00194   if (pdv.lb < mmin || pdv.ub > mmax)
00195   {
00196     cerr << "index out of bounds in dvector subvector operator" << endl;
00197     ad_exit(1);
00198   }
00199 
00200   shape = pdv.p->shape;
00201   if (shape)
00202   {
00203     (shape->ncopies)++;
00204   }
00205   else
00206   {
00207     cerr << "Waring: Taking a subvector of an unallocated dvector.\n";
00208   }
00209 
00210   v = pdv.p->v;
00211   index_min = pdv.lb;
00212   index_max = pdv.ub;
00213 }
00214 
00222 dvector& dvector::operator=(const double x)
00223  {
00224    #ifdef DIAG
00225      myheapcheck("Entering dvector =");
00226    #endif
00227 
00228    {
00229      for (int i=indexmin();i<=indexmax();i++)
00230      {
00231        elem(i)=x;
00232      }
00233    }
00234    #ifdef DIAG
00235      myheapcheck("Leaving dvector =");
00236    #endif
00237    return (*this);
00238  }
00239 
00249 dvector& dvector::operator=(const dvector& t)
00250  {
00251    if (!(*this))
00252    {
00253      allocatec(t);
00254    }
00255 #  if defined (AD_FAST_ASSIGN)
00256    else if (!(t.shape->ncopies))
00257    {
00258      deallocate();
00259      allocatec(t);
00260    }
00261 #  endif
00262    else
00263    {
00264      if (indexmin() != t.indexmin() ||  indexmax() != t.indexmax() )
00265      {
00266        cerr << "Index bounds do not match in "
00267        "dvector& operator = (const dvector&)\n";
00268        ad_exit(1);
00269      }
00270 
00271      if (v != t.v)
00272      {
00273          int min=indexmin();
00274          int max=indexmax();
00275          memcpy(&elem(min),&t.elem(min),(max-min+1)*sizeof(double));
00276      }
00277    }
00278    return (*this);
00279  }
00280 
00291  independent_variables& independent_variables::operator=(const dvector& t)
00292  {
00293    #ifdef DIAG
00294      myheapcheck("Entering dvector =");
00295    #endif
00296 
00297    if (indexmin() != t.indexmin() ||  indexmax() != t.indexmax() )
00298    {
00299      cerr << "Index bounds do not match in "
00300   "independent_variables& independent_variables::operator=(const dvector& t)\n";
00301      ad_exit(1);
00302    }
00303      //double tmp;
00304    // disconnect dvector  pointer  from old array
00305    if (v != t.address())
00306    {
00307      for (int i=indexmin();i<=indexmax();i++)
00308      {
00309        (*this)[i]=t[i];
00310      }
00311    }
00312    #ifdef DIAG
00313      myheapcheck("Leaving dvector =");
00314    #endif
00315    return (*this);
00316  }
00317 
00324  dvector::dvector( unsigned int sz, double * x )
00325  {
00326    allocate(0,sz-1);
00327 
00328    for (unsigned int i=0; i<sz; i++)
00329    {
00330      cout << "Doing the assignment in constructor\n";
00331      v[i] = x[i];
00332    }
00333  }
00334 
00335 /*
00336  dvector::operator double* ( )
00337  {
00338    return(v);
00339  }
00340 */
00341 
00349 dvector::dvector(const dvar_vector_position& dvp, const kkludge_object& kk)
00350  {
00351    allocate(dvp.indexmin(),dvp.indexmax());
00352  }
00353 
00359  dvector::dvector( int ncl,  int nch)
00360  {
00361    if (ncl>nch)
00362      allocate();
00363    else
00364      allocate(ncl,nch);
00365  }
00366 
00371  dvector::dvector(void)
00372  {
00373    allocate();
00374  }
00375 
00383  void dvector::safe_allocate(int ncl,int nch)
00384  {
00385    if (allocated())
00386    {
00387      cerr << "trying to allocate an already allocated dvector " << endl;
00388      ad_exit(1);
00389    }
00390    else
00391    {
00392      allocate(ncl,nch);
00393    }
00394  }
00395 
00402  void dvector::allocate(int ncl,int nch)
00403  {
00404    if (ncl>nch)
00405      allocate();
00406    else
00407    {
00408      int itemp=nch-ncl;
00409      if (itemp<0)
00410      {
00411        cerr << "Error in dvector constructor max index must be >= minindex\n"
00412             << "minindex = " << ncl << " maxindex = " << nch <<endl;
00413        ad_exit(1);
00414      }
00415      if ( (v = new double [itemp+2]) ==0)
00416      {
00417        cerr << " Error trying to allocate memory for dvector\n";
00418        ad_exit(21);
00419      }
00420 #if defined(THREAD_SAFE)
00421    if ( (shape=new ts_vector_shapex(ncl,nch,v)) == NULL)
00422 #else
00423    if ( (shape=new vector_shapex(ncl,nch,v)) == NULL)
00424 #endif
00425      {
00426        cerr << "Error trying to allocate memory for dvector\n";
00427        ad_exit(1);
00428      }
00429 
00430      //int align= ((int) v)%8 ;
00431      //if (align)
00432      //{
00433      //  int diff=(8-align)%8;
00434      //  v=(double*)( ((char*)v)+diff);
00435      //}
00436 
00437      index_min=ncl;
00438      index_max=nch;
00439      v -= indexmin();
00440      #ifdef SAFE_INITIALIZE
00441        for ( int i=indexmin(); i<=indexmax(); i++)
00442        {
00443          v[i]=0.;
00444        }
00445      #endif
00446    }
00447  }
00448 
00453 void dvector::allocate(const dvector& dv)
00454 {
00455   allocate(dv.indexmin(),dv.indexmax());
00456 }
00457 
00462 void dvector::allocate(const dvar_vector& dv)
00463 {
00464   allocate(dv.indexmin(),dv.indexmax());
00465 }
00466 
00472 void dvector::allocatec(const dvector& t)
00473 {
00474   if (!(*this))
00475   {
00476     if (t.shape)
00477     {
00478       shape=t.shape;
00479       (shape->ncopies)++;
00480     }
00481     v = t.v;
00482     index_min=t.index_min;
00483     index_max=t.index_max;
00484   }
00485   else
00486   {
00487     cerr << "Trying to alocate to an already allocated dvector" << endl;
00488   }
00489 }
00490 
00495  void dvector::allocate(void)
00496  {
00497    shape=NULL;
00498    v = NULL;
00499   index_min=1;
00500   index_max=0;
00501  }
00502 
00512 double operator*(const dvector& t1, const dvector& t2)
00513   {
00514      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00515      {
00516        cerr << "Index bounds do not match in "
00517        "dvector operator * (const dvector&, const dvector&)\n";
00518        ad_exit(1);
00519      }
00520      double tmp;
00521      tmp=0;
00522    #ifdef OPT_LIB
00523      const double * pt1=&t1[t1.indexmin()];
00524      const double * pt1m=&t1[t1.indexmax()];
00525      const double * pt2=&t2[t2.indexmin()];
00526      do
00527      {
00528        tmp+=*pt1++ * *pt2++;
00529      }
00530      while (pt1<=pt1m);
00531 
00532    #else
00533      #ifndef USE_ASSEMBLER
00534        int min=t1.indexmin();
00535        int max=t1.indexmax();
00536        for (int i=min; i<=max; i++)
00537        {
00538          tmp+=t1[i]*t2[i];
00539        }
00540      #else
00541        int min=t1.indexmin();
00542        int n=t1.indexmax()-min+1;
00543        dp_dotproduct(&tmp,&(t1(min)),&(t2(min)),n);
00544      #endif
00545    #endif
00546 
00547      return(tmp);
00548   }
00549 
00559 dvector operator+(const dvector& t1, const dvector& t2)
00560   {
00561      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00562      {
00563        cerr << "Index bounds do not match in "
00564        "dvector operator+(const dvector&, const dvector&)\n";
00565        ad_exit(1);
00566      }
00567      dvector tmp(t1.indexmin(),t1.indexmax());
00568    #ifdef OPT_LIB
00569      const double * pt1=&t1[t1.indexmin()];
00570      const double * pt1m=&t1[t1.indexmax()];
00571      const double * pt2=&t2[t2.indexmin()];
00572      double * ptmp=&tmp[t1.indexmin()];
00573      do
00574      {
00575        *ptmp++=*pt1++ + *pt2++;
00576      }
00577      while (pt1<=pt1m);
00578 
00579    #else
00580 
00581      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00582      {
00583        tmp[i]=t1[i]+t2[i];
00584      }
00585    #endif
00586      return(tmp);
00587   }
00588 
00598 dvector operator-(const dvector& t1, const dvector& t2)
00599   {
00600      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00601      {
00602        cerr << "Index bounds do not match in "
00603        "dvector operator - (const dvector&, const dvector&)\n";
00604        ad_exit(1);
00605      }
00606      dvector tmp(t1.indexmin(),t1.indexmax());
00607    #ifdef OPT_LIB
00608      const double * pt1=&t1[t1.indexmin()];
00609      const double * pt1m=&t1[t1.indexmax()];
00610      const double * pt2=&t2[t2.indexmin()];
00611      double * ptmp=&tmp[t1.indexmin()];
00612      do
00613      {
00614        *ptmp++=*pt1++ - *pt2++;
00615      }
00616      while (pt1<=pt1m);
00617 
00618    #else
00619 
00620      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00621      {
00622        tmp[i]=t1[i]-t2[i];
00623      }
00624    #endif
00625      return(tmp);
00626   }
00627 
00635 dvector operator*(const double x, const dvector& t1)
00636   {
00637      dvector tmp(t1.indexmin(),t1.indexmax());
00638    #ifdef OPT_LIB
00639      const double * pt1=&t1[t1.indexmin()];
00640      const double * pt1m=&t1[t1.indexmax()];
00641      double * ptmp=&tmp[t1.indexmin()];
00642      do
00643      {
00644        *ptmp++=x * *pt1++;
00645      }
00646      while (pt1<=pt1m);
00647 
00648    #else
00649 
00650      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00651      {
00652        tmp[i]=x*t1[i];
00653      }
00654    #endif
00655      return(tmp);
00656   }
00657 /*
00658 #ifdef __TURBOC__
00659    void myheapcheck(char * msg)
00660    {
00661      if( ::heapcheck() == _HEAPCORRUPT )
00662      {
00663        cerr << msg << "Heap is corrupted.\n";
00664      }
00665      else
00666      {
00667        cerr << msg << "Heap is OK.\n";
00668      }
00669    }
00670 #else
00671 */
00677    void myheapcheck(char * msg){}
00678 
00679 /*
00680 #endif
00681 */
00682 
00689 int max(int a,int b)
00690 {
00691   if (a>b)
00692     return a;
00693   else
00694     return b;
00695 }
00696 
00702 double cube(const double m)
00703 {
00704   return m*m*m;
00705 }
00706 
00712 double fourth(const double m)
00713 {
00714   double m2=m*m;
00715   return m2*m2;
00716 }
00717