dmat.cpp

Go to the documentation of this file.

#include "fvar.hpp"

 struct dvec_ptr_ptr
 {
   void ** m;
 };

dmatrix::dmatrix(int nrl, int nrh, int ncl, int nch)
{
  allocate(nrl, nrh, ncl, nch);
}
void dmatrix::allocate(int nrl, int nrh, int ncl, int nch)
{
  allocate(nrl, nrh);
  dvector* pm = m + nrl;
  for (int i = nrl; i <= nrh; ++i)
  {
    pm->allocate(ncl, nch);
    ++pm;
  }
}
 dmatrix dmatrix::sub(int nrl,int nrh)
 {
   if (allocated(*this))
   {
     dmatrix tmp(nrl,nrh);
     for (int i=nrl; i<=nrh; i++)
     {
       tmp[i].shallow_copy((*this)(i));
     }
     return tmp;
   }
   else
   {
     return *this;
   }
 }

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

void dmatrix::allocate(int nrl, int nrh)
{
  if (nrh < nrl)
  {
    allocate();
  }
  else
  {
    index_min = nrl;
    index_max = nrh;
    if ((m = new dvector[rowsize()]) == 0)
    {
      cerr << " Error allocating in dmatrix::allocate(int, int)\n";
      ad_exit(1);
    }
    if ((shape = new mat_shapex(m)) == 0)
    {
      cerr << " Error allocating in dmatrix::allocate(int, int)\n";
      ad_exit(1);
    }
    m -= rowmin();
  }
}
void dmatrix::allocate(ad_integer nrl, ad_integer nrh)
{
  allocate(static_cast<int>(nrl), static_cast<int>(nrh));
}
double min(const dmatrix& m)
{
  double minimum = min(m(m.rowmin()));
  for (int i = m.rowmin() + 1; i <= m.rowmax(); ++i)
  {
    double value = min(m(i));
    if (value < minimum) minimum = value;
  }
  return minimum;
}
double max(const dmatrix& m)
{
  double maximum = max(m(m.rowmin()));
  for (int i = m.rowmin() + 1; i <= m.rowmax(); ++i)
  {
    double value = max(m(i));
    if (value > maximum) maximum = value;
  }
  return maximum;
}
void dmatrix::allocate(const dmatrix& other)
{
  int min = other.rowmin();
  int max = other.rowmax();
  allocate(min, max);

  dvector* pm = m + min;
  const dvector* pother = &other(min);
  for (int i = min; i <= max; ++i)
  {
    pm->allocate(*pother);
    ++pother;
    ++pm;
  }
}
bool ivector_check(const ivector& vector, int indexmin, int indexmax)
{
  return !(vector.indexmin() == indexmin && vector.indexmax() == indexmax);
}
void dmatrix::allocate(int nrl, int nrh, const ivector& ncl, const ivector& nch)
{
#ifndef OPT_LIB
  if (nrl != ncl.indexmin() || nrh != ncl.indexmax()
      || nrl != nch.indexmin() || nrh != nch.indexmax())
  {
    cerr << "Incompatible array bounds in "
         << "dmatrix(int, int, const ivector&,const ivector&)\n";
    ad_exit(1);
  }
#endif
  allocate(nrl, nrh);
  dvector* pm = m + nrl;
  int* pncl = ncl.get_v() + nrl;
  int* pnch = nch.get_v() + nrl;
  for (int i = nrl; i <= nrh; ++i)
  {
    pm->allocate(*pncl, *pnch);
    ++pncl;
    ++pnch;
    ++pm;
  }
}
void dmatrix::allocate(int nrl, int nrh, int ncl, const ivector& nch)
{
#ifndef OPT_LIB
  if (nrl != nch.indexmin() || nrh != nch.indexmax())
  {
    cerr << "Incompatible array bounds in "
         << "dmatrix(int nrl, int nrh, int ncl, const ivector& nch)\n";
    ad_exit(1);
  }
#endif
  allocate(nrl, nrh);
  dvector* pm = m + nrl;
  const int* pnch = nch.get_v() + nrl;
  for (int i = nrl; i <= nrh; ++i)
  {
    pm->allocate(ncl, *pnch);
    ++pm;
    ++pnch;
  }
}
void dmatrix::allocate(int nrl, int nrh, const ivector& ncl, int nch)
{
#ifndef OPT_LIB
  if (nrl != ncl.indexmin() || nrh != ncl.indexmax())
  {
    cerr << "Incompatible array bounds in "
         << "dmatrix(int nrl, int nrh, int ncl, const ivector& nch)\n";
    ad_exit(1);
  }
#endif
  allocate(nrl, nrh);
  dvector* pm = m + nrl;
  const int* pncl = ncl.get_v() + nrl;
  for (int i = nrl; i <= nrh; ++i)
  {
    pm->allocate(*pncl, nch);
    ++pncl;
    ++pm;
  }
}

/*
 dmatrix::dmatrix(void)
 {
   #ifdef DIAG
     myheapcheck("Entering dmatrix(nrl,nrh,ncl,nch)" );
   #endif

   shape = NULL;
   m=NULL;
 }
*/

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

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

dmatrix::dmatrix(const dmatrix& other)
{
  allocate();
  shallow_copy(other);
}
void dmatrix::shallow_copy(const dmatrix& other)
{
  if (other.shape)
  {
    shape = other.shape;
    ++(shape->ncopies);
    m = other.m;

    index_min = other.index_min;
    index_max = other.index_max;
  }
#ifdef DEBUG
  else
  {
    cerr << "Warning -- Unable to shallow copy an unallocated dmatrix.\n";
  }
#endif
}
dmatrix::~dmatrix()
{
  deallocate();
}
dvector cube(const dvector& vec)
{
  dvector results;
  results.allocate(vec);
  for (int i = results.indexmin(); i <= results.indexmax(); ++i)
  {
    results(i) = cube(vec(i));
  }
  return results;
}
dmatrix cube(const dmatrix& mat)
{
  dmatrix results;
  results.allocate(mat);
  for (int i = results.rowmin(); i <= results.rowmax(); ++i)
  {
    results(i) = cube(mat(i));
  }
  return results;
}
void dmatrix::deallocate()
{
  if (shape)
  {
    if (shape->ncopies)
    {
      --(shape->ncopies);
    }
    else
    {
      m = static_cast<dvector*>(shape->get_pointer());
      delete [] m;
      delete shape;
    }
    allocate();
  }
#if defined(DIAG)
  else
  {
    cerr << "Warning -- Unable to deallocate an unallocated dmatrix.\n";
  }
#endif
}