---

mtaylor.hh

Go to the documentation of this file.

// -*- C++ -*-
#ifndef _MTAY_HH_
#define _MTAY_HH_

#include <iostream>
#include <vector>
#include "dataref"
#include "matrix.hh"
#include "std_except.hh"

class MTayElementIndex;

class MTaylor {
  typedef Matrix MType;

  friend ostream& operator << (ostream&, const MTaylor&);
  friend void uplus(MTaylor&, const MTaylor&);
  friend void uminus(MTaylor&, const MTaylor&);
  friend void mplus(MTaylor&, const MTaylor&,const MTaylor&);
  friend void mminus(MTaylor&, const MTaylor&,const MTaylor&);
  friend void times(MTaylor&, const MTaylor&,const MTaylor&);
  friend void mtimes(MTaylor&, const MTaylor&,const MTaylor&);
  friend void ldivide(MTaylor&, const MTaylor&,const MTaylor&);
  friend void mldivide(MTaylor&, const MTaylor&,const MTaylor&);
  friend void mldividefac(MTaylor&, const MTaylor&,const MTaylor&,
                          MTaylor::MType& fac, vector<int>& piv);
  friend void rdivide(MTaylor&, const MTaylor&,const MTaylor&);
  friend void mrdivide(MTaylor&, const MTaylor&,const MTaylor&);
  friend void transpose(MTaylor&, const MTaylor&);
  friend void ctranspose(MTaylor&, const MTaylor&);
  friend void horzcat(MTaylor&, const vector<MTaylor>&);
  friend void vertcat(MTaylor&, const vector<MTaylor>&);
  friend void horzcat(MTaylor&, const vector<const MTaylor*>&);
  friend void vertcat(MTaylor&, const vector<const MTaylor*>&);
  friend void cross(MTaylor&, const MTaylor&,const MTaylor&);
  friend void sincos(MTaylor&, MTaylor&, const MTaylor&);
  friend void mpower(MTaylor&, const MTaylor&, double);
  friend void sumall(MTaylor&, const MTaylor&);
  friend void exp(MTaylor&, const MTaylor&);
  friend void log(MTaylor&, const MTaylor&);

  friend void subsref(MTaylor&, const MTaylor&, const MTaylor::MType&);
  friend void subsref(MTaylor&, const MTaylor&, const MTaylor::MType&, 
                      const MTaylor::MType&);
  friend void subsasgn(MTaylor&, const MTaylor&, const MTaylor::MType&, 
                       const MTaylor&);
  friend void subsasgn(MTaylor&, const MTaylor&, const MTaylor::MType&, 
                       const MTaylor::MType&,const MTaylor&);

  friend MTaylor::MType size(const MTaylor&);
  friend int size(const MTaylor&,int);
  friend int length(const MTaylor&);
  friend bool isScalar(const MTaylor&);

  friend MTaylor uplus(const MTaylor&);
  friend MTaylor uminus(const MTaylor&);
  friend MTaylor mplus(const MTaylor&,const MTaylor&);
  friend MTaylor mminus(const MTaylor&,const MTaylor&);
  friend MTaylor times(const MTaylor&,const MTaylor&);
  friend MTaylor mtimes(const MTaylor&,const MTaylor&);
  friend MTaylor ldivide(const MTaylor&,const MTaylor&);
  friend MTaylor mldivide(const MTaylor&,const MTaylor&);
  friend MTaylor mldividefac(const MTaylor&,const MTaylor&,
                             MTaylor::MType& fac, vector<int>& piv);
  friend MTaylor rdivide(const MTaylor&,const MTaylor&);
  friend MTaylor mrdivide(const MTaylor&,const MTaylor&);
  friend MTaylor transpose(const MTaylor&);
  friend MTaylor ctranspose(const MTaylor&);
  friend MTaylor horzcat(const vector<MTaylor>&);
  friend MTaylor vertcat(const vector<MTaylor>&);
  friend MTaylor horzcat(const vector<const MTaylor*>&);
  friend MTaylor vertcat(const vector<const MTaylor*>&);
  friend MTaylor cross(const MTaylor&,const MTaylor&);
  friend MTaylor sin(const MTaylor&);
  friend MTaylor cos(const MTaylor&);
  friend MTaylor sincos(MTaylor&, const MTaylor&);
  friend MTaylor mpower(const MTaylor&, double);
  friend MTaylor sumall(const MTaylor&);
  friend MTaylor exp(const MTaylor&);
  friend MTaylor log(const MTaylor&);

  friend MTaylor subsref(const MTaylor&, const MTaylor::MType&);
  friend MTaylor subsref(const MTaylor&, const MTaylor::MType&, 
                         const MTaylor::MType&);
  friend MTaylor subsasgn(const MTaylor&, const MTaylor::MType&, 
                          const MTaylor&);
  friend MTaylor subsasgn(const MTaylor&, const MTaylor::MType&, 
                          const MTaylor::MType&, const MTaylor&);

  friend bool operator == (const MTaylor&, const MTaylor&);
  friend bool operator != (const MTaylor&, const MTaylor&);
  friend bool operator < (const MTaylor&, const MTaylor&);
  friend bool operator > (const MTaylor&, const MTaylor&);
  friend bool operator <= (const MTaylor&, const MTaylor&);
  friend bool operator >= (const MTaylor&, const MTaylor&);

  friend class MTayElementIndex;

public:
  MTaylor() 
    : ref(new MTaylorData) {}
  MTaylor(const MTaylor& mtay)
    : ref(mtay.ref) {}
  MTaylor(const MType& val);
  MTaylor(MType& val);
  MTaylor(const vector<MType>& cv);
  MTaylor(double val);
  MTaylor(int d, const MType& val);
  MTaylor(int rows, int cols, int deg, double *v)
    : ref(new MTaylorData(rows, cols, deg, v)) {}
  ~MTaylor() {}

  MTaylor& operator = (const MTaylor& tay) {
    ref = tay.ref; return *this; }

  bool IsEmpty() const { return ref->la == 0; }
  bool IsScalar() const { return ref->m * ref->n == 1; }
  int Degree() const { return ref->d; }
  void Reserve(int deg);
  void Resize(int deg);
  void Clear();
  void ExtractData(double*) const;

  vector<MType> Taylor2Vector() const;
  MTayElementIndex operator [] (int i);
  MTayElementIndex operator [] (int i) const;

  MTaylor& operator += (const MTaylor& tay) {
    *this = mplus(*this,tay); return *this; }
  MTaylor& operator -= (const MTaylor& tay) {
    *this = mminus(*this,tay); return *this; }

private:

  class MTaylorData {
  public:
    int m, n, d, da, la; // la=m*n*(da+1)
    double *re;

    MTaylorData()
      : m(0), n(0), d(-1), da(-1), la(0), re(NULL) {}
    MTaylorData(int rows, int cols, int deg, double *v);
    MTaylorData(const MTaylorData& mtay);
    MTaylorData& operator = (const MTaylorData& mtay);

    ~MTaylorData();
  };

  DataRef<MTaylorData> ref;

  void Realloc(int deg);
};

class MTayElementIndex {
public:
  MTayElementIndex(MTaylor *tayptr, int k)
    : tay(tayptr), ctay(tayptr), ind(k) {}
  MTayElementIndex(const MTaylor *tayptr, int k)
    : tay(NULL), ctay(tayptr), ind(k) {}
  operator MTaylor::MType() const;
  void operator = (const MTaylor::MType& val);
  void operator = (const MTayElementIndex&);
private:
  MTaylor *tay;
  const MTaylor *ctay;
  int ind;
};

class MTaylorException : public StandardException {
public:
  MTaylorException() : StandardException() {};
  MTaylorException(const string& fname, int ln, 
                   const string& msg = string())
    : StandardException(fname, ln, head()+msg) {}
  const char *head() const {return NULL;};
};

class MTaylorInvalidDimensionException : public StandardException {
public:
  MTaylorInvalidDimensionException() : StandardException() {};
  MTaylorInvalidDimensionException(const string& fname, int ln, 
                                   const string& msg = string())
    : StandardException(fname, ln, head()+msg) {}
  const char *head() const {return "Cannot allocate mtaylor:  Invalid dimension(s).";};
};

class MTaylorEmptyException : public StandardException {
public:
  MTaylorEmptyException() : StandardException() {};
  MTaylorEmptyException(const string& fname, int ln, 
                        const string& msg = string())
    : StandardException(fname, ln, head()+msg) {}
  const char *head() const {return "Illegal operation.  MTaylor is empty.";};
};

class MTaylorRowDimException : public StandardException {
public:
  MTaylorRowDimException() : StandardException() {};
  MTaylorRowDimException(const string& fname, int ln, 
                         const string& msg = string())
    : StandardException(fname, ln, head()+msg) {}
  const char *head() const {return "MTaylor row dimension mismatch";};
};

class MTaylorColDimException : public StandardException {
public:
  MTaylorColDimException() : StandardException() {};
  MTaylorColDimException(const string& fname, int ln, 
                         const string& msg = string())
    : StandardException(fname, ln, head()+msg) {}
  const char *head() const {return "MTaylor column dimension mismatch";};
};

inline MTayElementIndex
MTaylor::operator [] (int i) {
    return MTayElementIndex(this,i); 
}

inline MTayElementIndex
MTaylor::operator [] (int i) const {
    return MTayElementIndex(this,i); 
}

inline int
length(const MTaylor& a) {
  return a.ref->m * a.ref->n;
}

inline bool
isScalar(const MTaylor& a) {
  return (a.ref->m * a.ref->n == 1);
}

inline int
size(const MTaylor& a, int dim) {
  if (dim == 1)
    return a.ref->m;
  else
    return a.ref->n;
}

inline MTaylor
uplus(const MTaylor& a) {
  MTaylor c;
  uplus(c,a);
  return c;
}

inline MTaylor
uminus(const MTaylor& a) {
  MTaylor c;
  uminus(c,a);
  return c;
}

inline MTaylor
mplus(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  mplus(c,a,b);
  return c;
}

inline MTaylor
mminus(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  mminus(c,a,b);
  return c;
}

inline MTaylor
times(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  times(c,a,b);
  return c;
}

inline MTaylor
mtimes(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  mtimes(c,a,b);
  return c;
}

inline MTaylor
ldivide(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  ldivide(c,a,b);
  return c;
}

inline MTaylor
mldivide(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  mldivide(c,a,b);
  return c;
}

inline MTaylor
mldividefac(const MTaylor& a, const MTaylor& b, MTaylor::MType& fac, 
         vector<int>& piv) {
  MTaylor c;
  mldividefac(c,a,b,fac,piv);
  return c;
}

inline MTaylor
rdivide(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  rdivide(c,a,b);
  return c;
}

inline MTaylor
mrdivide(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  mrdivide(c,a,b);
  return c;
}

inline MTaylor
transpose(const MTaylor& a) {
  MTaylor c;
  transpose(c,a);
  return c;
}

inline MTaylor
ctranspose(const MTaylor& a) {
  MTaylor c;
  ctranspose(c,a);
  return c;
}

inline MTaylor
horzcat(const vector<MTaylor>& av) {
  MTaylor c;
  horzcat(c,av);
  return c;
}

inline MTaylor
vertcat(const vector<MTaylor>& av) {
  MTaylor c;
  vertcat(c,av);
  return c;
}

inline MTaylor
horzcat(const vector<const MTaylor*>& av) {
  MTaylor c;
  horzcat(c,av);
  return c;
}

inline MTaylor
vertcat(const vector<const MTaylor*>& av) {
  MTaylor c;
  vertcat(c,av);
  return c;
}

inline MTaylor
cross(const MTaylor& a, const MTaylor& b) {
  MTaylor c;
  cross(c,a,b);
  return c;
}

inline MTaylor
sin(const MTaylor& a) {
  MTaylor s,c;
  sincos(s,c,a);
  return s;
}

inline MTaylor
cos(const MTaylor& a) {
  MTaylor s,c;
  sincos(s,c,a);
  return c;
}

inline MTaylor
sincos(MTaylor& c, const MTaylor& a) {
  MTaylor s;
  sincos(s,c,a);
  return s;
}

inline MTaylor
mpower(const MTaylor& a, double b) {
  MTaylor c;
  mpower(c,a,b);
  return c;
}

inline MTaylor
sumall(const MTaylor& a) {
  MTaylor c;
  sumall(c,a);
  return c;
}

inline MTaylor
exp(const MTaylor& a) {
  MTaylor c;
  exp(c,a);
  return c;
}

inline MTaylor
log(const MTaylor& a) {
  MTaylor c;
  log(c,a);
  return c;
}

inline MTaylor
subsref(const MTaylor& a, const MTaylor::MType& ind) {
  MTaylor c;
  subsref(c,a,ind);
  return c;
}

inline MTaylor
subsref(const MTaylor& a, const MTaylor::MType& ind, 
        const MTaylor::MType& jnd) {
  MTaylor c;
  subsref(c,a,ind,jnd);
  return c;
}

inline MTaylor
subsasgn(const MTaylor& a, const MTaylor::MType& ind, const MTaylor& b) {
  MTaylor c;
  subsasgn(c,a,ind,b);
  return c;
}

inline MTaylor
subsasgn(const MTaylor& a, const MTaylor::MType& ind, 
        const MTaylor::MType& jnd, const MTaylor& b) {
  MTaylor c;
  subsasgn(c,a,ind,jnd,b);
  return c;
}

inline MTaylor
operator + (const MTaylor& a) {
  return uplus(a);
}

inline MTaylor
operator - (const MTaylor& a) {
  return uminus(a);
}

inline MTaylor
operator + (const MTaylor& a, const MTaylor& b) {
  return mplus(a,b);
}

inline MTaylor
operator - (const MTaylor& a, const MTaylor& b) {
  return mminus(a,b);
}

inline MTaylor
operator * (const MTaylor& a, const MTaylor& b) {
  return mtimes(a,b);
}

inline MTaylor
operator / (const MTaylor& a, const MTaylor& b) {
  return mrdivide(a,b);
}

inline bool
operator != (const MTaylor& a, const MTaylor& b) {
  return ! (a == b);
}

inline bool
operator <= (const MTaylor& a, const MTaylor& b) {
  return ! (a > b);
}

inline bool
operator >= (const MTaylor& a, const MTaylor& b) {
  return ! (a < b);
}

#endif

---