functions.hh

Go to the documentation of this file.

/*************************************************************************

Copyright Rice University, 2004 - 2015.
All rights reserved.

Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, provided that the above copyright notice(s) and this
permission notice appear in all copies of the Software and that both the
above copyright notice(s) and this permission notice appear in supporting
documentation.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY
RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS
NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL
DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.

Except as contained in this notice, the name of a copyright holder shall
not be used in advertising or otherwise to promote the sale, use or other
dealings in this Software without prior written authorization of the
copyright holder.

**************************************************************************/

#ifndef __RVL_FCNS
#define __RVL_FCNS

#include "utility.hh"
#include "local.hh"

namespace RVL {

  template<class Scalar>
  class RVLCopy: public BinaryLocalFunctionObject<Scalar> {
  private:
    RVLCopy(const RVLCopy<Scalar> &) {}
  public:
    RVLCopy() {}
    ~RVLCopy() {}
  
    using RVL::BinaryLocalEvaluation<Scalar>::operator();
    void operator()(LocalDataContainer<Scalar> & x,
            LocalDataContainer<Scalar> const & y) {
      if (x.getSize() > y.getSize()) {
    RVLException e; e<<"Error: RVLCopy\n";
    e<<"input shorter than output - copy does not cover\n";
    throw e;
      }
      else {
    size_t n = x.getSize();
    Scalar * px = x.getData();
    Scalar const * py = y.getData();
    for (size_t i=0;i<n;i++) {
      px[i]=py[i];
    }
      }
    }
    string getName() const  { return "RVLCopy"; }
  };

  template<>
  class RVLCopy<float>: public BinaryLocalFunctionObject<float> {
  private:
    RVLCopy(const RVLCopy<float> &) {}
  public:
    RVLCopy() {}
    ~RVLCopy() {}
  
    using RVL::BinaryLocalEvaluation<float>::operator();
    void operator()(LocalDataContainer<float> & x,
            LocalDataContainer<float> const & y) {
      if (x.getSize() > y.getSize()) {
    RVLException e; e<<"Error: RVLCopy\n";
    e<<"input shorter than output - copy does not cover\n";
    throw e;
      }
      else {
    float * px = x.getData();
    float const * py = y.getData();
    memcpy(px, py, x.getSize()*sizeof(float));
      }
    }
    string getName() const  { return "RVLCopy<float>"; }
  };

  template<>
  class RVLCopy<double>: public BinaryLocalFunctionObject<double> {
  private:
    RVLCopy(const RVLCopy<double> &) {}
  public:
    RVLCopy() {}
    ~RVLCopy() {}
  
    using RVL::BinaryLocalEvaluation<double>::operator();
    void operator()(LocalDataContainer<double> & x,
            LocalDataContainer<double> const & y) {
      if (x.getSize() > y.getSize()) {
    RVLException e; e<<"Error: RVLCopy\n";
    e<<"input shorter than output - copy does not cover\n";
    throw e;
      }
      else {
    double * px = x.getData();
    double const * py = y.getData();
    memcpy(px, py, x.getSize()*sizeof(double));
      }
    }
    string getName() const  { return "RVLCopy<double>"; }
  };

  template<class Scalar>
  class RVLScale: public UnaryLocalFunctionObject<Scalar> {
  private:
    Scalar c;
    RVLScale();
    RVLScale(const RVLScale<Scalar> &);
  public:
    RVLScale(Scalar _c): c(_c) {}
    ~RVLScale() {}

    using RVL::UnaryLocalEvaluation<Scalar>::operator();
    void operator()(LocalDataContainer<Scalar> & x) {
      size_t n = x.getSize();
      Scalar * px = x.getData();
      for (size_t i=0;i<n;i++) {
    px[i]=c*(px[i]);
      }
    }
    string getName() const  { return "RVLScale"; }
  };

  template<class Scalar>
  class RVLMax
    : public UnaryLocalFunctionObjectScalarRedn<Scalar,Scalar> {
  public:
    RVLMax()
      : UnaryLocalFunctionObjectScalarRedn<Scalar,Scalar>(-numeric_limits<Scalar>::max()) {}
    RVLMax(Scalar _res) 
      : UnaryLocalFunctionObjectScalarRedn<Scalar,Scalar>(_res) {}
    RVLMax(const RVLMax<Scalar> & m)
      : UnaryLocalFunctionObjectScalarRedn<Scalar,Scalar>(m) {}
    ~RVLMax() {}

    // override base class
    void setValue() { ScalarRedn<Scalar>::setValue(-numeric_limits<Scalar>::max()); }
    using RVL::UnaryLocalConstEval<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> const & x) {
      Scalar maxr = ScalarRedn<Scalar>::getValue();
      size_t n=x.getSize();
      Scalar const * px = x.getData();
      for (size_t i=0;i<n;i++) maxr = max<Scalar>(maxr,px[i]);
      ScalarRedn<Scalar>::setValue(maxr);
    }

    string getName() const  { return "RVLMax"; }
    
  };
  
  template<class Scalar>
  class RVLMin: public UnaryLocalFunctionObjectScalarRedn<Scalar,Scalar> {
  public:
    RVLMin()
      : UnaryLocalFunctionObjectScalarRedn<Scalar,Scalar>(numeric_limits<Scalar>::max()) {}
    RVLMin(Scalar _res)
      : UnaryLocalFunctionObjectScalarRedn<Scalar,Scalar>(_res) { }
    RVLMin(const RVLMin<Scalar> & m)
      : UnaryLocalFunctionObjectScalarRedn<Scalar,Scalar>(m) {}
    ~RVLMin() {}
  
    void setValue() { ScalarRedn<Scalar>::setValue(numeric_limits<Scalar>::max()); }
    using RVL::UnaryLocalConstEval<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> const & x) {
      Scalar minr = ScalarRedn<Scalar>::getValue();
      size_t n=x.getSize();
      Scalar const * px = x.getData();
      for (size_t i=0;i<n;i++) { minr=min<Scalar>(minr,px[i]); }
      ScalarRedn<Scalar>::setValue(minr);
    }

    string getName() const  { return "RVLMin"; }
  };

  template<class Scalar>
  class RVLL2innerProd: public BinaryLocalFunctionObjectScalarRedn<Scalar,Scalar> {

  private:
    mutable Scalar scale;
  public:
    RVLL2innerProd(Scalar _scale = ScalarFieldTraits<Scalar>::One(), 
           Scalar _init = ScalarFieldTraits<Scalar>::Zero())
      :  BinaryLocalFunctionObjectScalarRedn<Scalar,Scalar>(_init),
     scale(_scale) {}
    RVLL2innerProd(const RVLL2innerProd<Scalar> & ipfo) 
      :  BinaryLocalFunctionObjectScalarRedn<Scalar,Scalar>(ipfo),
     scale(ipfo.scale) {}
    ~RVLL2innerProd() {}

    void setValue() { ScalarRedn<Scalar>::setValue(ScalarFieldTraits<Scalar>::Zero()); }

    using RVL::BinaryLocalConstEval<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> const & v, 
             LocalDataContainer<Scalar> const & w) {
      try {
    size_t n=v.getSize();
    if (n != w.getSize()) {
      RVLException e; e<<"Error: RVLL2innerProd::operator()\n";
      e<<"operands do not have same dimension\n";
      e<<"\noperand 1:\n";
      v.write(e);
      e<<"\noperand 2:\n";
      w.write(e);
      throw e;
    }
    Scalar const * pv = v.getData();
    Scalar const * pw = w.getData();
    Scalar raw = ScalarFieldTraits<Scalar>::Zero();
    Scalar ip = this->getValue();
    for (size_t i=0;i<n;i++) {
      raw += pv[i]*pw[i];
    }
    ip += scale*raw;
    ScalarRedn<Scalar>::setValue(ip);
      }
      catch (RVLException & e) {
    e<<"\ncalled from RVLL2innerProduct::operator()\n";
    throw e;
      }
    }

    Scalar getScale() const { return scale; }

    void setScale(Scalar newscale) { scale = newscale; }

    string getName() const  { return "basic rvl L2innerProd"; }

  };

  template<class Scalar>
  class RVLL2innerProd<complex<Scalar> >: 
    public BinaryLocalFunctionObjectScalarRedn<complex<Scalar>, complex<Scalar> > {
  private:
    Scalar scale;
  public:
    RVLL2innerProd(Scalar _scale = ScalarFieldTraits<Scalar>::One(), 
           complex<Scalar> _init = complex<Scalar>(ScalarFieldTraits<Scalar>::Zero()))
      : BinaryLocalFunctionObjectScalarRedn<complex<Scalar>, complex<Scalar> > (_init),
     scale(_scale) {}
    RVLL2innerProd(const RVLL2innerProd<complex<Scalar> > & ipfo)
      : BinaryLocalFunctionObjectScalarRedn<complex<Scalar>, complex<Scalar> > (ipfo), 
     scale(ipfo.scale) {}
    ~RVLL2innerProd() {}

    void setValue() { ScalarRedn<complex<Scalar> >::setValue(complex<Scalar>(ScalarFieldTraits<Scalar>::Zero())); }

    using RVL::BinaryLocalConstEval< complex<Scalar> >::operator();
    void operator() (LocalDataContainer<complex<Scalar> > const & v, 
             LocalDataContainer<complex<Scalar> > const & w) {
      try {
    size_t n=v.getSize();
    if (n != w.getSize()) {
      RVLException e; e<<"Error: RVLL2innerProd<complex>::operator()\n";
      e<<"operands do not have same dimension\n";
      e<<"\noperand 1:\n";
      v.write(e);
      e<<"\noperand 2:\n";
      w.write(e);
      throw e;
    }
    complex<Scalar> const * pv = v.getData();
    complex<Scalar> const * pw = w.getData();
    complex<Scalar> raw = complex<Scalar>(ScalarFieldTraits<Scalar>::Zero());
    complex<Scalar> ip = BinaryLocalFunctionObjectScalarRedn<complex<Scalar>, complex<Scalar> >::getValue();
    for (size_t i=0;i<n;i++) {
      raw += pv[i]*conj(pw[i]);
    }
    ip += scale*raw;
    ScalarRedn<complex<Scalar> >::setValue(ip);
      }
      catch (RVLException & e) {
    e<<"\ncalled from RVLL2innerProduct::operator()\n";
    throw e;
      }
    }

    Scalar getScale() const { return scale; }

    void setScale(Scalar newscale) { scale = newscale; }

    string getName() const  { return "basic rvl L2innerProd <complex>"; }

  };

  template<class Scalar>
  class RVLAddAccumulate: public BinaryLocalFunctionObject<Scalar> {
    
  public:
    RVLAddAccumulate() {}
    RVLAddAccumulate(const RVLAddAccumulate<Scalar> &) {}
    virtual ~RVLAddAccumulate() {}

    using RVL::BinaryLocalEvaluation<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> & v,
             LocalDataContainer<Scalar> const & w) {
      size_t n=v.getSize();
      if (n != w.getSize()) {
    RVLException e;
    e<<"Error: RVLAddAccumulate::operator()\n";
    e<<"inputs of different sizes\n";
    e<<"first input:\n";
    v.write(e);
    e<<"second input:\n";
    w.write(e);
    throw e;
      }
      Scalar * pv = v.getData();
      Scalar const * pw = w.getData();
      for (size_t i=0;i<n;i++) {
    pv[i] += pw[i];
      }
    }
    string getName() const  { string sname =  "RVLAddAccumulate"; return sname; }
  };

  template<class Scalar>
  class RVLAssignConst: public UnaryLocalFunctionObject<Scalar> {
  private:
    Scalar c;
    RVLAssignConst();
    RVLAssignConst(const RVLAssignConst<Scalar> &);
  public:
    RVLAssignConst(Scalar _c): c(_c) {}
    ~RVLAssignConst() {}

    using RVL::UnaryLocalEvaluation<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> & v) {
      size_t n = v.getSize();
      Scalar * pv = v.getData();
      for (size_t i=0;i<n;i++) {
    pv[i]=c;
      }
    }

    string getName() const { return "RVLAssignConst"; }

  };

  template<class Scalar>
  class RVLRandomize: public UnaryLocalFunctionObject<Scalar> {
  private:
    Scalar a;
    Scalar w;

    RVLRandomize(const RVLRandomize<Scalar> &) {}
    
  public:
    RVLRandomize(): a(0), w(1) {}
    RVLRandomize( long seed, Scalar _a = 0, Scalar _b = 1 )
      : a(_a), w(_b-_a) { PlantSeeds(seed); }
    ~RVLRandomize() {}
  
    virtual bool readsData(size_t i=0) { return false; }

    using RVL::UnaryLocalEvaluation<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> & v) {
      size_t n = v.getSize();
      Scalar * pv = v.getData();
      if( a == ScalarFieldTraits<Scalar>::Zero()) {
    if( w == ScalarFieldTraits<Scalar>::One()) {
      for (size_t i=0;i<n;i++) {
        //      pv[i]=rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
        pv[i]=Random();
      }
    }
    else {
      for (size_t i=0;i<n;i++) {
        //      pv[i]=w*rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
        pv[i]=w*Random();
      }
    }
      } else if (w == 1) {
    for (size_t i=0;i<n;i++) {
      //      pv[i]=rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One()) + a;
      pv[i]=Random() + a;
    }
      } else {
    for (size_t i=0;i<n;i++) {
      //      pv[i]=w*rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One()) + a;
      pv[i]=w*Random() + a;
    }
      }
    }
    
    string getName() const { return "RVLRandomize"; }
    
  };

  template<class Scalar>
  class RVLRandomize<complex<Scalar> >: public UnaryLocalFunctionObject<complex<Scalar> >{
  private:
    Scalar a;
    Scalar w;

    RVLRandomize(const RVLRandomize<complex<Scalar> >&) {}
    
  public:
    RVLRandomize()
      : a(ScalarFieldTraits<Scalar>::Zero()), w(ScalarFieldTraits<Scalar>::One()) {}
    RVLRandomize( //unsigned int seed, 
         long seed, 
         Scalar _a = ScalarFieldTraits<Scalar>::Zero(), 
         Scalar _b = ScalarFieldTraits<Scalar>::One() )
      : a(_a), w(_b-_a) { //srand(seed); }
      PlantSeeds(seed); }
    ~RVLRandomize() {}
  
    using RVL::UnaryLocalEvaluation< complex<Scalar> >::operator();
    void operator() (LocalDataContainer<complex<Scalar> > & v) {
      size_t n = v.getSize();
      complex<Scalar> * pv = v.getData();
      Scalar rex;
      Scalar imx;
      if( a == ScalarFieldTraits<Scalar>::Zero()) {
    if( w == ScalarFieldTraits<Scalar>::One()) {
      for (size_t i=0;i<n;i++) {
        //rex = rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
        rex = Random();
        //imx = rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
        imx = Random();
        pv[i]=complex<Scalar>(rex,imx);
      }
    }
    else {
      for (size_t i=0;i<n;i++) {
        //rex = w*rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
        rex = w*Random();
        //imx = w*rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
        imx = w*Random();
        pv[i]=complex<Scalar>(rex,imx);
      }
    }
      } else if (w == 1) {
    for (size_t i=0;i<n;i++) {
      //rex = a + rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
      rex = a + Random();
      //imx = a + rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
      imx = a + Random();
      pv[i]=complex<Scalar>(rex,imx);
    }
      } else {
    for (size_t i=0;i<n;i++) {
      //rex = a + w*rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
      rex = a + w*Random();
      //imx = a + w*rand()/(RAND_MAX+ScalarFieldTraits<Scalar>::One());
      imx = a + w*Random();
      pv[i]=complex<Scalar>(rex,imx);
    }
      }
    }
    
    string getName() const { return "RVLRandomize<complex>"; }
    
  };

  template<>
  class RVLRandomize<int>: public UnaryLocalFunctionObject<int> {
  private:
    RVLRandomize(const RVLRandomize<int> &) {}
  public:
    RVLRandomize() {}
    ~RVLRandomize() {}

    using RVL::UnaryLocalEvaluation<int>::operator();
    virtual bool readsData(size_t i=0) { return false; }

    void operator() (LocalDataContainer<int> & v) {
      size_t n = v.getSize();
      int * pv = v.getData();
      for (size_t i=0;i<n;i++) {
    pv[i]=rand();
      }
    }

    string getName() const  { return "RVLRandomize<int>"; }
  };

  template<class Scalar>
  class ASCIIReader: public UnaryLocalFunctionObject<Scalar> {
  private:
    ifstream from;
    string file;
  
    ASCIIReader() {}
    ASCIIReader(const ASCIIReader &) {}
  public:
    ASCIIReader(string fname): from(fname.c_str(),ios::in),file(fname) {
      if (!from) {
    RVLException e; e<<"Error: ASCIIReader constructor\n";
    e<<"failed to open file "<<fname<<"\n";
    throw e;
      }
      // skip the first line as it is simply used to indicate dimension
      string l;
      getline(from,l);
    }
    ~ASCIIReader() { from.close(); }

    using RVL::UnaryLocalEvaluation<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> & v) {
      size_t n = v.getSize();
      Scalar * pv = v.getData();
      size_t i=0;
      while (i<n && from>>pv[i]) {
    // cout<<"i = "<<i<<" v = "<<pv[i]<<endl;
    i++;
      }
      if (i!=n) {
    RVLException e; e<<"Error: ASCIIReader::operator()\n";
    e<<"failed to read "<<n<<" scalars from file \""<<file<<"\"\n";
    throw e;
      }
    }

    string getName() const  { return "ASCIIReader"; }
  };

  template<class Scalar>
  class ASCIIWriter: public UnaryLocalFunctionObjectConstEval<Scalar> {
  private:
    ofstream to;
    string file;
    ASCIIWriter() {}
    ASCIIWriter(const ASCIIWriter &) {}
  public:
    ASCIIWriter(string fname): to(fname.c_str(),ios::out), file(fname) {
      if (!to) {
    RVLException e; e<<"Error: ASCIIWriter constructor\n";
    e<<"failed to open file "<<fname<<"\n";
    throw e;
      }
    }
    ~ASCIIWriter() { to.close(); }

    using RVL::UnaryLocalConstEval<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> const & v) {
      size_t n = v.getSize();
      to<<n<<endl;
      Scalar const * const pv = v.getData();
      to.setf(ios::right,ios::adjustfield);
      size_t i=0;
      while (i<n && to<<pv[i]<<endl) i++;
      if (i!=n) {
    RVLException e; e<<"Error: ASCIIWriter::operator()\n";
    e<<"failed to write "<<n<<" scalars to file \""<<file<<"\"\n";
    throw e;
      }
      to.flush();
    }

    string getName() const  { return "ASCIIWriter"; }
  };

  template<class Scalar>
  class BinaryReader: public UnaryLocalFunctionObjectConstEval<Scalar> {
  private:
    FILE * fp;
    string file;
    mutable long first;
    BinaryReader(){}
    BinaryReader(const BinaryReader<Scalar> &){}
  public:

    BinaryReader(char const * fname, long _first=0L)
      : file(fname), first(_first) {
      if (!(fp=fopen(fname,"r+"))) {
    RVLException e; e<<"Error: BinaryReader constructor\n";
    e<<"failed to open file "<<fname<<"\n";
    throw e;
      }
      if (fseek(fp,first,SEEK_SET)) {
    RVLException e; e<<"Error: BinaryReader constructor\n";
    e<<"failed to seek to specified first position = "<<first<<" in file "<<fname<<"\n";
    throw e;
      }
    }
    BinaryReader(const string fname, long _first=0L)
      : file(fname), first(_first) {
      if (!(fp=fopen(fname.c_str(),"r+"))) {
    RVLException e; e<<"Error: BinaryReader constructor\n";
    e<<"failed to open file "<<fname<<"\n";
    throw e;
      }
      if (fseek(fp,first,SEEK_SET)) {
    RVLException e; e<<"Error: BinaryReader constructor\n";
    e<<"failed to seek to specified first position = "<<first<<" in file "<<fname<<"\n";
    throw e;
      }
    }
    ~BinaryReader() {fclose(fp);}

    bool seek(size_t firstword) {
      first = firstword * sizeof(Scalar);
      if (fseek(fp,first,SEEK_SET)) return false;
      return true;
    }

    using RVL::UnaryLocalEvaluation<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> & v) {
      size_t n = v.getSize();
      Scalar * pv = v.getData();
      if (n != fread(pv,sizeof(Scalar),n,fp)) {
    RVLException e; e<<"Error: BinaryReader::operator()\n";
    e<<"failed to read "<<n<<" "<<sizeof(Scalar)<<"-byte words starting at position "<<first<<" from file "<<file<<"\n";
      throw e;
      }
    }

    string getName() const  { return "BinaryReader"; }

  };

  template<class Scalar>
  class BinaryWriter: public UnaryLocalConstEval<Scalar> {
  private:
    FILE * fp;
    string file;
    long first;
    BinaryWriter(){}
    BinaryWriter(const BinaryWriter<Scalar> &){}
  public:

    BinaryWriter(char const * fname, long _first=0L)
      : file(fname), first(_first) {
      if (!(fp=fopen(fname,"w+"))) {
    RVLException e; e<<"Error: BinaryWriter constructor\n";
    e<<"failed to open file "<<fname<<"\n";
    throw e;
      }
      if (fseek(fp,first,SEEK_SET)) {
    RVLException e; e<<"Error: BinaryWriter constructor\n";
    e<<"failed to seek to specified first position = "<<first<<" in file "<<fname<<"\n";
    throw e;
      }
    }

    BinaryWriter(const string fname, long _first=0L): 
      file(fname), first(_first) {
      if (!(fp=fopen(fname.c_str(),"w+"))) {
    RVLException e; e<<"Error: BinaryWriter constructor\n";
    e<<"failed to open file "<<fname<<"\n";
    throw e;
      }
      if (fseek(fp,first,SEEK_SET)) {
    RVLException e; e<<"Error: BinaryWriter constructor\n";
    e<<"failed to seek to specified first position = "<<first<<" in file "<<fname<<"\n";
    throw e;
      }
    }
    ~BinaryWriter() {fclose(fp);}

    bool seek(size_t firstword) {
      first = firstword * sizeof(Scalar);
      if (fseek(fp,first,SEEK_SET)) return false;
      return true;
    }

    using RVL::UnaryLocalConstEval<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> const & v) {
      size_t n = v.getSize();
      Scalar const * pv = v.getData();
      if (n != fwrite(pv,sizeof(Scalar),n,fp)) {
    RVLException e; e<<"Error: BinaryWriter::operator()\n";
    e<<"failed to write "<<n<<" "<<sizeof(Scalar)<<"-byte words starting at position "<<first<<" to file "<<file<<"\n";
    throw e;
      }
      fflush(fp);
    }

    string getName() const  { return "BinaryWriter"; }

  };

  template<class Scalar>
  class RVLBoxMaxStep: 
    public QuaternaryLocalFunctionObjectScalarRedn<Scalar, Scalar> {

  public:
    RVLBoxMaxStep() {}
    RVLBoxMaxStep(const RVLBoxMaxStep<Scalar> & b) {}
    ~RVLBoxMaxStep() {}

    using RVL::QuaternaryLocalConstEval<Scalar>::operator();    
    void operator()(LocalDataContainer<Scalar> const & x,
            LocalDataContainer<Scalar> const & dx,
            LocalDataContainer<Scalar> const & xmin,
            LocalDataContainer<Scalar> const & xmax) {
      try {
    size_t n = x.getSize();
    if ((n != dx.getSize()) ||
        (n != xmin.getSize()) ||
        (n != xmax.getSize())) {
      RVLException e;
      e<<"Error: GridMaxStep::operator()\n";
      e<<"incompatible inputs\n";
      throw e;
    }
    
    Scalar step = numeric_limits<Scalar>::max();
    Scalar const * px    = x.getData();
    Scalar const * pdx   = dx.getData();
    Scalar const * pxmin = xmin.getData();
    Scalar const * pxmax = xmax.getData();
    
    for (size_t i=0;i<n;i++) {
      Scalar tmp;
      if (pdx[i] > ScalarFieldTraits<Scalar>::Zero() && 
          !ProtectedDivision<Scalar>(pxmax[i]-px[i],pdx[i],tmp)) {
        if (tmp > ScalarFieldTraits<Scalar>::Zero()) { step = min(step,tmp); }
        else { step = ScalarFieldTraits<Scalar>::Zero(); }
      }
      if (pdx[i] < ScalarFieldTraits<Scalar>::Zero() &&
          !ProtectedDivision<Scalar>(pxmin[i]-px[i],pdx[i],tmp)) {
        if (tmp > ScalarFieldTraits<Scalar>::Zero()) { step = min(step,tmp); }
        else { step = ScalarFieldTraits<Scalar>::Zero(); }
      }
    }
    setValue(step);
      }
      catch (RVLException & e) {
    e<<"\ncalled from RVLBoxMaxStep::operator()\n";
    throw e;
      }
    }
   
    string getName() const { string tmp = "RVLBoxMaxStep"; return tmp; }
  };

  template<class Scalar>
  class ElementwiseMultiply: public TernaryLocalFunctionObject<Scalar> {

  public:
    ElementwiseMultiply() {}
    ~ElementwiseMultiply() {}
    
    using RVL::TernaryLocalEvaluation<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> & u, 
             LocalDataContainer<Scalar> const & v, 
             LocalDataContainer<Scalar> const & w) {
      Scalar * up = u.getData(); 
      Scalar const * vp = v.getData();
      Scalar const * wp = w.getData();
      size_t size = u.getSize();
      if( (size > v.getSize())||(size > w.getSize())) {
    RVLException e;
    e << "Error in ElementwiseMultiply::operator() - target LDC has more elements than source LDC.\n";
    throw e;
      }
      for(size_t i = 0; i < size; i++)
    up[i] = vp[i]*wp[i];
    }

    string getName() const  { string s = "ElementwiseMultiply"; return s; }
  };

  template<class Scalar>
  class ElementwiseDivision: public TernaryLocalFunctionObject<Scalar> {

  public:
    ElementwiseDivision() {}
    ~ElementwiseDivision() {}
    
    using RVL::TernaryLocalEvaluation<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> & u, 
             LocalDataContainer<Scalar> const & v, 
             LocalDataContainer<Scalar> const & w) {
      Scalar * up = u.getData();
      Scalar const * vp = v.getData(); 
      Scalar const * wp = w.getData();
      size_t size = u.getSize();
      if( (size > v.getSize())||(size > w.getSize())) {
    RVLException e;
    e << "Error: ElementwiseDivision::operator()\n";
    e << "target LDC has more elements than source LDC.\n";
    throw e;
      }
      for(size_t i = 0; i < size; i++) {
    if (ProtectedDivision<Scalar>(vp[i],wp[i],up[i])) {
      RVLException e;
      e<<"Error: ElementwiseDivision::operator()\n";
      e<<"zerodivide in element "<<i<<"\n";
      throw e;
    }
      }
    }

    string getName() const  { string s = "ElementwiseDivision"; return s; }
  };

  template<class Scalar>
  class ElementwiseSqrtAbs: public BinaryLocalFunctionObject<Scalar> {

  public:
    ElementwiseSqrtAbs() {}
    ~ElementwiseSqrtAbs() {}
    
    using RVL::BinaryLocalEvaluation<Scalar>::operator();
    void operator() (LocalDataContainer<Scalar> & u, 
             LocalDataContainer<Scalar> const & v) {
      Scalar * up = u.getData();
      Scalar const * vp = v.getData(); 
      size_t size = u.getSize();
      if (size > v.getSize()) {
    RVLException e;
    e << "Error: ElementwiseSqrtAbs::operator()\n";
    e << "target LDC has more elements than source LDC.\n";
    throw e;
      }
      for(size_t i = 0; i < size; i++) {
    up[i]=sqrt(abs(vp[i]));
      }
    }

    string getName() const  { string s = "ElementwiseSqrtAbs"; return s; }
  };

}
      
#endif



  

---