ADMB Documentation  11.6rc.3352
 All Classes Files Functions Variables Typedefs Friends Defines
dvector.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00012 #ifndef OPT_LIB
00013   #ifndef __SUNPRO_CC
00014     #include <algorithm>
00015     using std::fill_n;
00016   #endif
00017 #endif
00018 
00019 #include "fvar.hpp"
00020 
00021 #ifdef DOSX286
00022   int heapcheck(void){return 0;}
00023 #else
00024 
00025   int heapcheck(void);
00026 #endif
00027 
00028 #ifdef _MSC_VER
00029 #include <memory.h>
00030 #endif
00031 #ifndef OPT_LIB
00032   #include <cassert>
00033   #include <climits>
00034 #endif
00035 
00042 double avg(double x,double y)
00043 {
00044   return 0.5*(x+y);
00045 }
00046 
00052 dvector& dvector::shift(int min)
00053 {
00054   v += indexmin()-min;
00055   index_max=index_max-index_min+min;
00056   index_min=min;
00057   shape->shift(min);
00058   return *this;
00059 }
00060 
00066 dvector::~dvector()
00067 {
00068   if (shape)
00069   {
00070     if (shape->ncopies)
00071     {
00072       (shape->ncopies)--;
00073     }
00074     else
00075     {
00076 #ifdef SAFE_ALL
00077   #ifdef DIAG
00078       myheapcheck(" Entering ~dvector");
00079   #endif
00080       if (v == NULL)
00081       {
00082         cerr << " Trying to delete NULL pointer in ~dvector\n";
00083         ad_exit(21);
00084       }
00085 #endif
00086       deallocate();
00087     }
00088   }
00089 }
00093 void dvector::deallocate()
00094 {
00095   if (shape)
00096   {
00097     v = (double*)(shape->trueptr);
00098     if (v)
00099     {
00100       delete [] v;
00101       v = NULL;
00102     }
00103 
00104     delete shape;
00105     shape = NULL;
00106   }
00107 }
00111 void dvector::safe_deallocate()
00112 {
00113   if (shape)
00114   {
00115     if (shape->ncopies)
00116     {
00117       cerr << "trying to deallocate a dvector with copies" << endl;
00118       ad_exit(1);
00119     }
00120 
00121     deallocate();
00122   }
00123 }
00124 
00149 dvector::dvector(const dvector& t)
00150 {
00151 #ifdef DIAG
00152   cout << "starting out in dvector contructor\n";
00153 #endif
00154   shape=t.shape;
00155   if (shape)
00156   {
00157     (shape->ncopies)++;
00158   }
00159   else
00160   {
00161     cerr << "Making a copy of an unallocated dvector"<<endl;
00162   }
00163   v = t.v;
00164   index_min = t.index_min;
00165   index_max = t.index_max;
00166 }
00167 
00172 void dvector::shallow_copy(const dvector& t)
00173  {
00174    #ifdef DIAG
00175     // cout << "starting out in dvector contructor\n";
00176    #endif
00177    shape=t.shape;
00178    if (shape)
00179    {
00180      (shape->ncopies)++;
00181    }
00182    else
00183    {
00184      cerr << "Making a copy of an unallocated dvector"<<endl;
00185    }
00186    v = t.v;
00187    index_min=t.index_min;
00188    index_max=t.index_max;
00189  }
00190 
00197 dvector::dvector(const predvector& pdv)
00198 {
00199   int mmin = pdv.p->indexmin();
00200   int mmax = pdv.p->indexmax();
00201   if (pdv.lb < mmin || pdv.ub > mmax)
00202   {
00203     cerr << "index out of bounds in dvector subvector operator" << endl;
00204     ad_exit(1);
00205   }
00206 
00207   shape = pdv.p->shape;
00208   if (shape)
00209   {
00210     (shape->ncopies)++;
00211   }
00212   else
00213   {
00214     cerr << "Waring: Taking a subvector of an unallocated dvector.\n";
00215   }
00216 
00217   v = pdv.p->v;
00218   index_min = pdv.lb;
00219   index_max = pdv.ub;
00220 }
00221 
00229 dvector& dvector::operator=(const double x)
00230  {
00231    #ifdef DIAG
00232      myheapcheck("Entering dvector =");
00233    #endif
00234 
00235    {
00236      for (int i=indexmin();i<=indexmax();i++)
00237      {
00238        elem(i)=x;
00239      }
00240    }
00241    #ifdef DIAG
00242      myheapcheck("Leaving dvector =");
00243    #endif
00244    return (*this);
00245  }
00246 
00256 dvector& dvector::operator=(const dvector& t)
00257 {
00258   if (!(*this))
00259   {
00260     allocatec(t);
00261   }
00262 #if defined (AD_FAST_ASSIGN)
00263   else if (!(t.shape->ncopies))
00264   {
00265     deallocate();
00266     allocatec(t);
00267   }
00268 #endif
00269   else
00270   {
00271     if (indexmin() != t.indexmin() ||  indexmax() != t.indexmax() )
00272     {
00273        cerr << "Index bounds do not match in "
00274        "dvector& operator = (const dvector&)\n";
00275        ad_exit(1);
00276     }
00277 
00278     if (v != t.v)
00279     {
00280       int min=indexmin();
00281       int max=indexmax();
00282 #ifndef OPT_LIB
00283       assert(max >= min);
00284 #endif
00285       size_t size = (size_t)(max - min + 1);
00286       memcpy(&elem(min), &t.elem(min), size * sizeof(double));
00287     }
00288   }
00289   return *this;
00290 }
00291 
00302  independent_variables& independent_variables::operator=(const dvector& t)
00303  {
00304    #ifdef DIAG
00305      myheapcheck("Entering dvector =");
00306    #endif
00307 
00308    if (indexmin() != t.indexmin() ||  indexmax() != t.indexmax() )
00309    {
00310      cerr << "Index bounds do not match in "
00311   "independent_variables& independent_variables::operator=(const dvector& t)\n";
00312      ad_exit(1);
00313    }
00314      //double tmp;
00315    // disconnect dvector  pointer  from old array
00316    if (v != t.address())
00317    {
00318      for (int i=indexmin();i<=indexmax();i++)
00319      {
00320        (*this)[i]=t[i];
00321      }
00322    }
00323    #ifdef DIAG
00324      myheapcheck("Leaving dvector =");
00325    #endif
00326    return (*this);
00327  }
00328 
00335 dvector::dvector(unsigned int sz, double* x)
00336 {
00337 #ifndef OPT_LIB
00338   assert(sz > 0 && sz - 1 <= INT_MAX);
00339 #endif
00340   allocate(0, (int)(sz - 1));
00341   for (unsigned int i = 0; i < sz; i++)
00342   {
00343     //cout << "Doing the assignment in constructor\n";
00344     v[i] = x[i];
00345   }
00346 }
00347 
00348 /*
00349  dvector::operator double* ( )
00350  {
00351    return(v);
00352  }
00353 */
00354 
00362 dvector::dvector(const dvar_vector_position& dvp, const kkludge_object& kk)
00363  {
00364    allocate(dvp.indexmin(),dvp.indexmax());
00365  }
00366 
00372  dvector::dvector( int ncl,  int nch)
00373  {
00374    if (ncl>nch)
00375      allocate();
00376    else
00377      allocate(ncl,nch);
00378  }
00379 
00384  dvector::dvector(void)
00385  {
00386    allocate();
00387  }
00388 
00396  void dvector::safe_allocate(int ncl,int nch)
00397  {
00398    if (allocated())
00399    {
00400      cerr << "trying to allocate an already allocated dvector " << endl;
00401      ad_exit(1);
00402    }
00403    else
00404    {
00405      allocate(ncl,nch);
00406    }
00407  }
00408 
00415 void dvector::allocate(int ncl, int nch)
00416 {
00417   //\todo Should check if v and shape are already allocated.
00418 
00419   if (nch < ncl)
00420   {
00421     allocate();
00422   }
00423   else
00424   {
00425     //Originally +2
00426     //int size = nch - ncl + 2;
00427 
00428     int size = nch - ncl + 1;
00429 
00430     v = new double[size];
00431     if (v == NULL)
00432     {
00433       cerr << " Error trying to allocate memory for dvector\n";
00434       ad_exit(21);
00435     }
00436 #ifndef OPT_LIB
00437     assert(size >= 0);
00438     memset(v, 0, sizeof(double) * (unsigned int)size);
00439 #endif
00440 
00441 #if defined(THREAD_SAFE)
00442     shape = new ts_vector_shapex(ncl, nch, v);
00443 #else
00444     shape = new vector_shapex(ncl, nch, v);
00445 #endif
00446     if (shape == NULL)
00447     {
00448       cerr << "Error trying to allocate memory for dvector\n";
00449       ad_exit(1);
00450     }
00451 
00452     index_min = ncl;
00453     index_max = nch;
00454 
00455     //reset v(index_min) to v[0]
00456     v -= index_min;
00457   }
00458 }
00459 
00464 void dvector::allocate(const dvector& dv)
00465 {
00466   allocate(dv.indexmin(),dv.indexmax());
00467 }
00468 
00473 void dvector::allocate(const dvar_vector& dv)
00474 {
00475   allocate(dv.indexmin(),dv.indexmax());
00476 }
00477 
00483 void dvector::allocatec(const dvector& t)
00484 {
00485   if (!(*this))
00486   {
00487     if (t.shape)
00488     {
00489       shape=t.shape;
00490       (shape->ncopies)++;
00491     }
00492     v = t.v;
00493     index_min=t.index_min;
00494     index_max=t.index_max;
00495   }
00496   else
00497   {
00498     cerr << "Trying to alocate to an already allocated dvector" << endl;
00499   }
00500 }
00501 
00506 void dvector::allocate()
00507 {
00508   shape = NULL;
00509   v = NULL;
00510   index_min = 1;
00511   index_max = 0;
00512 }
00513 
00523 double operator*(const dvector& t1, const dvector& t2)
00524   {
00525      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00526      {
00527        cerr << "Index bounds do not match in "
00528        "dvector operator * (const dvector&, const dvector&)\n";
00529        ad_exit(1);
00530      }
00531      double tmp;
00532      tmp=0;
00533    #ifdef OPT_LIB
00534      const double * pt1=&t1[t1.indexmin()];
00535      const double * pt1m=&t1[t1.indexmax()];
00536      const double * pt2=&t2[t2.indexmin()];
00537      do
00538      {
00539        tmp+=*pt1++ * *pt2++;
00540      }
00541      while (pt1<=pt1m);
00542 
00543    #else
00544      #ifndef USE_ASSEMBLER
00545        int min=t1.indexmin();
00546        int max=t1.indexmax();
00547        for (int i=min; i<=max; i++)
00548        {
00549          tmp+=t1[i]*t2[i];
00550        }
00551      #else
00552        int min=t1.indexmin();
00553        int n=t1.indexmax()-min+1;
00554        dp_dotproduct(&tmp,&(t1(min)),&(t2(min)),n);
00555      #endif
00556    #endif
00557 
00558      return(tmp);
00559   }
00560 
00570 dvector operator+(const dvector& t1, const dvector& t2)
00571   {
00572      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00573      {
00574        cerr << "Index bounds do not match in "
00575        "dvector operator+(const dvector&, const dvector&)\n";
00576        ad_exit(1);
00577      }
00578      dvector tmp(t1.indexmin(),t1.indexmax());
00579    #ifdef OPT_LIB
00580      const double * pt1=&t1[t1.indexmin()];
00581      const double * pt1m=&t1[t1.indexmax()];
00582      const double * pt2=&t2[t2.indexmin()];
00583      double * ptmp=&tmp[t1.indexmin()];
00584      do
00585      {
00586        *ptmp++=*pt1++ + *pt2++;
00587      }
00588      while (pt1<=pt1m);
00589 
00590    #else
00591 
00592      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00593      {
00594        tmp[i]=t1[i]+t2[i];
00595      }
00596    #endif
00597      return(tmp);
00598   }
00599 
00609 dvector operator-(const dvector& t1, const dvector& t2)
00610   {
00611      if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
00612      {
00613        cerr << "Index bounds do not match in "
00614        "dvector operator - (const dvector&, const dvector&)\n";
00615        ad_exit(1);
00616      }
00617      dvector tmp(t1.indexmin(),t1.indexmax());
00618    #ifdef OPT_LIB
00619      const double * pt1=&t1[t1.indexmin()];
00620      const double * pt1m=&t1[t1.indexmax()];
00621      const double * pt2=&t2[t2.indexmin()];
00622      double * ptmp=&tmp[t1.indexmin()];
00623      do
00624      {
00625        *ptmp++=*pt1++ - *pt2++;
00626      }
00627      while (pt1<=pt1m);
00628 
00629    #else
00630 
00631      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00632      {
00633        tmp[i]=t1[i]-t2[i];
00634      }
00635    #endif
00636      return(tmp);
00637   }
00638 
00646 dvector operator*(const double x, const dvector& t1)
00647   {
00648      dvector tmp(t1.indexmin(),t1.indexmax());
00649    #ifdef OPT_LIB
00650      const double * pt1=&t1[t1.indexmin()];
00651      const double * pt1m=&t1[t1.indexmax()];
00652      double * ptmp=&tmp[t1.indexmin()];
00653      do
00654      {
00655        *ptmp++=x * *pt1++;
00656      }
00657      while (pt1<=pt1m);
00658 
00659    #else
00660 
00661      for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00662      {
00663        tmp[i]=x*t1[i];
00664      }
00665    #endif
00666      return(tmp);
00667   }
00668 /*
00669 #ifdef __TURBOC__
00670    void myheapcheck(char * msg)
00671    {
00672      if( ::heapcheck() == _HEAPCORRUPT )
00673      {
00674        cerr << msg << "Heap is corrupted.\n";
00675      }
00676      else
00677      {
00678        cerr << msg << "Heap is OK.\n";
00679      }
00680    }
00681 #else
00682 */
00688    void myheapcheck(char * msg){}
00689 
00690 /*
00691 #endif
00692 */
00693 
00700 int max(int a,int b)
00701 {
00702   if (a>b)
00703     return a;
00704   else
00705     return b;
00706 }
00707 
00713 double cube(const double m)
00714 {
00715   return m*m*m;
00716 }
00717 
00723 double fourth(const double m)
00724 {
00725   double m2=m*m;
00726   return m2*m2;
00727 }
00728