//
//  Matrix.h
//
//  Qdos2 - base numeric matrix class interface
//
//  Janet Barnes
//  Sept 1999
//

#ifndef QDOS_MATRIX_H
#define QDOS_MATRIX_H

#include <cmath>
#include <algorithm>
#include <functional>
#include <climits>
#include <memory>
#include <vector>
#include <valarray>
#include <string>
#ifdef HAVE_SSTREAM_H
#include <sstream>
#else
#include <strstream>
#endif
#include <iomanip>
#include <typeinfo>
#include <complex>
#include "Promote.h"
#include "Exception.h"
#include "Viostream.h"
#include "Var.h"
#include "MatrixException.h"
#include "Dim.h"
#include "Sequence.h"
#include "DummyArg.h"

using namespace std;

namespace QDOS {

template<class _ForwardIterator, class T>
void iota(_ForwardIterator __first, _ForwardIterator __last, T value)
{
  while(__first!=__last)
    *__first++=value++;
}

typedef enum MatrixCopyOption_e {DEEP,SHALLOW} MatrixCopyOption;    

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class Void

// Needs more work & thought...to overcome problems in mem_fun_ref where
// a result is required and compiler does not like return foo() where
// void foo(); i.e. cannot return voids at all
//-----------------------------------------------------------------------------
class Void
{
public:
    Void(){;}
};
// ----------- End of class Void -------------------
// --- Forward class declarations
template<class T> class MatrixRes; 
template<class T> class MatrixSeq; 
template<class T> class MatrixResSeq;
template<class T> class MatrixSubSeq;
template<class T> class ScalarSeq;
template<class T> class Scalar;
template<class T> class Matrix;
template<class T, class ContainerT> class MatrixBaseSeq;

// --- Matrix Non-member operators declarations
template<class T>
MatrixRes<T> operator!(const Matrix<T>& m);
template<class T>
MatrixRes<complex<T> > operator!(const Matrix<complex<T> >& m);
template<class T> 
Scalar<typename promote_trait<T,int>::T_promote> operator+(const int& i,\
    const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,int>::T_promote> operator+(const int& i,\
    const Matrix<T>& m);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> operator+(const double & d,\
    const Scalar<T>& s);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> \
    operator+(const double & d,\
    const Matrix<T>& m);
template<class U, class T>
Scalar<typename promote_trait<T,complex<U> >::T_promote> operator+\
    (const complex<U>& c,const Scalar<T>& s);
template<class U, class T>
MatrixRes<typename promote_trait<T,complex<U> >::T_promote> operator+\
    (const complex<U>& c,const Matrix<T>& m);
template<class T>
Scalar<complex<T> > operator+(const Scalar<complex<T> >& s,const double& t);
template<class T>
Scalar<complex<T> > operator+(const Scalar<complex<T> >& s,const int& i);
template<class T>
MatrixRes<complex<T> > operator+(const Matrix<complex<T> >& m,const double& t);
template<class T>
MatrixRes<complex<T> > operator+(const Matrix<complex<T> >& m,const int& i);

template<class T> 
Scalar<typename promote_trait<T,int>::T_promote> operator-(const int& i,\
    const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,int>::T_promote> operator-(const int& i,\
    const Matrix<T>& m);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> operator-(const double & d,\
    const Scalar<T>& s);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> \
    operator-(const double & d,\
    const Matrix<T>& m);
template<class U, class T>
Scalar<typename promote_trait<T,complex<U> >::T_promote> operator-\
    (const complex<U>& c,const Scalar<T>& s);
template<class U, class T>
MatrixRes<typename promote_trait<T,complex<U> >::T_promote> operator-\
    (const complex<U>& c,const Matrix<T>& m);
template<class T>
Scalar<complex<T> > operator-(const Scalar<complex<T> >& s,const double& t);
template<class T>
Scalar<complex<T> > operator-(const Scalar<complex<T> >& s,const int& i);
template<class T>
MatrixRes<complex<T> > operator-(const Matrix<complex<T> >& m,const double& t);
template<class T>
MatrixRes<complex<T> > operator-(const Matrix<complex<T> >& m,const int& i);

template<class U,class T>
MatrixRes<typename promote_trait<T,U>::T_promote> operator*(const vector<U>& v,\
    const Matrix<T>& m);
template<class T> Scalar<typename promote_trait<T,int>::T_promote> operator*\
(const int& i,const Scalar<T>& s);
template<class T> MatrixRes<typename promote_trait<T,int>::T_promote> operator*\
(const int& i,const Matrix<T>& m);
template<class T> Scalar<typename promote_trait<T,double>::T_promote> operator*\
(const double & d,const Scalar<T>& s);
template<class T> MatrixRes<typename promote_trait<T,double>::T_promote> operator*\
(const double & d,const Matrix<T>& m);
template<class U, class T> Scalar<typename promote_trait<T,complex<U> >::T_promote> \
operator*(const complex<U>& c,const Scalar<T>& s);
template<class U, class T> MatrixRes<typename promote_trait<T,complex<U> >::T_promote> \
operator*(const complex<U>& c,const Matrix<T>& m);
template<class T>
Scalar<complex<T> > operator*(const Scalar<complex<T> >& s,const double& t);
template<class T>
Scalar<complex<T> > operator*(const Scalar<complex<T> >& s,const int& i);
template<class T>
MatrixRes<complex<T> > operator*(const Matrix<complex<T> >& m,const double& t);
template<class T>
MatrixRes<complex<T> > operator*(const Matrix<complex<T> >& m,const int& i);

template<class T>
Scalar<typename promote_trait<T,int>::T_promote > operator/(const Scalar<T>&,const int&);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote > operator/(const Scalar<T>&,const double&);
template<class T>
Scalar<complex<T> > operator/(const Scalar<complex<T> >& s,const double& d);
template<class T>
Scalar<complex<T> > operator/(const Scalar<complex<T> >& m,const int& i);
template<class T>
MatrixRes<typename promote_trait<T,int>::T_promote > operator/(const Matrix<T>&,const int&);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote > operator/(const Matrix<T>&,const double&);
template<class T>
MatrixRes<complex<T> > operator/(const Matrix<complex<T> >& m,const double& d);
template<class T>
MatrixRes<complex<T> > operator/(const Matrix<complex<T> >& m,const int& i);

template<class T> ostream& operator<< (ostream& os,const Scalar<T>& s);  
template<class T> ostream& operator<< (ostream& os,const Matrix<T>& m);  
template<class T> istream& operator>> (istream& is,const Matrix<T>& m);
template<class T> Vostream& operator<< (Vostream&,\
  const Matrix<T>&);
template<class T> Vistream& operator>> (Vistream&,Matrix<T>&); 

// ---- Matrix Non-member maths functions declarations
template<class T> Scalar<T> abs(const Scalar<T>& s);
template<class T> MatrixRes<T> abs(const Matrix<T>& m);
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> acos(const Scalar<T>& m);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> acos(const Matrix<T>& m);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> asin(const Scalar<T>& s);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> asin(const Matrix<T>& m);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> atan(const Scalar<T>& s);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> atan(const Matrix<T>& m);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> atan2(const Scalar<T>& s1,\
     const Scalar<T>& s2);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> atan2(const Matrix<T>& m1,\
     const Matrix<T>& m2);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> atan2(const Matrix<T>& m,\
     const Scalar<T>& s); 
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> atan2(const Matrix<T>& m,\
     const T& t); 
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> atan2(const Scalar<T>& s,\
     const Matrix<T>& m);    
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> atan2(const T& t,\
     const Matrix<T>& m);    
template<class T>
Scalar<typename promote_trait<T,float_complex>::T_promote> conj(const Scalar<T>& s);    
template<class T>
MatrixRes<typename promote_trait<T,float_complex>::T_promote> conj(const Matrix<T>& m);    
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> cos(const Scalar<T>& s);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> cos(const Matrix<T>& m);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> cosh(const Scalar<T>& s);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> cosh(const Matrix<T>& m);
template<class T>
ScalarSeq<double> eigenvalues(const Matrix<T>& m);
template<class T> 
MatrixResSeq<double> eigenvectors(const Matrix<T>& m); 
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> exp(const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> exp(const Matrix<T>& m);
MatrixRes<double> hilbert(int n);
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> log(const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> log(const Matrix<T>& m);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> log10(const Scalar<T>& s);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> log10(const Matrix<T>& m);
template<class T> T max(const Matrix<T>& m);
template<class T> T min(const Matrix<T>& m);
template<class T> T norm_svd(const Matrix<T>& m); 
template<class T> T norm(const Scalar<T>& s);
template<class T> T norm(const Matrix<T>& m);
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> pow(const Scalar<T>& s,const T& t);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const Matrix<T>& m,\
  const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const Matrix<T>& m,const T& t);
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> pow(const T& t,const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const Scalar<T>& s,\
  const Matrix<T>& m);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const T& t,const Matrix<T>& m);
template<class U,class T>
Scalar<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> pow(const Scalar<T>& s1, \
    const Scalar<U>& s2);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const Matrix<T>& m1, \
    const Matrix<T>& m2);
template<class U, class T> 
MatrixRes<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const Scalar<U>& s,const Matrix<T>& m);
template<class U, class T> 
MatrixRes<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const Matrix<U>& m, const Scalar<T>& s);
template<class U, class T> 
MatrixRes<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const Matrix<U>& m1,const Matrix<T>& m2);
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> sin(const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> sin(const Matrix<T>& m);
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> sinh(const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> sinh(const Matrix<T>& m);
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> sqrt(const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> sqrt(const Matrix<T>& m);
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> tan(const Scalar<T>& s);
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> tan(const Matrix<T>& m);
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> tan(const Scalar<T>& s);
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> tanh(const Matrix<T>& m);

// MatrixBaseSeq Non-Member operator declarations
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator+\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<int>& v); 
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator+\
    (const vector<int>& v, const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator+\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<double>& v); 
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator+\
    (const vector<double>& v, const MatrixBaseSeq<T, ContainerT >& ms); 
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator+\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<complex<U> >& v); 
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator+\
    (const vector<complex<U> >& v, const MatrixBaseSeq<T, ContainerT >& ms); 
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator+\
    (const int& i,const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator+\
    (const double& d,const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator+\
    (const complex<U>& c,const MatrixBaseSeq<T, ContainerT >& ms);  
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,U>::T_promote> operator+\
(const Scalar<U>& s,const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,U>::T_promote> operator+\
(const Matrix<U>& u,const MatrixBaseSeq<T, ContainerT >& ms);

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator-\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<int>& v); 
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator-\
    (const vector<int>& v, const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator-\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<double>& v); 
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator-\
    (const vector<double>& v, const MatrixBaseSeq<T, ContainerT >& ms); 
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator-\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<complex<U> >& v); 
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator-\
    (const vector<complex<U> >& v, const MatrixBaseSeq<T, ContainerT >& ms); 
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator-\
    (const int& i,const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator-\
    (const double& d,const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator-\
    (const complex<U>& c,const MatrixBaseSeq<T, ContainerT >& ms);  
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,U>::T_promote> operator-\
(const Scalar<U>& u,const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,U>::T_promote> operator-\
(const Matrix<U>& u,const MatrixBaseSeq<T, ContainerT >& ms);

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator*\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<int>& v); 
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator*\
    (const vector<int>& v, const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator*\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<double>& v); 
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator*\
    (const vector<double>& v, const MatrixBaseSeq<T, ContainerT >& ms); 
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator*\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<complex<U> >& v); 
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator*\
    (const vector<complex<U> >& v, const MatrixBaseSeq<T, ContainerT >& ms); 
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator*\
    (const int& i,const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator*\
    (const double& d,const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator*\
    (const complex<U>& c,const MatrixBaseSeq<T, ContainerT >& ms);  
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,U>::T_promote> operator*\
(const Scalar<U>& s,const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,U>::T_promote> operator*\
(const Matrix<U>& u,const MatrixBaseSeq<T, ContainerT >& ms);

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator/\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<int>& v); 
template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator/\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<double>& v); 
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote> operator/\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<complex<U> >& v);  
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,U>::T_promote> operator/\
(const Scalar<U>& u,const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,U>::T_promote> operator/\
(const Matrix<U>& u,const MatrixBaseSeq<T, ContainerT >& ms);

template<class T, class ContainerT > 
MatrixResSeq<T> abs(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    acos(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    asin(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    atan(const MatrixBaseSeq<T, ContainerT >& ms); 
template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const T& t,\
    const MatrixBaseSeq<T,ContainerT>& ms2); 
template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const Scalar<T>& s,\
    const MatrixBaseSeq<T,ContainerT>& ms2); 
template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const Matrix<T>& m,\
    const MatrixBaseSeq<T,ContainerT>& ms2); 
template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const \
    MatrixBaseSeq<T, ContainerT>& ms,const T& t); 
template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const \
    MatrixBaseSeq<T, ContainerT>& ms,const Scalar<T>& s); 
template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const \
    MatrixBaseSeq<T, ContainerT>& ms,const Matrix<T>& m); 
template<class U,class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const \
    MatrixBaseSeq<T, ContainerT>& ms,const vector<U>& v); 
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    cos(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    cosh(const MatrixBaseSeq<T, ContainerT >& ms);   
template<class U, class T, class ContainerT >
MatrixResSeq<typename promote_trait<T,U>::T_promote> cross(const \
    MatrixBaseSeq<T,ContainerT >& ms, const vector<Matrix<U> >& v);
template<class U, class T, class ContainerT >
MatrixResSeq<typename promote_trait<T,U>::T_promote> cross(const \
    MatrixBaseSeq<T,ContainerT >& ms, const Matrix<U>& m);
template<class U, class T, class ContainerT >
MatrixResSeq<typename promote_trait<T,U>::T_promote> cross(const Matrix<U>& m,\
    MatrixBaseSeq<T,ContainerT >& ms); 
template<class U, class T, class ContainerT >
MatrixResSeq<typename promote_trait<T,U>::T_promote> dot(const \
    MatrixBaseSeq<T,ContainerT >& ms, const vector<Matrix<U> >& v);
template<class U, class T, class ContainerT >
MatrixResSeq<typename promote_trait<T,U>::T_promote> dot(const \
    MatrixBaseSeq<T,ContainerT >& ms, const Matrix<U>& m);
template<class U, class T, class ContainerT >
MatrixResSeq<typename promote_trait<T,U>::T_promote> dot(const Matrix<U>& m,\
    MatrixBaseSeq<T,ContainerT >& ms); 
template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    exp(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    pinverse(const MatrixBaseSeq<T, ContainerT >& ms);  
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
log(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
log10(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<T> norm(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT >  
MatrixResSeq<T> pow(const MatrixBaseSeq<T, ContainerT >& ms, const T& t);
template<class U, class T, class ContainerT >  
MatrixResSeq<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const MatrixBaseSeq<T, ContainerT >& ms, const Scalar<U>& s);  
template<class U, class T, class ContainerT >  
MatrixResSeq<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const MatrixBaseSeq<T, ContainerT >& ms, const Matrix<U>& m);  
template<class T, class ContainerT >   
MatrixResSeq<T> pow(const T& t, const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const Scalar<U>& s, const MatrixBaseSeq<T, ContainerT >& ms);
template<class U, class T, class ContainerT >   
MatrixResSeq<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const Matrix<U>& m, const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT >  
MatrixResSeq<T> pow(const MatrixBaseSeq<T, ContainerT >& ms1, \
    const vector<T>& v);
template<class U, class T, class ContainerT >  
MatrixResSeq<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const MatrixBaseSeq<T, ContainerT >& ms, \
    const vector<Matrix<U> >& v);
template<class U, class T, class ContainerT >  
MatrixResSeq<typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote> \
    pow(const MatrixBaseSeq<T, ContainerT >& ms, \
    const vector<Scalar<U> >& v); 
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    sin(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> sinh(\
    const MatrixBaseSeq<T, ContainerT >& ms);   
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> sqrt(const \
MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    tan(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    tanh(const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
ostream& operator<<(ostream& os, const MatrixBaseSeq<T, ContainerT >& ms);
template<class T, class ContainerT > 
istream& operator>> (istream& is, const MatrixBaseSeq<T, ContainerT >& ms);

// +++++++++++++++++++++++++++++++++++++++++++++++++++++ class ValarrayObj
// "Cover" class for valarray<T>. Enables smart pointers to valarray<T>
// objects to be set up
//
//------------------------------------------------------------------------------
template<class T>
class ValarrayObj:public VarTarget, public valarray<T>
{
public:
// --- Constructors     
     ValarrayObj():valarray<T>(){;}
     ValarrayObj(size_t sz):valarray<T>(sz){;}
     ValarrayObj(const T& val,size_t sz):valarray<T>(val,sz){;}
     ValarrayObj(const T* parr, size_t sz):valarray<T>(parr,sz){;}
     ValarrayObj(const valarray<T>& v):valarray<T>(v){;}
     ValarrayObj(const slice_array<T>& s):valarray<T>(s){;}
     ValarrayObj(const gslice_array<T>& g):valarray<T>(g){;}
     ValarrayObj(const mask_array<T>& m):valarray<T>(m){;}
     ValarrayObj(const indirect_array<T>& i):valarray<T>(i){;}
    ~ValarrayObj(){;} 
};
// ----------- End of class ValarrayObj -------------------

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix

// Numeric Matrix class of n-dimensions n>=0
// Matrix with rank=0 corresponds to a T
// Matrix with rank=1 corresponds to a vector<T>
//------------------------------------------------------------------------------
template<class T>
class Matrix 
{
    typedef Var<ValarrayObj<T> > ValarrayObj_var;
    typedef T value_type;
    
    friend class MatrixRes<T>; 
    
 protected:
  
    mutable ValarrayObj_var _val; // data - ptr to a Valarray_objT
    mutable Dim_var _diminfo;  // ptr to dimension info for Matrix
    bool _subarray;               // 1 if subarray, 0 otherwise(default) 
    size_t _dimno;                // which dimension for [] style indexing
    MatrixCopyOption _copyoption; // Deep or shallow copies - deep default
    mutable bool _not_empty;      // true if matrix initialised, false othewise
             
    void createval()
    {
    // Create a new ValarrayObj_var and grab space for a valarray object
         _val=new ValarrayObj<T>(start()+size());
         if (_val.is_nil())
             throw MatrixValueError();
    }
    
    void initialise()
    {
         _subarray=0;
         createval();
         _dimno=1;
         set_copyoption(DEEP);
    }        
    
    // Constructor below is called when indexing a Matrix object -
    // another Matrix object of a lower rank, which accesses the same set 
    // of data - or creating a "Lazy Matrix"
    Matrix(ValarrayObj_var v,Dim_var dim,\
        size_t idim=1,bool sub=false,MatrixCopyOption copyopt=DEEP,\
        bool not_empty=true):_diminfo(dim),_val(v),_dimno(idim),_subarray(sub),\
        _copyoption(copyopt),_not_empty(not_empty)
        {;}
    
    // cover function for "lazy" Matrix construction
    Matrix<T> matrix_result() const
    { return  Matrix<T>(data_objvar(),_diminfo,cdimno(),false,SHALLOW,\
          !empty());}
  
    Matrix<T>& clone(const Matrix<T>& m)
    // make current matrix a clone of m - i.e. points to same data
    {
       _val=m._val;
       _subarray=m._subarray;
       _dimno=m._dimno;
       _not_empty=m._not_empty;
       set_copyoption(DEEP); // futures copies of m will be deep copies
       _diminfo=m.dimvar();
       return *this;
    }

    // Force matrix object to be stored in contiguous memory
    void make_contiguous()
    {
    // force current matrix object to be stored in contiguous memory
        if (!dimvar()->is_contiguous())
        {
            // create a DEEP copy of matrix - will save in contiguous memory
            Matrix<T> m=deepcopy();
            // force current object to point to this new copy i.e. SHALLOW 
            // copy!!
            *this=m.shallowcopy();
        }
    }    
    
    // Matrix multiplication routines
    void matmult(const Matrix<T>& m,Matrix<T>& mres) const;
    
public:
    // --- Constructor Functions 
    Matrix():_diminfo(new Dim())
    { 
        #ifdef DIAG
        cout << "constructing rank 0 matrix" << endl;
        #endif
        initialise();
	T val;
        assign_val(val);
    }

    Matrix<T>(const T& val):_diminfo(new Dim())
    {
    // create a rank 0 Matrix with value val
        #ifdef DIAG
        cout << "constructing rank 0 matrix of " << val << endl;
        #endif
        initialise();
        assign_val(val);
    }
    
    template<class U>
    Matrix(const vector<U>& v):\
             _diminfo(new Dim(v.size()))
    {
    // create a rank 1 Matrix from vector
        initialise();
        assign_vector(v);
    } 
    
    Matrix(const Dim& d):_diminfo(new Dim(d.dims(),d.storageorder()))
    {
    // Dimension info provided
        #ifdef DIAG
        cout << "Constructing empty matrix..allocating space" << endl ;
        #endif
        initialise();
	T val;
        assign_val(val);
    }
    Matrix(const Dim& d,const T& val):\
       _diminfo(new Dim(d.dims(),d.storageorder()))
    {
    // Diminfo and fill value 
         #ifdef DIAG
         cout << "Constructing matrix of " << val << endl;
         #endif                
         initialise(); 
         assign_val(val);
    }  
       
    Matrix(const Dim& d,size_t n,const T* p):\
        _diminfo(new Dim(d.dims(),d.storageorder()))
    {
    // Diminfo and array of values
         #ifdef DIAG
         cout << "Constructing matrix from array" << endl;
         #endif
         initialise();
         assign_ptr(n,p); 
    }
    
    template<class U>
    Matrix(const Dim& d,const vector<U>& v):\
          _diminfo(new Dim(d.dims(),d.storageorder()))
    {
    // Diminfo and vector of values
        if (size()!=v.size())
            throw MatrixLengthError();
        initialise();
        assign_vector(v);
    }

    Matrix(const Dim& d,const valarray<T>& v):\
        _diminfo(new Dim(d.dims(),d.storageorder()))
    {
    // Diminfo and valarray of values
        initialise();
        assign_array(v);
    }
    
    Matrix(const Matrix<T>& m)
    {
      if (m.copyoption()==SHALLOW)
       _diminfo=m.dimvar();
      else 
       _diminfo=new Dim(m.dims());
      initialise();assign_matrix(m);
    }
 
    template<class U>
    Matrix(const Matrix<U>& m):_diminfo(new Dim(m.dims()))
    { initialise();assign_matrix(m);}
      
    Matrix(Dim_var d):_diminfo(d)
    { initialise();
      T val;
      assign_val(val);
    }
    
    class iterator
    {
        const Matrix<T>* _matrix;
        size_t _eltno; // element  no.
        
        T& ref() const 
        {
            if(!_matrix)
                throw MatrixValueError();
            valarray<size_t> v(_matrix->data_ptr()->size());
	    for(register int i=0;i<v.size();i++)
	      v[i]=i;
	    const gslice& g=_matrix->getgslice();
            valarray<size_t> vv(v[g]);
            return (*_matrix->data_ptr())[vv[_eltno%_matrix->size()]];
        }
        
    public:
        iterator(){;}
        iterator(const Matrix<T>* m,size_t eltno=0):_matrix(m),_eltno(eltno){;}
        operator T*(){ return &ref();}  
        iterator operator+(int i)
        {
            iterator itr(*this);
            itr._eltno+=i;
            return itr;
        }           
    
        iterator& operator++()
        {
             _eltno++;
             return *this;
        }
    
        iterator operator++(int)
        {
            iterator itr(*this);
            _eltno++;
            return itr;
        }

        T& operator*(){ return ref();}
        T* operator->(){return &ref();}
    
        bool operator==(const iterator& itr) const
        {
            if(_matrix!=itr._matrix)
                return false;
            return (_eltno==itr._eltno);
        }
    
        bool operator!=(const iterator& itr) const
        { return !(itr==*this);}
    
        bool operator<(const iterator& itr) const
        { 
            if(_matrix!=itr._matrix)
                throw MatrixDomainError();
            return (_eltno<itr._eltno);    
        }  
    
        size_t eltno() const {return _eltno;}        
       
    };
    
    iterator begin() const
    { return iterator(this);}
    
    iterator end() const
    { return iterator(this,this->size());} 
    
    class ListInitializer
    {
        iterator _itr;
    
    public:
        ListInitializer(const iterator&  itr):_itr(itr){;}
        
        ListInitializer operator,(const T& v)
        {
             *_itr=v;
             return ListInitializer(_itr+1);
        }
        
    };
    
    class ListInitialization
    {
        const Matrix<T>& _array;
        T _value;
        mutable bool _singleval;
        
    public:
        ListInitialization(const Matrix<T>& array,const T& value):\
            _array(array),_value(value),_singleval(true){;}
        ~ListInitialization()
        { 
           if(_singleval)
             _array.assign_val(_value);
        }
        
        ListInitializer operator,(const T& v)
        {
            _singleval=false;
            iterator itr=_array.begin();
            *itr=_value;
            itr++;
            *itr=v;
            itr++;
            return ListInitializer(itr);
        }
    };
    
    // --- Assign functions
    void assign_val(const T& t) const
    // sets Matrix values to t
    {
        valarray<T> v(t,size());
        assign_array(v);
    }            
    
    void assign_array(const valarray<T>& v) const
    // set Matrix values to v
    {
        if (_val.is_nil())
            throw MatrixValueError();
        if (size()!=v.size())
            throw MatrixLengthError();
        (*_val)[getgslice()]=v;
        _not_empty=(v.size()!=0);
    }
    
    template<class U>
    void assign_array(const valarray<U>& v) const
    // set Matrix values to v
    {
        if (_val.is_nil())
            throw MatrixValueError();
        if (size()!=v.size())
            throw MatrixLengthError();
	
	valarray<T> v_data(size());
	for(register int i=0;i<size();i++)
	  v_data[i]=v[i];
	_val.ptr_()->operator[](getgslice())=v_data;

        _not_empty=(v.size()!=0);
    }

    void assign_ptr(size_t n,const T* p) const
    // sets values in array pointed to by p to Matrix object
    {
        if(p==NULL)
            throw MatrixValueError();
        if(n!=size())
            throw MatrixLengthError();
        valarray<T> v(p,size());
        assign_array(v);
    }
     
    template<class U>
    void assign_vector(const vector<U>& v) const
    {
    // Assign a vector<U> to a Matrix 
        if (size()!=v.size())
            throw MatrixLengthError();
	valarray<T> v_data(size());
	for(register int i=0;i<size();i++)
	  v_data[i]=v[i];
	_val.ptr_()->operator[](getgslice())=v_data;
        _not_empty=(v.size()!=0);
    }
    
    void copy(const Matrix<T>& m) const
    // Copy data from Matrix m to current Matrix object
    {
        if(this!=&m) // different objects
        {
            // check that there is data to copy
            if (m._val.is_nil())
                throw MatrixValueError();
	    const gslice& g(m.getgslice());
	    valarray<T> v((size_t)product(g.size()));
	    v=(*m.data_ptr())[g];
            assign_array(v);
        }
    }

    template<class U>
    void copy(const Matrix<U>& m) const
    // Copy data from Matrix m to current Matrix object
    {
        // check that there is data to copy
        if (m.data_objvar().is_nil())
            throw MatrixValueError();
        if(size()!=m.size())
           throw MatrixLengthError();
	const gslice& g(m.getgslice());
	valarray<U> v((size_t)product(g.size()));
	v=(*m.data_ptr())[g];
        assign_array(v);
    }
    
    void assign_matrix(const Matrix<T>& m)
    {
    // force a deepcopy if copying into a subarray           
        if (m.copyoption()==SHALLOW && !_subarray)
        {
            #ifdef DIAG
            cout << "Copy assignment - Shallow copy" << endl;
            #endif
            // cannot perform a shallow copy into a subarray
            if (this!=&m) // test for self-assignment              
                 clone(m);
        }
        else 
        { 
            #ifdef DIAG
            cout << "Copy assignment - Deep copy" << endl;
            #endif
            if(m.rank()==0)
                assign_val(m.scalar());
            else
            {               
                if (!samesize(m))
                    resize(m.dimdata());
                                
                if(this==&m)  // self assignment - force a deep copy
                {
                    Matrix<T> tmp(dimdata());
                    tmp.copy(m);
                    clone(tmp);
                }
                else     
                    copy(m);
            }
        }
    }
    
    template<class U>
    void assign_matrix(const Matrix<U>& m)
    {
    // --- Assignment of a Matrix object to another Matrix object
        #ifdef DIAG
        cout << "Copy assignment - Deep copy" << endl;
        #endif
        if(m.rank()==0)
        {
            T t=(T)m.scalar();
            assign_val(t);
        }
        else
        {
            if (!samesize(m))
                resize(m.dimdata());              
            copy(m);
        }
    }
    
    template<class U>
    void assign_matrix(const Matrix<U>& m) const
    {
        #ifdef DIAG
        cout << "Copy assignment - Deep copy" << endl;
        #endif
        copy(m);    
    }
        
    // --- Member Conversion functions
    T scalar() const
    { 
    //converts a rank 0 Matrix to a T
        #ifdef DIAG
        cout << "Scalar conversion" << endl;
        #endif
        if (rank()!=0)
            throw MatrixRankError();
        T result;
        result=(*_val)[start()];
        return result;
    }
    
    operator T() const
    {return scalar();} 
    
    operator vector<T>() const
    {
    // convert a rank 1 Matrix to a vector.. should we allow this for any
    // rank??
        #ifdef DIAG
        cout << "Vector conversion" << endl;
        #endif
        if (rank()!=1)
            throw MatrixRankError();       
        vector<T> v(size());
	valarray<T> v_data(_val.ptr_()->operator[](getgslice()));
	for(register int i=0;i<size();i++)
	  v[i]=v_data[i];
        return v;
    }
    
    T* ptr() 
    {
    // Convert a Matrix of any rank into a T* - NB care must be taken
    // with memory allocation i.e. the resultant T* will need to be
    // deleted in the calling function if copy_to_ptr used
        #ifdef DIAG
        cout << "Conversion to a T*" << endl;
        #endif
        make_contiguous();
        return begin();
        //return copy_to_ptr();
    }                  
    
    T* copy_to_ptr() const
    {
    // copy all elements of a Matrix out to a T* i.e. an array of  T   
        //T* parr=new T[size()];
        size_t sz=_val->size(); // size of valarray
        valarray<size_t> ind(sz);
        iota(&ind[0],&ind[sz],0);    // generate a valarray of indices
        valarray<size_t> sub_ind(ind[getgslice()]); // get gslice indices
        sort(&sub_ind[0],&sub_ind[size()]); // sort these indices in storage order
        valarray<T>* val=new valarray<T>((*_val)[sub_ind]);
        return &((*val)[0]);
        //for(int i=0;i<size();i++)
        //    *(parr+i)=(*_val)[sub_ind[i]];   // output to T*
    }   
    
    //operator T*() const
    //{ return ptr();}  

    // --- Operator = (Copy assignment)
    ListInitialization operator=(const T& t)
    { return ListInitialization(*this,t);}
    ListInitialization operator=(const T& t) const
    { return ListInitialization(*this,t);}  
    template<class U>
    Matrix<T>& operator=(const vector<U>& v)
    { assign_vector(v); return *this;}
    template<class U>
    const Matrix<T>& operator=(const vector<U>& v) const
    { assign_vector(v); return *this;}
    Matrix<T>& operator=(const Matrix<T>& m)
    { assign_matrix(m); return *this;}
    const Matrix<T>& operator=(const Matrix<T>& m) const
    { assign_matrix(m); return *this;}
    template<class U>
    Matrix<T>& operator=(const Matrix<U>& m)
    { assign_matrix(m); return *this;}
    template<class U>
    const Matrix<T>& operator=(const Matrix<U>& m) const
    { assign_matrix(m); return *this;}
    
    // --- Operators for indexing []
    Matrix<T> operator[](const DummyArg& arg) const;
    Matrix<T> operator[](const size_t i) const;
    Matrix<T> operator[](const slice& sl) const;
    Matrix<T> operator[](const gslice& g) const;
    
    // --- Operators for () indexing
    // --- Rank 1
    T& operator()(const size_t i) const;
    Matrix<T> operator()(const slice& sl) const;
    // --- Rank 2
    T& operator()(const size_t i,const size_t j) const;
    Matrix<T> operator()(const slice& sl1,const size_t j) const;
    Matrix<T> operator()(const size_t i,const slice& sl2) const;
    Matrix<T> operator()(const slice& sl1,const slice& sl2) const;
    Matrix<T> operator()(const DummyArg& arg,const size_t j) const;
    Matrix<T> operator()(const DummyArg& arg,const slice& sl) const;
    Matrix<T> operator()(const size_t i,const DummyArg& arg) const;
    Matrix<T> operator()(const slice& sl,const DummyArg& arg) const;
    // --- Rank 3
    T& operator()(const size_t i,const size_t j,const size_t k) const;
    Matrix<T> operator()(const slice& sl1,const size_t j,const size_t k) const;
    Matrix<T> operator()(const size_t i,const slice& sl2,const size_t k) const;
    Matrix<T> operator()(const size_t i,const size_t j,const slice& sl3) const;
    Matrix<T> operator()(const slice& sl1,const slice& sl2,const size_t k) const;
    Matrix<T> operator()(const slice& sl1,const size_t j,const slice& sl3) const;
    Matrix<T> operator()(const size_t i,const slice& sl1,const slice& sl2) const;
    Matrix<T> operator()(const slice& sl1,const slice& sl2,const slice& sl3) const;
    Matrix<T> operator()(const DummyArg& arg,const size_t j,const size_t k)\
      const;
    Matrix<T> operator()(const DummyArg& arg,const slice& sl,const size_t k)\
      const;
    Matrix<T> operator()(const DummyArg& arg,const size_t j,const slice& sl)\
      const;
    Matrix<T> operator()(const DummyArg& arg,const slice& sl1,const slice& sl2)\
      const;      
    Matrix<T> operator()(const size_t i,const DummyArg& arg,const size_t k)\
      const;
    Matrix<T> operator()(const slice& sl,const DummyArg& arg,const size_t k)\
      const;
    Matrix<T> operator()(const size_t i,const DummyArg& arg,const slice& sl)\
      const;  
    Matrix<T> operator()(const slice& sl1,const DummyArg& arg,const slice& sl2)\
      const;            
    Matrix<T> operator()(const size_t i,const size_t j,const DummyArg& arg)\
      const;
    Matrix<T> operator()(const slice& sl,const size_t j,const DummyArg& arg)\
      const;
    Matrix<T> operator()(const size_t i,const slice& sl,const DummyArg& arg)\
      const;
    Matrix<T> operator()(const slice& sl1,const slice& sl2,const DummyArg& arg)\
      const;
    Matrix<T> operator()(const DummyArg& arg1,const DummyArg& arg2,\
                           const size_t k) const;
    Matrix<T> operator()(const DummyArg& arg1,const DummyArg& arg2,\
                           const slice& sl) const;
    Matrix<T> operator()(const DummyArg& arg1,const size_t j,\
                           const DummyArg& arg2) const;
    Matrix<T> operator()(const DummyArg& arg1,const slice& sl,\
                           const DummyArg& arg2) const;                           
    Matrix<T> operator()(const size_t i,const DummyArg& arg1,\
                           const DummyArg& arg2) const;
    Matrix<T> operator()(const slice& sl,const DummyArg& arg1,\
                           const DummyArg& arg2) const;
    Matrix<T> operator()(const gslice& g) const;
    
    // --- Operator for negation
    MatrixRes<T> operator-() const;
    // --- Operator for identity
    MatrixRes<T> operator+() const; 
    // --- Operator for transpose
    MatrixRes<T> operator~() const;  
    
    // --- Operator for power 
    MatrixRes<typename promote_trait<T,double>::T_promote> operator^(const double& n) const;
    template<class U>
    MatrixResSeq<typename promote_trait<T,double>::T_promote> \
        operator^(const vector<U>& v) const;
    
    // --- Operator for addition and subtraction
    MatrixRes<T> operator+(const T& t) const;  // Scalar addition   
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator+(const Scalar<U>& s) const;  // Scalar addition   
    MatrixRes<T> operator+(const Matrix<T>& m) const; // Matrix addition
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator+(const Matrix<U>& m) const;
   
    MatrixRes<T> operator-(const T& t) const; // Scalar subtraction 
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator-(const Scalar<U>& s) const; // Scalar subtraction 
    MatrixRes<T> operator-(const Matrix<T>& m) const; // Matrix subtraction
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote > operator-(const Matrix<U>& m) const;
   
    // --- Operators for computed +- assignments
    Matrix<T>& operator+=(const T& t);
    template<class U>
    Matrix<T>& operator+=(const Scalar<U>& s);
    const Matrix<T>& operator+=(const T& t) const;
    template<class U>
    const Matrix<T>& operator+=(const Scalar<U>& s) const;
    Matrix<T>& operator+=(const Matrix<T>& m);
    const Matrix<T>& operator+=(const Matrix<T>& m) const;
    template<class U>
    Matrix<T>& operator+=(const Matrix<U>& m);
    template<class U>
    const Matrix<T>& operator+=(const Matrix<U>& m) const;
   
    Matrix<T>& operator-=(const T& t);
    template<class U>
    Matrix<T>& operator-=(const Scalar<U>& s);
    const Matrix<T>& operator-=(const T& t) const;
    template<class U>
    const Matrix<T>& operator-=(const Scalar<U>& s) const;
    Matrix<T>& operator-=(const Matrix<T>& m);
    const Matrix<T>& operator-=(const Matrix<T>& m) const;
    template<class U>
    Matrix<T>& operator-=(const Matrix<U>& m);
    template<class U>
    const Matrix<T>& operator-=(const Matrix<U>& m) const;
  
    // --- Operators for multiplication and division
    MatrixRes<T> operator*(const T& t) const;
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator*(const Scalar<U>& s) const;
    MatrixRes<T> operator*(const Matrix<T>& m) const;
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator*(const Matrix<U>& m) const;
      
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator/(const Scalar<U>& t) const;
    MatrixRes<T> operator/(const Matrix<T>& m) const;
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator/(const Matrix<U>& m) const;
    
    //--- Computed scalar multiplication/division
    Matrix<T>& operator*=(const T& t);
    template<class U>
    Matrix<T>& operator*=(const Scalar<U>& s);
    const Matrix<T>& operator*=(const T& t) const;
    template<class U>
    const Matrix<T>& operator*=(const Scalar<U>& s) const;
    Matrix<T>& operator*=(const Matrix<T>& t);
    template<class U>
    Matrix<T>& operator*=(const Matrix<U>& m);
   
    Matrix<T>& operator/=(const T& t);
    template<class U>
    Matrix<T>& operator/=(const Scalar<U>& s);
    const Matrix<T>& operator/=(const T& t) const;
    template<class U>
    const Matrix<T>& operator/=(const Scalar<U>& s) const;
    Matrix<T>& operator/=(const Matrix<T>& t);
    const Matrix<T>& operator/=(const Matrix<T>& t) const;
    template<class U>
    Matrix<T>& operator/=(const Matrix<U>& m);
    template<class U>
    const Matrix<T>& operator/=(const Matrix<U>& m) const;    
    
    // --- Logical operators
    bool operator==(const Matrix<T>& m) const;
    bool operator!=(const Matrix<T>& m) const;
     
    // ---  Misc functions    
    size_t cdimno() const{return _dimno;}
    void changerank(const size_t rk);
    Matrix<T> col(const size_t i) const;  // access a col
    Matrix<T> col(const slice& sl) const;          
    MatrixCopyOption copyoption() const
    {return _copyoption;}
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> cross(const Matrix<U>& m) const;    
    valarray<T>& data(){return *_val;}
    const valarray<T>& data() const{return *_val;}
    ValarrayObj_var data_objvar() const{return _val;} // ValarrayObj_var        
    valarray<T>* data_ptr() const{return _val.ptr_();} // ptr to valarray    
    Matrix<T>& deepcopy(){set_copyoption(DEEP);return *this;}        
    Matrix<T> diagonal() const;
    size_t& dimno(){return _dimno;}
    Void display() const; // display matrix object    
    MatrixRes<T> duplicate(MatrixCopyOption option=DEEP) const;
    bool empty() const { return !_not_empty;}    
    MatrixRes<T> fortranarray() const; 
    Matrix<T> indexdim(const size_t idim,const size_t i) const; // access any
    Matrix<T> indexdim(const size_t idim,const slice& sl) const;// dim
    template<class U>  
    typename promote_trait<T,U>::T_promote dot(const Matrix<U>& m) const;   
    Matrix<T> layer(const size_t i) const; // access third dim
    Matrix<T> layer(const slice& sl) const;     
    T& matrix_indexelt(size_t *ind,size_t nel) const;
    T& matrix_indexelt(valarray<size_t> v) const;
    T& matrix_indexelt(const vector<int>& v) const;
    MatrixResSeq<T> matrix_to_seq(size_t dimno=1) const;    
    MatrixRes<T> mult_element(const Matrix<T>& m) const;
    void null_matrix() const;
    MatrixRes<typename promote_trait<T,double>::T_promote> pinverse() const;
    MatrixRes<T> plusreduce(const size_t idim=1) const;
    void plusreduce(const size_t idim,Matrix<T>& mres) const;    
    MatrixRes<typename promote_trait<T,double>::T_promote> power(const double& n) const;    
    void reshape(const Dim& d); 
    void resize(const Dim& d); 
    Matrix<T> row(const size_t i) const;  // access a row
    Matrix<T> row(const slice& sl) const;    
    void set_copyoption(MatrixCopyOption copyopt=DEEP)
    {_copyoption=copyopt;}      
    void set_storageorder(DimStorageOption order); 
    Matrix<T>& shallowcopy(){set_copyoption(SHALLOW);return *this;}
    MatrixRes<typename promote_trait<T,double>::T_promote> solve(const Matrix<T>& B) const;           
    void storage(DimStorageOption order=ROW_MAJOR_ORDER);
    bool subarray() const{return _subarray;}
    T sum() const {return (*_val).sum();} // return sum of all elts       
    T trace() const;
    void transpose(); 
    void unit_matrix() const;

    // --- Cover functions for Dim Object
    virtual Dim_var dimvar() const{return _diminfo;}
    Dim& dimdata() const{return *dimvar();}
    valarray<size_t> dims() const{return dimvar()->dims();} // dimension info
    gslice& getgslice() const{return dimvar()->getgslice();}  // return gslice of info
    //void get_dims(valarray<size_t>& dms) const{dimvar()->get_dims(dms);}
    void get_dims(vector<int>& dms) const {dimvar()->get_dims(dms);}
    void get_strides(vector<int>& str) const{dimvar()->get_strides(str);}
    bool is_symmetric() const;   
    size_t rank() const{return dimvar()->rank();}  // no. of dimensions 
    template<class U>
    bool samesize(const Matrix<U>& m) const
    { return dimdata().samesize(m.dimdata());}           
    size_t size() const{return dimvar()->size();}  // total no. elts
    size_t start() const{return dimvar()->start();} // start index 
    DimStorageOption storageorder() const{return dimvar()->storageorder();}         
    valarray<size_t> strides() const{return dimvar()->strides();} // strides    

};
// ----------- End of class Matrix -------------------
template<class T> Vostream& operator<< (Vostream& v_os, \
  const Matrix<T>& m)
{  const Dim& d=m.dimdata();
   v_os << d;
   const valarray<T>& v=m.data();
   v_os << v;
   return v_os;
}

template<class T> Vistream& operator>> (Vistream& v_is,Matrix<T>& m)
{ 
  Dim d;
  v_is >> d;
  valarray<T> v;
  v_is >> v;
  m=Matrix<T>(d,v);
  return v_is;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixRes


// Derived class of Matrix: class for returning QNmatrices by value from
// function calls ...avoids unecessary deep copying
// MatrixRes objects are Matrices with flag set for shallow copying..
// so result of a matrix expression is shallow copied into matrix object
// on assignment/construction
// Avoids data being copied around unecessarily when temporaries are created
// Utility class only ...transparent to user
//------------------------------------------------------------------------------
template<class T>
class MatrixRes : public Matrix<T>
{
 public:
     MatrixRes():Matrix<T>(){;}
     MatrixRes(const Matrix<T>& m):Matrix<T>(m.data_objvar(), m._diminfo,\
         m.cdimno(),m._subarray,SHALLOW, !m.empty())
     {;}
 }; 
// ----------- End of class MatrixRes -------------------
template<class T>
Vostream& operator<<(Vostream&,const Scalar<T>&);
template<class T>
Vistream& operator>>(Vistream&,Scalar<T>&);
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class Scalar


template<class T>
class Scalar
{
  T _val;
public:
    Scalar():_val(0.){;}
    Scalar(const Matrix<T>& mat):_val(mat.scalar()){;}
    template<class U>
    Scalar(const Matrix<U>& mat):_val((T)mat.scalar()){;}
    Scalar(const T& tval):_val(tval){;}
    Scalar(const Scalar<T>& sca):_val(sca.scalar()){;}
    
    Scalar& operator=(const Matrix<T>& mat)
    {
      if(mat.rank()!=0)
          throw MatrixRankError();  
      _val=mat.scalar(); return *this;} 
      
    template<class U>
    Scalar& operator=(const Matrix<U>& mat)
    {
      if(mat.rank()!=0)
          throw MatrixRankError();  
      _val=mat.scalar(); return *this;} 
      
    Scalar& operator=(const T& tval)
    { _val=tval; return *this;}
    
    operator T() const
    {return scalar();} 

    T scalar() const
    { return _val;}
    
    size_t rank() const
    { return 0;}

    Dim dimdata() const
    { return Dim();}
    
    DimStorageOption storageorder() const{return ROW_MAJOR_ORDER;}         
    // --- Operator for power
    Scalar<typename promote_trait<T,double>::T_promote> operator^(const double& n) const;
    template<class U>
    ScalarSeq<typename promote_trait<T,double>::T_promote> \
        operator^(const vector<U>& v) const;
    
    // --- Operator for negation
    Scalar<T> operator-() const
    { return Scalar<T>(-_val);}

    // --- Operator for identity
    Scalar<T> operator+() const 
    { return Scalar<T>(_val);}
    
    // --- Operator for addition and subtraction
    Scalar<T> operator+(const T& t) const
    { return Scalar<T>(_val+t);}  // Scalar addition   
    
    Scalar<T> operator+(const Scalar<T>& s) const
    { return Scalar<T>(_val+s.scalar());}  // Scalar addition 
   
    template<class U>
    Scalar<typename promote_trait<T,U>::T_promote> operator+(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return Scalar<T_promote>(val1+val2);
    } 
  
    MatrixRes<T> operator+(const Matrix<T>& m) const
    { return m+_val;} // Matrix addition
    
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator+(const Matrix<U>& m) const
    { return add_scalar(m,_val);}

    Scalar<T> operator-(const T& t) const
    { return Scalar<T>(_val-t);} // Scalar subtraction 
    
    Scalar<T> operator-(const Scalar<T>& s) const
    { return Scalar<T>(_val-s.scalar());}  // Scalar subtraction 
   
    template<class U>
    Scalar<typename promote_trait<T,U>::T_promote> operator-(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return Scalar<T_promote>(val1-val2);
    } 

    MatrixRes<T> operator-(const Matrix<T>& m) const
    { return (-m)+_val;} // Matrix subtraction
    
    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote > operator-(const Matrix<U>& m) const
    { return sub_scalar(m,_val);}
   
    // --- Operators for computed +- assignments
    Scalar<T>& operator+=(const T& t)
    { _val+=t;return *this;}
    
    Scalar<T>& operator+=(const Scalar<T>& s) 
    { _val+=s.scalar(); return *this;} 
   
    template<class U>
    Scalar<T>& operator+=(const Scalar<U>& s) 
    {  _val+=(T)s.scalar(); return *this;}

    Scalar<T>& operator+=(const Matrix<T>& m)
    { if(m.rank()!=0)
        throw MatrixRankError();
      _val+=m.scalar();
      return *this;
    }

    template<class U>
    Scalar<T>& operator+=(const Matrix<U>& m)
    { if(m,rank()!=0)
        throw MatrixRankError();
      _val+=(T)m.scalar();
      return *this;
    }
   
    Scalar<T>& operator-=(const T& t)
    { _val-=t; return *this;}
    
    Scalar<T>& operator-=(const Scalar<T>& s) 
    { _val-=s.scalar(); return *this;} 
   
    template<class U>
    Scalar<T>& operator-=(const Scalar<U>& s) 
    {  _val-=(T)s.scalar(); return *this;}

    Scalar<T>& operator-=(const Matrix<T>& m)
    { if(m.rank()!=0)
        throw MatrixRankError();
      _val-=m.scalar();
      return *this;
    }

    template<class U>
    Scalar<T>& operator-=(const Matrix<U>& m)
    { if(m.rank()!=0)
        throw MatrixRankError();
      _val-=m.scalar();
      return *this;
    }
  
    // --- Operators for multiplication and division
    Scalar<T> operator*(const T& t) const
    {return Scalar<T>(_val*t);}
    
    Scalar<T> operator*(const Scalar<T>& s) const
    { return Scalar<T>(_val*s.scalar());}  // Scalar addition 
   
    template<class U>
    Scalar<typename promote_trait<T,U>::T_promote> operator*(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return Scalar<T_promote>(val1*val2);
    } 

    MatrixRes<T> operator*(const Matrix<T>& m) const
    { return m*_val;}

    template<class U>
    MatrixRes<typename promote_trait<T,U>::T_promote> operator*(const Matrix<U>& m) const
    {return mult_scalar(m,_val);}
      
    Scalar<T> operator/(const Scalar<T>& s) const
    { return Scalar<T>(_val/s.scalar());}  // Scalar addition 
   
    template<class U>
    Scalar<typename promote_trait<T,U>::T_promote> operator/(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return Scalar<T_promote>(val1/val2);
    } 

    Scalar<T> operator/(const Matrix<T>& m) const
    { return Scalar<T>(_val/m.scalar());}

    template<class U>
    Scalar<typename promote_trait<T,U>::T_promote> operator/(const Matrix<U>& m) const
    {
       Scalar<typename promote_trait<T,U>::T_promote> s1(_val);
       Scalar<typename promote_trait<T,U>::T_promote> s2(m.scalar());
       return Scalar<typename promote_trait<T,U>::T_promote>(s1.scalar()/s2.scalar());
    }

    //--- Computed scalar multiplication/division
    Scalar<T>& operator*=(const T& t)
    { _val*=t; return *this;}
    
    Scalar<T>& operator*=(const Scalar<T>& s) 
    { _val*=s.scalar(); return *this;} 
   
    template<class U>
    Scalar<T>& operator*=(const Scalar<U>& s) 
    {  _val*=(T)s.scalar(); return *this;}
     
    Scalar<T>& operator*=(const Matrix<T>& m)
    { if(m.rank()!=0)
        throw MatrixRankError();
      _val*=m.scalar(); return *this;
    }
    
    template<class U>
    Scalar<T>& operator*=(const Matrix<U>& m)
    { if(m.rank()!=0)
        throw MatrixRankError();
      _val*=(T)m.scalar(); return *this;
    }
   
    Scalar<T>& operator/=(const T& t)
    { _val/=t; return *this;}

    Scalar<T>& operator/=(const Scalar<T>& s) 
    { _val/=s.scalar(); return *this;} 
   
    template<class U>
    Scalar<T>& operator/=(const Scalar<U>& s) 
    {  _val/=(T)s.scalar(); return *this;}

    Scalar<T>& operator/=(const Matrix<T>& m)
    { if(m.rank()!=0)
        throw MatrixRankError();
      _val/=m.scalar(); return *this;
    }
    template<class U>
    Scalar<T>& operator/=(const Matrix<U>& m)
    { if(m.rank()!=0)
        throw MatrixRankError();
      _val/=(T)m.scalar(); return *this;
    }
    // -- relational operators 
    bool operator>(const T& t) const
    { return (_val>t);}  
    
    bool operator>(const Scalar<T>& s) const
    { return (_val>s.scalar());}  // Scalar subtraction 
   
    template<class U>
    bool operator>(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return (val1>val2);
    }
    
    bool operator>=(const T& t) const
    { return (_val>=t);}  
    
    bool operator>=(const Scalar<T>& s) const
    { return (_val>=s.scalar());}  // Scalar subtraction 
   
    template<class U>
    bool operator>=(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return (val1>=val2);
    }

    bool operator<(const T& t) const
    { return (_val<t);}  
    
    bool operator<(const Scalar<T>& s) const
    { return (_val<s.scalar());}  // Scalar subtraction 
   
    template<class U>
    bool operator<(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return (val1<val2);
    }
    
    bool operator<=(const T& t) const
    { return (_val<=t);}  
    
    bool operator<=(const Scalar<T>& s) const
    { return (_val<=s.scalar());}  // Scalar subtraction 
    
    template<class U>
    bool operator<=(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return (val1<=val2);
    }
    
    bool operator==(const T& t) const
    { return (_val==t);}  
    
    bool operator==(const Scalar<T>& s) const
    { return (_val==s.scalar());}  // Scalar subtraction 
   
    template<class U>
    bool operator==(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return (val1==val2);
    }
    
    bool operator!=(const T& t) const
    { return (_val!=t);}  
    
    bool operator!=(const Scalar<T>& s) const
    { return (_val!=s.scalar());}  // Scalar subtraction 
   
    template<class U>
    bool operator!=(const Scalar<U>& s) const
    { typedef typename promote_trait<T,U>::T_promote T_promote;
      T_promote val1=_val;
      T_promote val2=s.scalar();
      return (val1>=val2);
    }
   
    Scalar<typename promote_trait<T,double>::T_promote> pinverse() const;
    
};

// ----------- End of class Scalar -------------------
template<class T>
Vostream& operator<<(Vostream& v_os,const Scalar<T>& scalar)
{  T val=scalar.scalar(); v_os << val; return v_os;}

template<class T>
Vistream& operator>>(Vistream& v_is,Scalar<T>& scalar)
{ T val; v_is >> val; scalar=val; return v_is;}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq
     
// Base Class for a sequence of Matrix Containers i.e. Matrix<T>, Scalar<T>,
// and Matrixres<T>
//------------------------------------------------------------------------------
template<class T, class ContainerT > 
class MatrixBaseSeq : public NumericSequence<ContainerT > 
{
           
 protected:
      // --- constructors
      MatrixBaseSeq(int num=0):\
          NumericSequence<ContainerT >(num,ContainerT()){;}
      MatrixBaseSeq(int num,const T& t):\
          NumericSequence<ContainerT >(num,ContainerT(t)){;}
      MatrixBaseSeq(int num,const ContainerT& m):\
        NumericSequence<ContainerT >(num,m){;}
      MatrixBaseSeq(const MatrixBaseSeq<T, ContainerT >& ms):\
          NumericSequence<ContainerT >(ms){;}
      template<class U>
      MatrixBaseSeq(const vector<U>& v):NumericSequence<ContainerT >()
      {assign(v);}

 public:     
      // --- computed assignments
      MatrixBaseSeq<T, ContainerT >& operator+=(const T& tval)
      { add_assign(tval);return *this;}
      MatrixBaseSeq<T, ContainerT >& operator+=(const vector<T>& v)
      { add_assign(v); return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator+=(const Scalar<U>& s)
      {add_assign(s); return*this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator+=(const Matrix<U>& m)
      {add_assign(m); return*this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator+=(const vector<Matrix<U> >& ms)
      { add_assign(ms);return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator+=(const vector<Scalar<U> >& ms)
      { add_assign(ms);return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator+=(const vector<MatrixRes<U> >& ms)
      { add_assign(ms);return *this;}
    
      MatrixBaseSeq<T, ContainerT >& operator-=(const T& tval)
      { sub_assign(tval);return *this;}
      MatrixBaseSeq<T, ContainerT >& operator-=(const vector<T>& v)
      { sub_assign(v); return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator-=(const Scalar<U>& s)
      {sub_assign(s); return*this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator-=(const Matrix<U>& m)
      {sub_assign(m); return*this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator-=(const vector<Matrix<U> >& ms)
      { sub_assign(ms);return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator-=(const vector<Scalar<U> >& ms)
      { sub_assign(ms);return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator-=(const vector<MatrixRes<U> >& ms)
      { sub_assign(ms);return *this;}
     
      MatrixBaseSeq<T, ContainerT >& operator*=(const T& tval)
      { mul_assign(tval);return *this;}
      MatrixBaseSeq<T, ContainerT >& operator*=(const vector<T>& v)
      { mul_assign(v); return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator*=(const Scalar<U>& s)
      {mul_assign(s); return*this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator*=(const Matrix<U>& m)
      {mul_assign(m); return*this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator*=(const vector<Matrix<U> >& ms)
      { mul_assign(ms);return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator*=(const vector<Scalar<U> >& ms)
      { mul_assign(ms);return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator*=(const vector<MatrixRes<U> >& ms)
      { mul_assign(ms);return *this;}
      
      MatrixBaseSeq<T, ContainerT >& operator/=(const T& tval)
      { div_assign(tval);return *this;}
      MatrixBaseSeq<T, ContainerT >& operator/=(const vector<T>& v)
      { div_assign(v); return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator/=(const Scalar<U>& s)
      {div_assign(s); return*this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator/=(const Matrix<U>& m)
      {div_assign(m); return*this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator/=(const vector<Matrix<U> >& ms)
      { div_assign(ms);return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator/=(const vector<Scalar<U> >& ms)
      { div_assign(ms);return *this;}
      template<class U>
      MatrixBaseSeq<T, ContainerT >& operator/=(const vector<MatrixRes<U> >& ms)
      { div_assign(ms);return *this;}
      
      // --- operator+
      MatrixResSeq<T> operator+(const T& t) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator+(const Scalar<U>& s) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator+(const Matrix<U>& m) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator+\
          (const vector<Matrix<U> >& ms) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator+\
          (const vector<Scalar<U> >& ms) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator+\
          (const vector<MatrixRes<U> >& ms) const;          
          
      // --- operator-
      MatrixResSeq<T> operator-() const;
      MatrixResSeq<T> operator-(const T& t) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator-(const Scalar<U>& s) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator-(const Matrix<U>& m) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator-\
          (const vector<Matrix<U> >& ms) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator-\
          (const vector<Scalar<U> >& ms) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator-\
          (const vector<MatrixRes<U> >& ms) const;
 
      // --- operator*
      MatrixResSeq<T> operator*(const T& t) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator*(const Scalar<U>& s) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator*(const Matrix<U>& m) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator*\
          (const vector<Matrix<U> >& ms) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator*\
          (const vector<Scalar<U> >& ms) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator*\
          (const vector<MatrixRes<U> >& ms) const;          
          
      // --- operator/
      MatrixResSeq<T> operator/(const T& t) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator/(const Scalar<U>& s) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator/(const Matrix<U>& m) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator/\
          (const vector<Matrix<U> >& ms) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator/\
          (const vector<Scalar<U> >& ms) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,U>::T_promote> operator/\
          (const vector<MatrixRes<U> >& ms) const;
          
       // --- operator^
      MatrixResSeq<typename promote_trait<T,double>::T_promote> operator^(\
          const double& t) const;
      template<class U>
      MatrixResSeq<typename promote_trait<T,double>::T_promote> operator^(\
          const vector<U>& v) const;
                
      MatrixRes<T> sum() const;
      MatrixResSeq<T> plus_scan() const;
      void display() const; 
      Matrix<T> seq_to_matrix(size_t dimno=1) const;
      bool is_homogeneous_seq() const;
      virtual bool is_scalar_seq() const;
      MatrixResSeq<T> get_element_seq(const vector<int>& v) const;
      MatrixResSeq<T> get_element_seq(valarray<size_t> ind) const;
      MatrixResSeq<T> get_submatrix_seq(const gslice& g) const;
};
// ----------- End of class MatrixBaseSeq -------------------

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixSeq
     
// Class for a sequence of Matrices
//------------------------------------------------------------------------------
template<class T>
class MatrixSeq : public MatrixBaseSeq<T, Matrix<T> > 
{
public:
    // do we need to replace Matrix<U> with U below? Probably?
    // --- constructors
    MatrixSeq(int num=0):MatrixBaseSeq<T, Matrix<T> >()
    {
		Dim_var d=new Dim();
        Matrix<T> m(d);
        m.set_copyoption(DEEP);
		this->assign(num, m);
    }
    MatrixSeq(int num,const Dim& d):MatrixBaseSeq<T,Matrix<T> >()
    {
	  Dim_var dv=new Dim(d);
	  Matrix<T> m(dv); 
	  m.set_copyoption(DEEP);
      this->assign(num, m);
    }
    template<class U>
    MatrixSeq(int num,const Matrix<U>& m):MatrixBaseSeq<T, Matrix<T> >()
    {
		Dim_var d=new Dim(m.dims());

        Matrix<T> m2(d);
        m2.assign_matrix(m);
        m2.set_copyoption(DEEP);
		this->assign(num, m2);
    }
    MatrixSeq(const MatrixSeq<T>& ms):MatrixBaseSeq<T, Matrix<T> >()
    {
		int sz=ms.size();
		if(sz!=0)
		{  
			Dim_var d=new Dim(ms.elt(0).dims());
			Matrix<T> m(d);
			m.assign_matrix(ms[0]);
			m.set_copyoption(SHALLOW);
			this->assign(sz, m);
			for(register int i=1;i<sz;i++)
			{
				Matrix<T> m(d);
				m.assign_matrix(ms[i]);
				m.set_copyoption(SHALLOW);
				this->at(i) = m;
			}
		}
    }

    template<class U>
    MatrixSeq(const vector<U>& v):MatrixBaseSeq<T, Matrix<T> >()
    {
       int sz=v.size();
       if(sz!=0)
       {
          Matrix<T> m0(v[0]);
          m0.set_copyoption(SHALLOW);
		  this->assign(sz, m0);
          Dim_var d=this->elt(0).dimvar();
          for(register int i=1;i<sz;i++)
          {
            this->at(i) = v[i];
		  }
       }
    }

    // --- assignments
    MatrixSeq& operator=(const MatrixSeq<T>& seq)
    { assign(seq); return *this;}
    template<class U>
    MatrixSeq& operator=(const vector<U>& v)
    { assign(v); return *this;}
    template<class U>
    MatrixSeq& operator=(const Matrix<U>& s)
    { assign(this->size(),s); return *this;}
    MatrixSeq& operator=(const T& t)
    { assign(this->size(),t); return *this;}
    
    MatrixSubSeq<T> operator[](const DummyArg& arg) const;
    Matrix<T>& operator[](int i) 
    { return this->elt(i);}
    const Matrix<T>& operator[](int i) const
    { return this->elt(i);}
        
         
};
// ----------- End of class MatrixSeq -------------------

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixResSeq

// Class for returning sequences of Matrices by value from function calls ...
// avoids unecessary deep copying
// Utility class only. Transparent to users
//------------------------------------------------------------------------------
template<class T>
class MatrixResSeq:public MatrixBaseSeq<T, MatrixRes<T> >
{
 public:
     template<class ContainerT>
     MatrixResSeq(const MatrixBaseSeq<T, ContainerT >& ms):\
       MatrixBaseSeq<T, MatrixRes<T> >()
     {
         #ifdef DIAG
         cout << "MatrixResSeq constructor" << endl;
         #endif
		assign(ms.size(), ms[0]);
         for(register int i=0;i<ms.size();i++)
         {
             MatrixRes<T> m(ms[i]);
			 (*this)[i] = m;
             (*this)[i].set_copyoption(SHALLOW);
         }
     }
     
}; 

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixSubSeq

// Helper class for obtaining matrix sub sequences from a matrix seq by using
// the operator[](all) on the MatrixSeq object and returning an object of this
// sub-type
// Utility class only. Transparent to users
//------------------------------------------------------------------------------

template<class T>
class MatrixSubSeq:public MatrixSeq<T>
{
private:
  MatrixSubSeq(const MatrixSeq<T>& ms):MatrixSeq<T>(ms.size())
  {
     for(register int i=0;i<ms.size();i++)
       this->elt(i)=MatrixRes<T>(ms[i]);
  } 
   
  friend class MatrixSeq<T>;

public:
  MatrixSubSeq<T> operator[](size_t idx) const;
  MatrixSubSeq<T> operator[](const DummyArg& arg) const;
  MatrixSubSeq<T> operator[](const slice& sl) const;
  
  // --- assignments
  MatrixSubSeq& operator=(const MatrixSeq<T>& seq)
  { assign(seq); return *this;}
  template<class U>
  MatrixSubSeq& operator=(const vector<U>& v)
  { assign(v); return *this;}
  template<class U>
  MatrixSubSeq& operator=(const Matrix<U>& s)
  { assign(this->size(),s); return *this;}
  MatrixSubSeq& operator=(const T& t)
  { assign(this->size(),t); return *this;}
};
// ----------- End of class MatrixResSeq -------------------

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class ScalarSeq

// Class for a sequence of Scalars
//------------------------------------------------------------------------------
template<class T>
class ScalarSeq : public MatrixBaseSeq<T, Scalar<T> > 
{
public:
    // --- constructors
    ScalarSeq(int num=0):MatrixBaseSeq<T, Scalar<T> >(num){;}
    ScalarSeq(int num,const T& t):\
      MatrixBaseSeq<T, Scalar<T> >(num,t){;}
    template<class U>
    ScalarSeq(int num,const Scalar<U>& s):\
      MatrixBaseSeq<T, Scalar<T> >( num,s){;}
    ScalarSeq(const ScalarSeq<T>& seq):MatrixBaseSeq<T, Scalar<T> >(seq){;}
    template<class U>
    ScalarSeq(const vector<U>& v):MatrixBaseSeq<T, Scalar<T> >(v){;}
      
    // --- assignments
    ScalarSeq& operator=(const ScalarSeq<T>& seq)
    { assign(seq); return *this;}
    template<class U>
    ScalarSeq& operator=(const vector<U>& v)
    { assign(v); return *this;}
    template<class U>
    ScalarSeq& operator=(const Scalar<U>& s)
    { assign(this->size(),s); return *this;}
    ScalarSeq& operator=(const T& t)
    { assign(this->size(),t); return *this;}
        
    virtual bool is_scalar_seq() const
    { return true;}
};
// ----------- End of class ScalarSeq -------------------

}  // end namespace QDOS
    
#endif
