d4arr.cpp

Go to the documentation of this file.

#include "fvar.hpp"
#include "admb_messages.h"

d4_array::d4_array(int nrl, int nrh)
{
  allocate(nrl, nrh);
}
four_array_shape::four_array_shape(int hsl, int hsu):
  ncopies(0),
  hslice_min(hsl),
  hslice_max(hsu)
{
  //,int sl,int su,int rl,
  // int ru,int cl,int cu)
  //slice_min=sl;
  //slice_max=su;
  //row_min=rl;
  //row_max=ru;
  //col_min=cl;
  //col_max=cu;
}
double sum(const d4_array& darray)
{
  double total = 0.0;

  int min = darray.indexmin();
  int max = darray.indexmax();
  const d3_array* pdarrayi = &darray(min);
  for (int i = min; i <= max; ++i)
  {
    total += sum(*pdarrayi);
    ++pdarrayi;
  }
  return total;
}

 d4_array d4_array::sub(int nrl,int nrh)
 {
   if (allocated(*this))
   {
     d4_array tmp(nrl,nrh);
     for (int i=nrl; i<=nrh; i++)
     {
       tmp[i].shallow_copy((*this)(i));
     }
     return tmp;
   }
   else
   {
     return *this;
   }
 }
d4_array::d4_array(const d4_array& other)
{
  shallow_copy(other);
}
void d4_array::shallow_copy(const d4_array& other)
{
  if (other.shape)
  {
    shape = other.shape;
    ++(shape->ncopies);
    t = other.t;
  }
  else
  {
#ifdef DEBUG
    cerr << "Warning -- Unable to shallow copy an unallocated d4_array.\n";
#endif
    allocate();
  }
}
void d4_array::deallocate()
{
  if (shape)
  {
    if (shape->ncopies > 0)
    {
      --(shape->ncopies);
    }
    else
    {
      t += indexmin();
      delete [] t;
      delete shape;
    }
    allocate();
  }
#if defined(DIAG)
  else
  {
    cerr << "Warning -- Unable to deallocate an unallocated d4_array.\n";
    //ad_exit(1);
  }
#endif
}
d4_array::~d4_array()
{
  deallocate();
}
d4_array& d4_array::operator=(const d4_array& other)
{
  int min = hslicemin();
  int max = hslicemax();
  if (min != other.hslicemin() || max != other.hslicemax())
  {
    cerr << "Incompatible bounds in"
         << " d4_array& d4_array::operator=(const d4_array&).\n";
    ad_exit(1);
  }
  d3_array* pti = t + min;
  const d3_array* potheri = &other(min);
  for (int i = min; i <= max; ++i)
  {
    *pti = *potheri;

    ++pti;
    ++potheri;
  }
  return *this;
}
void d4_array::allocate(const d4_array& other)
{
  int min = other.hslicemin();
  int max = other.hslicemax();
  allocate(min, max);
  d3_array* pti = t + min;
  const d3_array* potheri = &other(min);
  for (int i = min; i <= max; ++i)
  {
    pti->allocate(*potheri);

    ++pti;
    ++potheri;
  }
}

  #ifndef OPT_LIB

    d3_array& d4_array::operator () (int i)
    {
      if (i < hslicemin() || i > hslicemax())
      {
        ADMB_ARRAY_BOUNDS_ERROR("hslice index out of bounds",
        "d3_array& d4_array::operator() (int i)", hslicemin(), hslicemax(), i);
      }
      return t[i];
    }

    d3_array& d4_array::operator [] (int i)
    {
      if (i < hslicemin() || i > hslicemax())
      {
        ADMB_ARRAY_BOUNDS_ERROR("hslice index out of bounds",
        "d3_array& d4_array::operator[] (int i)",
        hslicemin(), hslicemax(), i);
      }
      return t[i];
    }

    dmatrix& d4_array::operator ( ) (int i ,int j)
    {
      if (i < hslicemin() || i > hslicemax())
      {
        ADMB_ARRAY_BOUNDS_ERROR("hslice index out of bounds",
        "dmatrix& d4_array::operator() (int i, int j)",
        hslicemin(), hslicemax(), i);
      }
      return ((*this)(i))(j);
    }

    dvector& d4_array::operator ( ) (int i,int j,int k)
    {
      if (i < hslicemin() || i > hslicemax())
      {
        ADMB_ARRAY_BOUNDS_ERROR("hslice index out of bounds",
        "dvector& d4_array::operator() (int i, int j, int k)",
        hslicemin(), hslicemax(), i);
      }
      return (((*this)(i,j))(k));
    }

    double& d4_array::operator ( ) (int i,int j,int k,int l)
    {
      if (i < hslicemin() || i > hslicemax())
      {
        ADMB_ARRAY_BOUNDS_ERROR("hslice index out of bounds",
        "double& d4_array::operator() (int i, int j, int k, int l)",
        hslicemin(), hslicemax(), i);
      }
      return ( ((*this)(i,j,k))(l));
    }

const d3_array& d4_array::operator()(int i) const
    {
        if (i<hslicemin()||i>hslicemax())
        { cerr << "Error hslice index out of bounds in\n"
            "d3_array& d4_array::operator ( )" << endl;
          ad_exit(1);
        }
      return t[i];
    }

const d3_array& d4_array::operator[](int i) const
    {
        if (i<hslicemin()||i>hslicemax())
        { cerr << "Error hslice index out of bounds in\n"
            "d3_array& d4_array::operator ( )" << endl;
          ad_exit(1);
        }
      return t[i];
    }

const dmatrix& d4_array::operator()(int i, int j) const
    {
        if (i<hslicemin()||i>hslicemax())
        { cerr << "Error hslice index out of bounds in\n"
            "dmatrix& d4_array::operator ( )" << endl;
          ad_exit(1);
        }
      return ((*this)(i))(j);
    }

const dvector& d4_array::operator()(int i, int j, int k) const
    {
        if (i<hslicemin()||i>hslicemax())
        { cerr << "Error hslice index out of bounds in\n"
          "dvector& d4_array::operator ( )" << endl;
          ad_exit(1);
        }
      return (((*this)(i,j))(k));
    }

const double& d4_array::operator()(int i, int j, int k, int l) const
    {
        if (i<hslicemin()||i>hslicemax())
        { cerr << "Error hslice index out of bounds in\n"
            "double& d4_array::operator ( )"  << endl;
          ad_exit(1);
        }
      return ( ((*this)(i,j,k))(l));
    }

#endif

void d4_array::allocate(
  int hsl, int hsu,
  int sl, int sh,
  int nrl, int nrh,
  int ncl, int nch)
{
  if ((shape = new four_array_shape(hsl, hsu)) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ( (t = new d3_array[hslicesize()]) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  t -= hslicemin();
  d3_array* pti = t + hsl;
  for (int i = hsl; i <= hsu; ++i)
  {
    pti->allocate(sl, sh, nrl, nrh, ncl, nch);
    ++pti;
  }
}
void d4_array::allocate(
  int hsl, int hsu,
  int sl, int sh,
  int nrl, int nrh,
  const ivector& ncl, const ivector& nch)
{
  if ((shape = new four_array_shape(hsl, hsu)) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ((t = new d3_array[hslicesize()]) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  t -= hslicemin();
  d3_array* pti = t + hsl;
  for (int i = hsl; i <= hsu; ++i)
  {
    pti->allocate(sl, sh, nrl, nrh, ncl, nch);
    ++pti;
  }
}
void d4_array::allocate(
  int hsl, int hsu,
  int sl, int sh,
  const ivector& nrl, const ivector& nrh,
  const ivector& ncl, const ivector& nch)
{
  if ( (shape = new four_array_shape(hsl, hsu)) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ( (t = new d3_array[hslicesize()]) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  t -= hslicemin();
  d3_array* pti = t + hsl;
  const int* pnrli = nrl.get_v() + hsl;
  const int* pnrhi = nrh.get_v() + hsl;
  const int* pncli = ncl.get_v() + hsl;
  const int* pnchi = nch.get_v() + hsl;
  for (int i = hsl; i <= hsu; ++i)
  {
    pti->allocate(sl, sh, *pnrli, *pnrhi, *pncli, *pnchi);
    ++pti;
    ++pnrli;
    ++pnrhi;
    ++pncli;
    ++pnchi;
  }
}
d4_array::d4_array(int hsl, int hsu, int sl, const ivector& sh,
  int nrl, const imatrix& nrh, int ncl, int nch)
{
  allocate(hsl,hsu,sl,sh,nrl,nrh,ncl,nch);
}

void d4_array::allocate(
  int hsl, int hsu,
  int sl, const ivector& sh,
  int nrl, const imatrix& nrh,
  int ncl, int nch)
{
  if ((shape = new four_array_shape(hsl, hsu)) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ((t = new d3_array[hslicesize()]) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  t -= hslicemin();
  d3_array* pti = t + hsl;
  const int* pshi = sh.get_v() + hsl;
  const ivector* pnrhi = &nrh(hsl);
  for (int i = hsl; i <= hsu; ++i)
  {
    pti->allocate(sl, *pshi, nrl, *pnrhi, ncl, nch);
    ++pti;
    ++pshi;
    ++pnrhi;
  }
}

 d4_array::d4_array(int hsl,int hsu,int sl,int sh,int nrl,
   int nrh,int ncl,int nch)
 {
   allocate(hsl,hsu,sl,sh,nrl,nrh,ncl,nch);
 }

d4_array::d4_array(int hsl,int hsu, int sl,int sh,ivector nrl,ivector nrh,
  ivector ncl,ivector nch)
{
  allocate(hsl,hsu,sl,sh,nrl,nrh,ncl,nch);
}

void d4_array::initialize()
{
  if (!(!(*this)))  // only initialize allocated objects
  {
    int min = hslicemin();
    int max = hslicemax();
    d3_array* pti = t + min;
    for (int i = min; i <= max; ++i)
    {
      pti->initialize();
      ++pti;
    }
  }
}

d4_array::d4_array(int hsl, int hsu, int sl, const ivector& sh, int nrl,
  const imatrix& nrh, int ncl, const i3_array& nch)
{
  allocate(hsl,hsu,sl,sh,nrl,nrh,ncl,nch);
}

void d4_array::allocate(
  int hsl, int hsu,
  int sl, const ivector& sh,
  int nrl, const imatrix& nrh,
  int ncl, const i3_array& nch)
{
  if ((shape = new four_array_shape(hsl, hsu)) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ((t = new d3_array[hslicesize()]) == 0)
  {
    cerr << " Error: d4_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  t -= hslicemin();
  d3_array* pti = t + hsl;
  const int* pshi = sh.get_v() + hsl;
  const ivector* pnrhi = &nrh(hsl);
  const imatrix* pnchi = &nch(hsl);
  for (int i = hsl; i <= hsu; ++i)
  {
    pti->allocate(sl, *pshi, nrl, *pnrhi, ncl, *pnchi);

    ++pti;
    ++pshi;
    ++pnchi;
  }
}

d4_array::d4_array(int hsl,int hsu,const index_type& sl,
    const index_type& sh,const index_type& nrl,const index_type& nrh,
    const index_type& ncl,const index_type& nch)
{
  allocate(hsl,hsu,sl,sh,nrl,nrh,ncl,nch);
}

void d4_array::allocate(
  ad_integer hsl, ad_integer hsu,
  const index_type& sl, const index_type& sh,
  const index_type& nrl, const index_type& nrh,
  const index_type& ncl, const index_type& nch)
{
  unsigned int ss =
    static_cast<unsigned int>(hsu < hsl ? 0 : hsu - hsl + 1);
  if (ss > 0)
  {
    if ( (t = new d3_array[ss]) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
      ad_exit(1);
    }
    if ( (shape=new four_array_shape(hsl,hsu)) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
      ad_exit(1);
    }
    t -= int(hsl);
    d3_array* pti = t + hsl;
    for (int i=hsl; i<=hsu; i++)
    {
      const index_type& rsl=sl[i];
      const index_type& rsh=sh[i];
      const index_type& rnrl=nrl[i];
      const index_type& rnrh=nrh[i];
      const index_type& rncl=ncl[i];
      const index_type& rnch=nch[i];
      const ad_integer& aa=ad_integer(rsl);
      const ad_integer& bb=ad_integer(rsh);
      pti->allocate(aa,bb,rnrl,rnrh,rncl,rnch);
      ++pti;
      //(*this)(i).allocate(sl[i],sh[i],nrl[i],nrh[i],ncl[i],nch[i]);
    }
  }
  else
  {
    allocate();
  }
}
void d4_array::allocate(
  int hsl, int hsu,
  const index_type& sl, const index_type& sh,
  const index_type& nrl, const index_type& nrh,
  const index_type& ncl,const index_type& nch)
{
  unsigned int ss =
    static_cast<unsigned int>(hsu < hsl ? 0 : hsu - hsl + 1);
  if (ss > 0)
  {
    if ((t = new d3_array[ss]) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
    }
    if ((shape = new four_array_shape(hsl, hsu)) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
    }
    t -= hsl;
    d3_array* pti = t + hsl;
    for (int i=hsl; i<=hsu; i++)
    {
      const index_type& rsl=sl[i];
      const index_type& rsh=sh[i];
      const index_type& rnrl=nrl[i];
      const index_type& rnrh=nrh[i];
      const index_type& rncl=ncl[i];
      const index_type& rnch=nch[i];
      const ad_integer& aa=ad_integer(rsl);
      const ad_integer& bb=ad_integer(rsh);
      pti->allocate(aa,bb,rnrl,rnrh,rncl,rnch);
      ++pti;
      //(*this)(i).allocate(sl[i],sh[i],nrl[i],nrh[i],ncl[i],nch[i]);
    }
  }
  else
  {
    allocate();
  }
}
void d4_array::allocate(
  ad_integer hsl, ad_integer hsu,
  const index_type& sl, const index_type& sh,
  const index_type& nrl, const index_type& nrh)
{
  unsigned int ss =
    static_cast<unsigned int>(hsu < hsl ? 0 : hsu - hsl + 1);
  if (ss > 0)
  {
    if ((t = new d3_array[ss]) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
    }
    if ((shape=new four_array_shape(hsl, hsu)) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
    }
    t -= int(hsl);
    d3_array* pti = t + hsl;
    for (int i = hsl; i <= hsu; ++i)
    {
      const index_type& rsl = sl[i];
      const index_type& rsh = sh[i];
      const index_type& rnrl = nrl[i];
      const index_type& rnrh = nrh[i];
      const ad_integer& aa = ad_integer(rsl);
      const ad_integer& bb = ad_integer(rsh);
      pti->allocate(aa, bb, rnrl, rnrh);
      ++pti;
    }
  }
  else
  {
    allocate();
  }
}

void d4_array::allocate(
  ad_integer hsl,ad_integer hsu,
  const index_type& sl, const index_type& sh)
{
  unsigned int ss =
    static_cast<unsigned int>(hsu < hsl ? 0 : hsu - hsl + 1);
  if (ss > 0)
  {
    if ( (t = new d3_array[ss]) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
      ad_exit(21);
    }
    if ( (shape=new four_array_shape(hsl,hsu)) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
      ad_exit(21);
    }
    t -= int(hsl);
    d3_array* pti = t + hsl;
    for (int i=hsl; i<=hsu; i++)
    {
      const index_type& rsl=sl[i];
      const index_type& rsh=sh[i];
      const ad_integer& aa=ad_integer(rsl);
      const ad_integer& bb=ad_integer(rsh);
      pti->allocate(aa,bb);
      ++pti;
    }
  }
  else
  {
    allocate();
  }
}
void d4_array::allocate(ad_integer hsl, ad_integer hsu)
{
  if (hsl > hsu)
  {
    allocate();
  }
  else
  {
    unsigned int size = static_cast<unsigned int>(hsu - hsl + 1);
    if ((t = new d3_array[size]) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
      ad_exit(1);
    }
    if ((shape = new four_array_shape(hsl, hsu)) == 0)
    {
      cerr << " Error: d4_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
      ad_exit(1);
    }
    t -= static_cast<int>(hsl);
    d3_array* pti = t + hsl;
    for (int i = hsl; i <= hsu; ++i)
    {
      pti->allocate();
      ++pti;
    }
  }
}