d3arr.cpp

Go to the documentation of this file.

#if (__cplusplus > 199711L)
  #include <iterator>
  #include <algorithm>
#endif

#include "fvar.hpp"
#ifdef DEBUG
  #include <cassert>
#endif

double sum(const d3_array& darray)
{
  double total = 0.0;
#if (__cplusplus <= 199711L)
  for (int i = darray.indexmin(); i <= darray.indexmax(); ++i)
  {
    total += sum(darray.elem(i));
  }
#else
  std::for_each(darray.begin(), darray.end(), [&total](dmatrix& matrix) {
    total += sum(matrix);
  });
#endif
  return total;
}
 d3_array d3_array::sub(int nrl,int nrh)
 {
   if (allocated(*this))
   {
     d3_array tmp(nrl,nrh);
     for (int i=nrl; i<=nrh; i++)
     {
       tmp[i].shallow_copy((*this)(i));
     }
     return tmp;
   }
   else
   {
     return *this;
   }
 }
d3_array::d3_array(const d3_array_position& position)
#if defined(__INTEL_COMPILER) || defined(__OPENCC__)
  { allocate(position.indexmin(), position.indexmax()); }
#else
  :d3_array(position.indexmin(), position.indexmax()) { }
#endif

d3_array::d3_array(int nrl, int nrh)
{
  allocate(nrl, nrh);
}

 d3_array::d3_array(int sl,int sh,int nrl,int nrh,int ncl,int nch)
 {
   allocate(sl,sh,nrl,nrh,ncl,nch);
#ifndef OPT_LIB
   initialize();
#endif
 }

d3_array::d3_array(int sl, int sh, int nrl, int nrh, const ivector& ncl,
  int nch)
 {
   allocate(sl,sh,nrl,nrh,ncl,nch);
#ifndef OPT_LIB
   initialize();
#endif
 }

d3_array::d3_array(int sl, int sh, int nrl, int nrh, const ivector& ncl,
  const ivector& nch)
 {
   allocate(sl,sh,nrl,nrh,ncl,nch);
#ifndef OPT_LIB
   initialize();
#endif
 }

void d3_array::allocate(int sl, int sh, int nrl, int nrh,
  const ivector& ncl, int nch)
{
#ifdef DIAG
  myheapcheck("Entering d3_array matrix(sl,sh,nrl,nrh,ncl,nch)" );
#endif
#ifndef OPT_LIB
  if (sl != ncl.indexmin() || sh != ncl.indexmax())
  {
     cerr << "Incompatible d3_array bounds in "
         << __FILE__ << ':' << __LINE__ << '\n';
     ad_exit(1);
  }
#endif
  allocate(sl, sh);

  dmatrix* pti = t + sl;
  int* pncli = ncl.get_v() + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(nrl, nrh, *pncli, nch);
    ++pti;
    ++pncli;
  }
}
void d3_array::allocate(
  int sl, int sh,
  int nrl, int nrh,
  int ncl, const ivector& nch)
{
#ifdef DIAG
  myheapcheck("Entering d3_array matrix(sl,sh,nrl,nrh,ncl,nch)" );
#endif
#ifndef OPT_LIB
  if (sl != nch.indexmin() || sh != nch.indexmax())
  {
     cerr << "Incompatible d3_array bounds in "
         << __FILE__ << ':' << __LINE__ << '\n';
     ad_exit(1);
  }
#endif
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  int* pnchi = nch.get_v() + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(nrl, nrh, ncl, *pnchi);
    ++pti;
    ++pnchi;
  }
}
void d3_array::allocate(int sl, int sh, int nrl, int nrh, int ncl, int nch)
{
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(nrl, nrh, ncl, nch);
    ++pti;
  }
}
void d3_array::allocate(int sl, int sh, int nrl, int nrh)
{
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(nrl, nrh);
    ++pti;
  }
}
void d3_array::allocate(
  int sl, int sh,
  const index_type& nrl,const index_type& nrh)
{
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(nrl(i), nrh(i));
    ++pti;
  }
}
void d3_array::allocate(int sl, int sh)
{
  if (sl > sh)
  {
    return allocate();
  }
  if ((shape = new three_array_shape(sl, sh)) == 0)
  {
    cerr << " Error: d3_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ((t = new dmatrix[slicesize()]) == 0)
  {
    cerr << " Error: d3_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  t -= slicemin();
  dmatrix* pti = t + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate();
    ++pti;
  }
}
void d3_array::allocate(const d3_array& other)
{
  int min = other.slicemin();
  int max = other.slicemax();
  allocate(min, max);
  dmatrix* pti = t + min;
  const dmatrix* potheri = &other(min);
  for (int i = min; i <= max; ++i)
  {
    pti->allocate(*potheri);
    ++pti;
    ++potheri;
  }
}
void d3_array::allocate(
  int sl, int sh,
  const ivector& nrl, const ivector& nrh,
  const imatrix& ncl, const imatrix& nch)
{
  if (sl !=nrl.indexmin() || sh !=nrl.indexmax()
      || sl !=nrh.indexmin() || sh !=nrh.indexmax())
  {
     cerr << "Incompatible d3_array bounds in "
         << __FILE__ << ':' << __LINE__ << '\n';
     ad_exit(1);
  }
  if ((shape = new three_array_shape(sl, sh)) == 0)
  {
    cerr << " Error: d3_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ((t = new dmatrix[slicesize()]) == 0)
  {
    cerr << " Error: d3_array unable to allocate memory in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  t -= sl;
  dmatrix* pti = t + sl;
  const int* pnrhi = nrh.get_v() + sl;
  const int* pnrli = nrl.get_v() + sl;
  const ivector* pnchi = &nch(sl);
  const ivector* pncli = &ncl(sl);
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(*pnrli, *pnrhi, *pncli, *pnchi);
    ++pti;
    ++pnrhi;
    ++pnrli;
    ++pnchi;
    ++pncli;
  }
}
void d3_array::allocate(
  int sl, int sh,
  const ivector& nrl, const ivector& nrh,
  int ncl, int nch)
{
  if (sl != nrl.indexmin() || sh != nrl.indexmax()
      || sl != nrh.indexmin() || sh != nrh.indexmax())
  {
     cerr << "Incompatible d3_array bounds in "
         << __FILE__ << ':' << __LINE__ << '\n';
     ad_exit(1);
  }
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  const int* pnrli = nrl.get_v() + sl;
  const int* pnrhi = nrh.get_v() + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(*pnrli, *pnrhi, ncl, nch);
    ++pti;
    ++pnrli;
    ++pnrhi;
  }
}
void d3_array::allocate(
  int sl, int sh,
  const ivector& nrl, int nrh,
  int ncl, int nch)
{
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  const int* pnrli = nrl.get_v() + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(*pnrli, nrh, ncl, nch);
    ++pti;
    ++pnrli;
  }
}
void d3_array::allocate(
  int sl, int sh,
  int nrl, const ivector& nrh,
  int ncl, int nch)
{
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  const int* pnrhi = nrh.get_v() + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(nrl, *pnrhi, ncl, nch);
    ++pti;
    ++pnrhi;
  }
}
void d3_array::allocate(int sl, int sh, const ivector& nrl, const ivector& nrh,
  int ncl, const imatrix& nch)
 {
   if (sl !=nrl.indexmin() || sh !=nrl.indexmax() ||
       sl !=nrh.indexmin() || sh !=nrh.indexmax())
   {
     cerr << "Incompatible d3_array bounds in "
          << __FILE__ << ':' << __LINE__ << '\n';
     ad_exit(1);
   }
   if ((shape=new three_array_shape(sl, sh)) == 0)
   {
     cerr << " Error: d3_array unable to allocate memory in "
          << __FILE__ << ':' << __LINE__ << '\n';
     ad_exit(1);
   }
   if ( (t = new dmatrix[slicesize()]) == 0)
   {
     cerr << " Error: d3_array unable to allocate memory in "
          << __FILE__ << ':' << __LINE__ << '\n';
     ad_exit(1);
   }
   t -= sl;
   dmatrix* pti = t + sl;
   const int* pnrhi = nrh.get_v() + sl;
   const int* pnrli = nrl.get_v() + sl;
   const ivector* pnchi = &nch(sl);
   for (int i = sl; i <= sh; ++i)
   {
     pti->allocate(*pnrli, *pnrhi, ncl, *pnchi);
     ++pti;
     ++pnrhi;
     ++pnrli;
     ++pnchi;
   }
 }

void d3_array::allocate(
  int sl, int sh,
  int nrl, const ivector& nrh,
  int ncl, const imatrix& nch)
{
  if (sl != nrh.indexmin() || sh != nrh.indexmax())
  {
    cerr << "Incompatible d3_array bounds in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ((shape = new three_array_shape(sl, sh)) == 0)
  {
    cerr << " Error: d3_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  if ( (t = new dmatrix[slicesize()]) == 0)
  {
    cerr << " Error: d3_array unable to allocate memory in "
           << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  t -= slicemin();
  dmatrix* pti = t + sl;
  const int* pnrhi = nrh.get_v() + sl;
  const ivector* pnchi = &nch(sl);
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(nrl, *pnrhi, ncl, *pnchi);
    ++pti;
    ++pnrhi;
    ++pnchi;
  }
}
void d3_array::allocate(int sl, int sh,
  const ivector& nrl, const ivector& nrh,
  const ivector& ncl, const ivector& nch)
{
  if (sl != nrl.indexmin() || sh != nrl.indexmax()
      || sl != nrh.indexmin() || sh != nrh.indexmax())
  {
    cerr << "Incompatible d3_array bounds in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  const int* pnrhi = nrh.get_v() + sl;
  const int* pnrli = nrl.get_v() + sl;
  const int* pnchi = nch.get_v() + sl;
  const int* pncli = ncl.get_v() + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(*pnrli, *pnrhi, *pncli, *pnchi);
    ++pti;
    ++pnrhi;
    ++pnrli;
    ++pnchi;
    ++pncli;
  }
}
void d3_array::allocate(
  int sl, int sh,
  int nrl, int nrh,
  const ivector& ncl, const ivector& nch)
{
  if (sl != ncl.indexmin() || sh != ncl.indexmax()
      || sl != nch.indexmin() || sh != nch.indexmax())
  {
    cerr << "Incompatible d3_array bounds in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  const int* pnchi = nch.get_v() + sl;
  const int* pncli = ncl.get_v() + sl;
  for (int i = slicemin(); i <= slicemax(); ++i)
  {
    pti->allocate(nrl, nrh, *pncli, *pnchi);
    ++pti;
    ++pnchi;
    ++pncli;
  }
}
void d3_array::allocate(int sl, int sh,
  int nrl, const ivector& nrh,
  int ncl, const ivector& nch)
{
  if (sl != nrh.indexmin() || sh != nrh.indexmax()
      || sl != nch.indexmin() || sh != nch.indexmax())
  {
    cerr << "Incompatible d3_array bounds in "
         << __FILE__ << ':' << __LINE__ << '\n';
    ad_exit(1);
  }
  allocate(sl, sh);
  dmatrix* pti = t + sl;
  int* pnrhi = nrh.get_v() + sl;
  int* pnchi = nch.get_v() + sl;
  for (int i = sl; i <= sh; ++i)
  {
    pti->allocate(nrl, *pnrhi, ncl, *pnchi);
    ++pti;
    ++pnrhi;
    ++pnchi;
  }
}

d3_array::d3_array(int sl, int sh, int nrl, const ivector& nrh, int ncl,
  const ivector& nch)
 {
   allocate(sl,sh,nrl,nrh,ncl,nch);
 }

d3_array::d3_array(int sl, int sh, const ivector& nrl, const ivector& nrh,
  const imatrix& ncl, const imatrix& nch)
 {
   allocate(sl,sh,nrl,nrh,ncl,nch);
 }

d3_array::d3_array(int sl, int sh, const ivector& nrl, const ivector& nrh,
  int ncl, const imatrix& nch)
 {
   allocate(sl,sh,nrl,nrh,ncl,nch);
 }

d3_array::d3_array(int sl, int sh, int nrl, const ivector& nrh, int ncl,
  const imatrix& nch)
 {
   allocate(sl,sh,nrl,nrh,ncl,nch);
 }

d3_array::d3_array(int sl, int sh, int nrl, const ivector& nrh, int ncl,
  int nch)
 {
   allocate(sl,sh,nrl,nrh,ncl,nch);
 }

d3_array::d3_array(const d3_array& other)
{
  shallow_copy(other);
}
void d3_array::shallow_copy(const d3_array& other)
{
  if (other.shape)
  {
    shape = other.shape;
    ++(shape->ncopies);
    t = other.t;
  }
  else
  {
#ifdef DEBUG
    cerr << "Warning -- Unable to shallow copy an unallocated d3_array.\n";
#endif
    allocate();
  }
}
void d3_array::initialize()
{
  // only initialize allocated objects
  if (t)
  {
#if (__cplusplus <= 199711L)
    for (int i = slicemin(); i <= slicemax(); ++i)
    {
      t[i].initialize();
    }
#else
    std::for_each(begin(), end(), [](dmatrix& matrix) {
      matrix.initialize();
    });
#endif
  }
}
void d3_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 d3_array.\n";
    //ad_exit(1);
  }
#endif
}
d3_array::~d3_array()
{
  deallocate();
}
three_array_shape::three_array_shape(int sl, int su)
: ncopies(0), slice_min(sl), slice_max(su)
{
}