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