//
//  Matrix.cc
//
//  Qdos - member functions and operators for Numeric Matrix Class
//
//  Janet Barnes
//  Sept 1999


#ifndef QDOS_MATRIX_CC
#define QDOS_MATRIX_CC

#include "Matrix.h"
#include "MatrixSVD.h"
#include "MatrixEig.h"

namespace QDOS{


// ++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::matrix_to_seq
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------- matrix_to_seq
template<class T>
MatrixResSeq<T> Matrix<T>::matrix_to_seq(size_t dimno) const
{
    if(dimno<=0 || dimno>rank())
        throw MatrixIndexError();
    size_t sz=dims()[dimno-1];
    MatrixSeq<T> ms(sz);
    for(register int i=0;i<sz;i++)
        ms[i]=indexdim(dimno,i);
    return MatrixResSeq<T>(ms);
}       
// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::matrix_indexelt
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ matrix_indexelt
template<class T>
inline T& Matrix<T>::matrix_indexelt(size_t *ind,size_t nel) const
{
    size_t curr=inner_product(&ind[0],&ind[nel],&strides()[0],start());
    return (*_val)[curr];
}

template<class T>
inline T& Matrix<T>::matrix_indexelt(valarray<size_t> v) const
{
    int sz=rank();
    if (v.size()!=sz)
        throw MatrixLengthError(\
            "[Matrix<T>::matrix_indexelt(valarray<size_t>)]");
    valarray<size_t> str=strides();
    size_t st=start();
    size_t curr=inner_product(&v[0],&v[sz],&str[0],st);
    return (*_val)[curr];
}

template<class T>
inline T& Matrix<T>::matrix_indexelt(const vector<int>& v) const
{
    int sz=rank();
    if (v.size()!=sz)
        throw MatrixLengthError(\
            "[Matrix<T>::matrix_indexelt(vector<int>)]");
    valarray<size_t> str=strides();
    int st=start();
    int curr=inner_product(v.begin(),v.end(),&str[0],st);
    return (*_val)[curr];
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::diagonal
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- diagonal
template<class T>
inline Matrix<T> Matrix<T>::diagonal() const
{
// Return diagonal elements of Matrix i.e (0,0,..0),(1,1,...1) ...
// (n,n,...,n) where n is size of smallest dimension
    DimStorageOption order=storageorder();
    valarray<size_t> sz(dims().min(),1);
    valarray<size_t> st(strides().sum(),1);
    Dim_var dim=new Dim(start(),sz,st,order);
    return Matrix<T>(data_objvar(),dim,cdimno(),true,DEEP,!empty());
}
// ++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::null_matrix
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ---------------------------------------------------------------- null_matrix
template<class T>
inline void Matrix<T>::null_matrix() const
{
// Set current matrix to all zero
    T zero=0;
    assign(zero);
}
// ++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::unit_matrix
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ---------------------------------------------------------------- unit_matrix
template<class T>
inline void Matrix<T>::unit_matrix() const
{
// Make current matrix object identity matrix i.e. set all values
// to zero and diagonal entries to 1
    T zero=0;
    T one=1;
    assign_val(zero);  // set all entries to zero
    diagonal()=one;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::transpose
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- transpose
template<class T>
inline void Matrix<T>::transpose() 
{
// Transpose current Matrix i.e. reverse axes of array
     valarray<size_t> dm(dims());
     valarray<size_t> st(strides());
     reverse(&dm[0],&dm[dm.size()]);
     reverse(&st[0],&st[st.size()]);
     DimStorageOption order=\
         (storageorder()==ROW_MAJOR_ORDER)?COL_MAJOR_ORDER:ROW_MAJOR_ORDER;
     if(rank()!=0)
       _diminfo=new Dim(start(),dm,st,order);
}
              
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::trace
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ---------------------------------------------------------------------- trace
template<class T>
inline T Matrix<T>::trace() const
{
// return the trace of a matrix
    return diagonal().sum();
}
// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::duplicate
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- duplicate
template<class T>
inline MatrixRes<T> Matrix<T>::duplicate(MatrixCopyOption option) const
{
    if (option==SHALLOW)
        return MatrixRes<T>(*this);
    else
    {
        // DEEP copy required
        Matrix<T> m(dimdata());
        m.copy(*this);
        return MatrixRes<T>(m);
    }
 }
// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::is_symmetric
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- is_symmetric
template<class T> 
inline bool Matrix<T>::is_symmetric() const
{
// is matrix identical to its transpose?
    Matrix<T> m;
    m.clone(*this); 
    m.transpose();
    return (m==*this);
}       
// ++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::norm_svd
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ norm_svd
template<class T> 
inline T norm_svd(const Matrix<T>& m)  
{
// calculate the largest singular value of a 2D matrix    
    if(m.rank()==2)
    {    
         MatrixSVD<T> svd(m);
         const Matrix<double>& w=svd.get_w();
         double max=*(max_element(w.begin(),w.end()));
         return (T)max;
    }
    else
        throw MatrixRankError(); 
} 
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ abs 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- abs  
template<class T> 
Scalar<T> abs(const Scalar<T>& s)
{ return Scalar<T>(std::abs(s.scalar()));}
 
template<class T> 
MatrixRes<T> abs(const Matrix<T>& m)
{  
    Matrix<T> mres(m);   
    (*mres.data_ptr())=abs(*m.data_ptr());
    return MatrixRes<T>(mres);
}

template<class T> 
MatrixRes<complex<T> > abs(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m);   
    for(register int i=0;i<m.size();i++)
        (*mres.data_ptr())[i]=std::abs((*m.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<T> abs(const MatrixBaseSeq<T, ContainerT >& ms)
{
    MatrixSeq<T> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=abs(ms[i]);
    return MatrixResSeq<T>(msres);   
}
 
// ++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::norm
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ norm
template<class T> 
inline T norm(const Scalar<T>& s)  
{
  return abs(s).scalar();
}
          
template<class T> 
inline T norm(const Matrix<T>& m)  
{
// calculate the norm of a vector or 2D matrix 

    T res;    
    switch(m.rank())
    {
        case 0: // trivial if scalars
            res=abs(m).scalar(); 
            break;              
          
        case 1:
            res=sqrt((Matrix<T>)m.dot(m));
            break;
            
        case 2: 
           res=norm_svd(m);
           break;
                        
        default:
        // currently not defined for matrices > 2D...
            throw MatrixDomainError();   
    }
        
    return res;
}
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::conj
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------------ conj
template<class T>
Scalar<typename promote_trait<T,float_complex>::T_promote> conj(const Scalar<T>& s)
{
  typedef typename promote_trait<T,float_complex>::T_promote T_promote;
  return Scalar<T_promote>(std::conj(s.scalar()));
}

template<class T>
MatrixRes<typename promote_trait<T,float_complex>::T_promote> conj(const Matrix<T>& m)
{
    typedef typename promote_trait<T,float_complex>::T_promote T_promote;
    Matrix<T_promote> mres=m.duplicate();
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=conj((*mres.data_ptr())[i]);
    return MatrixRes<T_promote>(mres);
}
// +++++++++++++++++++++++++++++++++++++++++++++ class Matrix::set_storageorder
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------- set_storageorder
template<class T>
inline void Matrix<T>::set_storageorder(DimStorageOption order)
{if(rank()!=0) _diminfo=new Dim(getgslice(),order);}

// +++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::storage
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- storage
template<class T> 
inline void Matrix<T>::storage(DimStorageOption order)
{
// Force matrix to be stored in specified order (row-major or col-major)
    if (storageorder()!=order)
    {
        transpose();
        make_contiguous();
        transpose();
        set_storageorder(order);
     }
}
// ++++++++++++++++++++++++++++++++++++++++++++ class Matrix::fortranarray
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ---------------------------------------------------------- fortranarray
template<class T> 
MatrixRes<T> Matrix<T>::fortranarray() const
{
// return copy of matrix stored in COL_MAJOR_ORDER
    Matrix<T> m(Dim(this->dims(),COL_MAJOR_ORDER));
    m.copy(*this);
    return MatrixRes<T>(m);
}              
// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator[]
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator[]

template<class T>
Matrix<T> Matrix<T>::operator[](const DummyArg& arg) const
{
    return Matrix<T>(data_objvar(),new Dim(getgslice(),storageorder()),\
        cdimno()+1,true,DEEP,!empty());
}

template<class T>
Matrix<T> Matrix<T>::operator[](const size_t i) const
{
     return indexdim(cdimno(),i);   
}    

template<class T>
Matrix<T> Matrix<T>::operator[](const slice& sl) const
{
     Matrix<T> m;
     m.clone(*this);
     m.dimno()+=1;
     return m.indexdim(cdimno(),sl);
}

template<class T>
Matrix<T> Matrix<T>::operator[](const gslice& g) const
{
// --- returns a sub-matrix of 
    Dim_var dim_var=new Dim(g,storageorder());                          
    return Matrix<T>(data_objvar(),dim_var,1,true,DEEP,!empty());
} 
// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator()
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator()

// --- Rank 1 Matrices
template<class T>
T& Matrix<T>::operator()(const size_t i) const
{
// --- returns elt(i) of a 1 D Matrix
    if (rank()!=1)
         throw MatrixRankError();
     size_t ind[1]={i};
    return matrix_indexelt(ind,1);
}

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl) const
{
// --- returns a sub-array of a 1D Matrix
    if (rank()!=1)
        throw MatrixRankError();

}
 
// --- Rank 2 Matrices
template<class T>
T& Matrix<T>::operator()(const size_t i,const size_t j) const
{
// --- returns elt (i,j) of Matrix

    if(rank()!=2)
        throw MatrixRankError();
    size_t ind[2]={i,j};
    return matrix_indexelt(ind,2);

}   
   
template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const slice& sl) const
{
// --- returns a 1D sub-matrix of Matrix[i;sl.start..sl.start+sl.size*sl.stride]
    if(rank()!=2)
        throw MatrixRankError();
    return row(i).row(sl);    

}

template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const DummyArg& arg) const
{
// --- returns a 1D sub-matrix of Matrix[i;]
    if(rank()!=2)
        throw MatrixRankError();
    return row(i);    

}

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl,const DummyArg& arg) const
{
// --- returns a 1D sub-matrix of Matrix[i;]
    if(rank()!=2)
        throw MatrixRankError();
    return row(sl);    

}

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl,const size_t j) const
{
 // --- returns a 1D sub-matrix of Matrix[sl.start..sl.start+sl.size*sl.stride;
 //     j]
    if(rank()!=2)
        throw MatrixRankError();
    return row(sl).col(j);    

}

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg,const size_t j) const
{
 // --- returns a 1D sub-matrix of Matrix[;j]
    if(rank()!=2)
        throw MatrixRankError();
    return col(j);    
}

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg,const slice& sl) const
{
 // --- returns a 1D sub-matrix of Matrix[;j]
   if(rank()!=2)
       throw MatrixRankError();
   return col(sl);    

}

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl1,const slice& sl2) const
{
// --- returns a 2D sub-matrix of Matrix[sl1.start..sl1.start+sl1.size*
//       sl1.stride;sl2.start..sl2.start+sl2.size*sl2.stride]
   if(rank()!=2)
       throw MatrixRankError();
    return row(sl1).col(sl2);

}

// --- Rank 3 Matrices
template<class T>
T& Matrix<T>::operator()(const size_t i,const size_t j,const size_t k) const
{
// --- returns elt (i,j,k) of Matrix
   if(rank()!=3)
       throw MatrixRankError();
   size_t ind[3]={i,j,k};
   return matrix_indexelt(ind,3); 

}

template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const size_t j,\
                                    const slice& sl) const
{
// --- returns a 1D sub-matrix of Matrix[i;j;sl.start..sl.start+
//     sl.size*sl.stride]
   if(rank()!=3)
       throw MatrixRankError();
   return row(i).row(j).row(sl);

}        

template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const size_t j,\
                                const DummyArg& arg) const
{
// --- returns a 1D sub-matrix of Matrix[i;j;]
   if(rank()!=3)
       throw MatrixRankError();
   return row(i).row(j);

}     

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl,const size_t j,\
                                const DummyArg& arg) const
{
// --- returns a 2D sub-matrix of Matrix[sl;j;]
   if(rank()!=3)
       throw MatrixRankError();
   return row(sl).col(j);

}     

template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i, const slice& sl,\
                                const DummyArg& arg) const
{
// --- returns a 2D sub-matrix of Matrix[i;sl;]
   if(rank()!=3)
       throw MatrixRankError();
   return row(i).row(sl);

}     

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl1, const slice& sl2,\
                                const DummyArg& arg) const
{
// --- returns a 3D sub-matrix of Matrix[sl1;sl2;]
   if(rank()!=3)
       throw MatrixRankError();
   return row(sl1).col(sl2);

}     
      
template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const slice& sl,const size_t k) const
{
// --- returns a 1D sub-matrix of Matrix[i;sl.start..sl.start+sl.size*sl.stride;
//     k]
   if(rank()!=3)
       throw MatrixRankError();
   return row(i).row(sl).col(k);

}
      
template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const DummyArg& arg,\
                                const size_t k) const
{
// --- returns a 1D sub-matrix of Matrix[i;;k]
   if(rank()!=3)
       throw MatrixRankError();
   return row(i).col(k);

}      

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl,const DummyArg& arg,\
                                const size_t k) const
{
// --- returns a 2D sub-matrix of Matrix[sl;;k]
   if(rank()!=3)
       throw MatrixRankError();
   return row(sl).layer(k);

}

template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const DummyArg& arg,\
                                const slice& sl) const
{
// --- returns a 2D sub-matrix of Matrix[i;;sl]
   if(rank()!=3)
       throw MatrixRankError();
   return row(i).col(sl);

}            

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl1,const DummyArg& arg,\
                                const slice& sl2) const
{
// --- returns a 3D sub-matrix of Matrix[sl1;;sl2]
   if(rank()!=3)
       throw MatrixRankError();
   return row(sl1).layer(sl2);

}            

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl,const size_t j,const size_t k) const
{
// --- returns a 1D sub-matrix of Matrix[sl.start..sl.start+sl.size*sl.stride;
//     j;k]
   if(rank()!=3)
       throw MatrixRankError();
   return MatrixRes<T>(row(sl).col(j).col(k));

}  

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg,const size_t j,\
                                const size_t k) const
{
// --- returns a 1D sub-matrix of Matrix[;j;k]
   if(rank()!=3)
       throw MatrixRankError();
   return col(j).col(k);

} 

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg,const slice& sl,\
                                const size_t k) const
{
// --- returns a 2D sub-matrix of Matrix[;sl;k]
   if(rank()!=3)
       throw MatrixRankError();
   return col(sl).layer(k);

} 

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg,const size_t j,\
                                const slice& sl) const
{
// --- returns a 2D sub-matrix of Matrix[;j;sl]
   if(rank()!=3)
       throw MatrixRankError();
   return col(j).col(sl);

}

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg,const slice& sl1,\
                                const slice& sl2) const
{
// --- returns a 3D sub-matrix of Matrix[;sl1;sl2]
   if(rank()!=3)
       throw MatrixRankError();
   return col(sl1).layer(sl2);

}   

template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const slice& sl1,const slice& sl2) const
{
// --- returns a 2D sub-matrix of 
//     Matrix[i;sl1.start..sl1.start+sl1.size*sl1.stride;
//     sl2.start..sl2.start+sl2.size*sl2.stride]
   if(rank()!=3)
       throw MatrixRankError();
   return row(i).row(sl1).col(sl2);

}  

template<class T>
Matrix<T> Matrix<T>::operator()(const size_t i,const DummyArg& arg1,\
                                const DummyArg& arg2) const
{
// --- returns a 2D sub-matrix of Matrix[i;;]
   if(rank()!=3)
       throw MatrixRankError();
   return row(i);

}  

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl,const DummyArg& arg1,\
                                const DummyArg& arg2) const
{
// --- returns a 3D sub-matrix of Matrix[sl;;]
   if(rank()!=3)
       throw MatrixRankError();
   return row(sl);

}  

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl1,const size_t j,const slice& sl2) const
{
// --- returns a 2D sub-matrix of
//     Matrix[sl1.start..sl1.start+sl1.size*sl1.stride;j;
//     sl2.start..sl2.start+sl2.size*sl2.stride]
   if(rank()!=3)
       throw MatrixRankError();
   return row(sl1).col(j).col(sl2);

}    

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg1,const size_t j,\
                                const DummyArg& arg2) const
{
// --- returns a 2D sub-matrix of Matrix[;j;]
   if(rank()!=3)
       throw MatrixRankError();
   return col(j);

}

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg1,const slice& sl,\
                                const DummyArg& arg2) const
{
// --- returns a 3D sub-matrix of Matrix[;sl;]
   if(rank()!=3)
       throw MatrixRankError();
   return col(sl);

}

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl1,const slice& sl2,const size_t k) const
{
// --- returns a 2D sub-matrix of
//     Matrix[sl1.start..sl1.start+sl1.size*sl1.stride;
//     sl2.start..sl2.start+sl2.size*sl2.stride;k]
   if(rank()!=3)
       throw MatrixRankError();
   return row(sl1).col(sl2).layer(k);

} 

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg1,const DummyArg& arg2,\
                                const size_t k) const
{
// --- returns a 2D sub-matrix of Matrix[;;k]
   if(rank()!=3)
       throw MatrixRankError();
   return layer(k);

}

template<class T>
Matrix<T> Matrix<T>::operator()(const DummyArg& arg1,const DummyArg& arg2,\
                                const slice& sl) const
{
// --- returns a 3D sub-matrix of Matrix[;;sl]
   if(rank()!=3)
       throw MatrixRankError();
    return layer(sl);

}

template<class T>
Matrix<T> Matrix<T>::operator()(const slice& sl1,const slice& sl2,\
    const slice& sl3) const
{
// --- returns a 3D sub-matrix of 
//     Matrix[sl1.start..sl1.start+sl1.size*sl1.stride;
//     sl2.start..sl2.start+sl2.size*sl2.stride;sl3.start..
//     sl3.start+sl3.size*sl3.stride]
   if(rank()!=3)
       throw MatrixRankError();
   return row(sl1).col(sl2).layer(sl3);

} 

template<class T>
Matrix<T> Matrix<T>::operator()(const gslice& g) const
{
// --- returns a sub-matrix of 
    Dim_var dim_var=new Dim(g,storageorder());                          
    return Matrix<T>(data_objvar(),dim_var,1,true,DEEP,!empty());
} 
// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator!
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator!
// operator for conjugate transpose
template<class T>
MatrixRes<T> operator!(const Matrix<T>& m)
{   return ~m;}

template<class T>
MatrixRes<complex<T> > operator!(const Matrix<complex<T> >& m)
{ return conj(~m);}
// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator~
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator~
// operator for transpose
template<class T>
MatrixRes<T> Matrix<T>::operator~() const
{
    Matrix<T> mres=duplicate();
    mres.transpose();
    return MatrixRes<T>(mres);
}
// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator^
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator^

template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> \
Scalar<T>::operator^(const double& n) const
{
     typedef typename promote_trait<T,double>::T_promote T_promote;   
     T_promote p(n);      
     return Scalar<T_promote>(std::pow((T_promote)scalar(),p));
}

// operator for power - what about complex nos.
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> \
Matrix<T>::operator^(const double& n) const
{ return power(n);}

template<class T> template<class U>
ScalarSeq<typename promote_trait<T,double>::T_promote> \
    Scalar<T>::operator^(const vector<U>& v) const
{
    typedef typename promote_trait<T,double>::T_promote T_promote;   
    ScalarSeq<T_promote> ssres(v.size());
    for(register int i=0;i<v.size();i++)
        ssres[i]=this->operator^((double)v[i]);
    return ssres;
}   

template<class T> template<class U>
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    Matrix<T>::operator^(const vector<U>& v) const
{
    typedef typename promote_trait<T,double>::T_promote T_promote;   
    MatrixSeq<T_promote> msres(v.size());
    for(register int i=0;i<v.size();i++)
        msres[i]=power((double)v[i]);
    return MatrixResSeq<T_promote>(msres);   
}
// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator+
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator+
// operator for identity
template<class T>    
// --- operator for identity 
MatrixRes<T> Matrix<T>::operator+() const
{
    Matrix<T> m=duplicate();
    return MatrixRes<T>(m);
}

template<class U, class T>
MatrixRes<typename promote_trait<T,U>::T_promote> add_scalar(const Matrix<T>& m,const U& u)
{
// Scalar addition
    typedef typename promote_trait<T,U>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());
    *mres.data_ptr()+=(T_promote)u;
    return MatrixRes<T_promote>(mres);
}

// operator for addition
template<class T>
MatrixRes<T> Matrix<T>::operator+(const T& t) const
{ return add_scalar(*this,t);}

template<class T> template<class U>
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::operator+(const Scalar<U>& s) const
{ return add_scalar(*this,s.scalar());}

template<class T>
MatrixRes<T> Matrix<T>::operator+(const Matrix<T>& m) const
{
// --- Matrix addition
// Current Matrix being pointed to by this and Matrix m1 must both be of same 
// rank and dims must be the same
   if(rank()==0)
        return add_scalar(m,scalar());
    if(m.rank()==0)
        return add_scalar(*this,m.scalar());
    Matrix<T> mres=duplicate();
    mres+=m;
    return MatrixRes<T>(mres);
}

template<class T> template<class U>
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::operator+(const Matrix<U>& m) const
{
    Matrix<typename promote_trait<T,U>::T_promote> m1(*this);
    Matrix<typename promote_trait<T,U>::T_promote> m2(m);
    return m1+m2;
}

template<class T>
Scalar<typename promote_trait<T,int>::T_promote> operator+(const int& i,\
     const Scalar<T>& s)
{ typedef typename promote_trait<T,int>::T_promote T_promote;
  T_promote s1(s);
  T_promote s2(i);
  return Scalar<T_promote>(s1+s2);
}

template<class T>
MatrixRes<typename promote_trait<T,int>::T_promote> operator+(const int& i,\
     const Matrix<T>& m)
{ return add_scalar(m,i);}

template<class T>
Scalar<typename promote_trait<T,double>::T_promote> operator+(const double& d,\
     const Scalar<T>& s)
{ typedef typename promote_trait<T,double>::T_promote T_promote;
  T_promote s1(s);
  T_promote s2(d);
  return Scalar<T_promote>(s1+s2);
}
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> operator+(const double& d,\
    const Matrix<T>& m)
{ return add_scalar(m,d);}

template<class U, class T>
Scalar<typename promote_trait<T,complex<U> >::T_promote> operator+(const complex<U>& c,\
     const Scalar<T>& s)
{ typedef typename promote_trait<T,complex<U> >::T_promote T_promote;
  T_promote s1(s);
  T_promote s2(c);
  return Scalar<T_promote>(s1+s2);
}

template<class U, class T>
MatrixRes<typename promote_trait<T,complex<U> >::T_promote> operator+(const complex<U>& c,\
    const Matrix<T>& m)
{ return add_scalar(m,c);}

template<class T>
Scalar<complex<T> > operator+(const Scalar<complex<T> >& s, const double& d)
{
  return Scalar<complex<T> >(s.scalar()+d);
}

template<class T>
MatrixRes<complex<T> > operator+(const Matrix<complex<T> >& m,const double& d)
{ return add_scalar(m,d);}

template<class T>
Scalar<complex<T> > operator+(const Scalar<complex<T> >& s, const int& i)
{
  return Scalar<complex<T> >(s.scalar()+i);
}

template<class T>
MatrixRes<complex<T> > operator+(const Matrix<complex<T> >& m,const int& i)
{ return add_scalar(m,i);}
// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator-
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator-

// operator for negation
template<class T>   
MatrixRes<T> Matrix<T>::operator-() const
{
    Matrix<T> m(dimdata(),-*_val);
    return MatrixRes<T>(m);
}

template<class U, class T>
MatrixRes<typename promote_trait<T,U>::T_promote> sub_scalar(const Matrix<T>& m,const U& u)
{
// Scalar subtraction
    typedef typename promote_trait<T,U>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());
    *mres.data_ptr()-=(T_promote)u;
    return MatrixRes<T_promote>(mres);
}

// Scalar subtraction
template<class T>
MatrixRes<T> Matrix<T>::operator-(const T& t) const
{ return sub_scalar(*this,t);}

template<class T> template<class U>
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::operator-(const Scalar<U>& s) const
{ return sub_scalar(*this,s.scalar());}

template<class T> 
MatrixRes<T> Matrix<T>::operator-(const Matrix<T>& m) const
{
// --- Matrix subtraction
// Current Matrix being pointed to by this and Matrix m1 must both be of same 
// rank and dims must be the same
   if(rank()==0)
       return add_scalar(-m,scalar());
   if(m.rank()==0)
       return sub_scalar((*this),m.scalar());
   Matrix<T> mres=duplicate();
   mres-=m;
   return MatrixRes<T>(mres);  
}

template<class T> template<class U>
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::operator-(const Matrix<U>& m) const
{
    Matrix<typename promote_trait<T,U>::T_promote> m1(*this);
    Matrix<typename promote_trait<T,U>::T_promote> m2(m);
    return m1-m2;
}

template<class T>
Scalar<typename promote_trait<T,int>::T_promote> operator-(const int& i,\
     const Scalar<T>& s)
{ typedef typename promote_trait<T,int>::T_promote T_promote;
  T_promote s2(i);
  T_promote s1(s);
  return Scalar<T_promote>(s2.scalar()-s1.scalar());
}

template<class T>
MatrixRes<typename promote_trait<T,int>::T_promote> operator-(const int& i,\
     const Matrix<T>& m)
{ return add_scalar(-m,i);}

template<class T>
Scalar<typename promote_trait<T,double>::T_promote> operator-(const double& d,\
     const Scalar<T>& s)
{ typedef typename promote_trait<T,double>::T_promote T_promote;
  T_promote s2(d);
  T_promote s1(s);
  return Scalar<T_promote>(s2.scalar()-s1.scalar());
}

template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> operator-(const double& d,\
    const Matrix<T>& m)
{ return add_scalar(-m,d);}

template<class U, class T>
Scalar<typename promote_trait<T,complex<U> >::T_promote> operator-(const complex<U>& c,\
     const Scalar<T>& s)
{ typedef typename promote_trait<T,complex<U> >::T_promote T_promote;
  T_promote s2(c);
  T_promote s1(s);
  return Scalar<T_promote>(s2.scalar()-s1.scalar());
}

template<class U, class T>
MatrixRes<typename promote_trait<T,complex<U> >::T_promote> operator-(const complex<U>& c,\
    const Matrix<T>& m)
{ return add_scalar(-m,c);}

template<class T>
Scalar<complex<T> > operator-(const Scalar<complex<T> >& s, const double& d)
{
  return Scalar<complex<T> >(s.scalar()-d);
}

template<class T>
MatrixRes<complex<T> > operator-(const Matrix<complex<T> >& m,const double& d)
{ return sub_scalar(m,d);}

template<class T>
Scalar<complex<T> > operator-(const Scalar<complex<T> >& s, const int& i)
{
  return Scalar<complex<T> >(s.scalar()-i);
}

template<class T>
MatrixRes<complex<T> > operator-(const Matrix<complex<T> >& m,const int& i)
{ return sub_scalar(m,i);}
// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator+=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator+=

template<class T>   
Matrix<T>& Matrix<T>::operator+=(const T& t)
{
// --- Matrix+=Scalar
   // uses gslice_array computed assignment...rhs must be a valarray
    valarray<T> v(t,size());
    (*_val)[getgslice()]+=v;
    return *this;
}

template<class T> template<class U>  
Matrix<T>& Matrix<T>::operator+=(const Scalar<U>& s)
{ return this->operator+=((T)s.scalar());}
    
template<class T>   
const Matrix<T>& Matrix<T>::operator+=(const T& t) const
{
// --- Matrix+=Scalar
   // uses gslice_array computed assignment...rhs must be a valarray
    valarray<T> v(t,size());
    (*_val)[getgslice()]+=v;
    return *this;
}

template<class T> template<class U>  
const Matrix<T>& Matrix<T>::operator+=(const Scalar<U>& s) const
{ return this->operator+=((T)s.scalar());}

template<class T> 
Matrix<T>& Matrix<T>::operator+=(const Matrix<T>& m) 
{
// --- Matrix+=Matrix addition
// Current Matrix being pointed to by this and Matrix m must both be of same 
// rank and dims
   if(m.rank()==0)
       return (*this)+=m.scalar();
   if (!samesize(m))
       throw MatrixLengthError();
   // uses gslice_array +=
   (*_val)[getgslice()]+=(*m._val)[m.getgslice()]; 
   return *this; 
}

template<class T>
const Matrix<T>& Matrix<T>::operator+=(const Matrix<T>& m) const
{
// --- Matrix+=Matrix addition
// Current Matrix being pointed to by this and Matrix m must both be of same 
// rank and dims
   if(m.rank()==0)
       return (*this)+=m.scalar();
   if (!samesize(m))
       throw MatrixLengthError();
   // uses gslice_array +=
   (*_val)[getgslice()]+=(*m._val)[m.getgslice()]; 
   return *this;
}

template<class T> template<class U>
Matrix<T>& Matrix<T>::operator+=(const Matrix<U>& m)
{ Matrix<T> mu(m); return this->operator+=(mu);}

template<class T> template<class U>
const Matrix<T>& Matrix<T>::operator+=(const Matrix<U>& m) const
{ Matrix<T> mu(m); return this->operator+=(mu);}

// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator-=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator-=

template<class T>
Matrix<T>& Matrix<T>::operator-=(const T& t) 
{
// --- Matrix-=Scalar
    // uses gslice_array computed assignment...rhs must be a valarray
    valarray<T> v(t,size());
    (*_val)[getgslice()]-=v;
    return *this;
}

template<class T> template<class U>
Matrix<T>& Matrix<T>::operator-=(const Scalar<U>& s) 
{ return this->operator-=((T)s.scalar());}

template<class T>
const Matrix<T>& Matrix<T>::operator-=(const T& t) const
{
// --- Matrix-=Scalar
    // uses gslice_array computed assignment...rhs must be a valarray
    valarray<T> v(t,size());
    (*_val)[getgslice()]-=v;
    return *this;
}

template<class T> template<class U>
const Matrix<T>& Matrix<T>::operator-=(const Scalar<U>& s) const
{ return this->operator-=((T)s.scalar());}

template<class T>
Matrix<T>& Matrix<T>::operator-=(const Matrix<T>& m)
{
// --- Matrix-=Matrix 
// Current Matrix being pointed to by this and Matrix m1 must both be of same 
// rank and dims
   if(m.rank()==0)
   {
       (*this)-=m.scalar();
       return *this;
   }
   if (!samesize(m))
       throw MatrixLengthError();
   // uses gslice_array -=
   (*_val)[getgslice()]-=(*m._val)[m.getgslice()]; 
   return *this;      
}

template<class T>
const Matrix<T>& Matrix<T>::operator-=(const Matrix<T>& m) const
{
// --- Matrix-=Matrix 
// Current Matrix being pointed to by this and Matrix m1 must both be of same 
// rank and dims
   if(m.rank()==0)
       return (*this)-=m.scalar();
   if (!samesize(m))
       throw MatrixLengthError(); 
   // uses gslice_array -=
   (*_val)[getgslice()]-=(*m._val)[m.getgslice()]; 
    return *this;
}

template<class T> template<class U>
Matrix<T>& Matrix<T>::operator-=(const Matrix<U>& m)
{ Matrix<T> mu(m); return this->operator-=(mu);}

template<class T> template<class U>
const Matrix<T>& Matrix<T>::operator-=(const Matrix<U>& m) const
{ Matrix<T> mu(m); return this->operator-=(mu);}

// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator*
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator*
template<class U, class T>
MatrixRes<typename promote_trait<T,U>::T_promote> mult_scalar(const Matrix<T>& m,const U& u)
{
// Scalar multiplication
    typedef typename promote_trait<T,U>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());
    *mres.data_ptr()*=(T_promote)u;
    return MatrixRes<T_promote>(mres);
}

template<class T>   
MatrixRes<T> Matrix<T>::operator*(const T& t) const
{ return mult_scalar(*this,t);}

template<class T> template<class U>  
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::operator*(const Scalar<U>& s) const
{ return mult_scalar(*this,s.scalar());}

template<class T>
MatrixRes<T> Matrix<T>::operator*(const Matrix<T>& m) const
{
// --- Matrix multiplication - will also perform Matrix*
// (const vector<T,allocator<T> > since latter will be converted 
// to a matrix by Matrix(Vector<T,allocator<T> >& v)
 // invoke scalar multiplication if m is rank 0
    if (m.rank()==0)
        return mult_scalar(*this,m.scalar()); 
    else if (rank()==0)
    {   
        // invoke scalar multiplication if this is rank 0             
        return mult_scalar(m,scalar()); 
    }
    else                
    {
        // change to vector?
        valarray<size_t> d(dims());
        valarray<size_t> md(m.dims());
        // check that matrices are compatible
        if (d[rank()-1]!=md[0])
            throw MatrixLengthError();
        const size_t nrank=rank()+m.rank()-2;
        size_t *dm=new size_t[nrank];
        register size_t i,j;
            
        for(i=0;i<rank()-1;i++) dm[i]=d[i];
        for(i=1;i<m.rank();i++) dm[rank()-2+i]=md[i];
          
        Matrix<T> mres(Dim(dm,nrank),(T)0);
        
        delete [] dm;     
        matmult(m,mres);   // ... call generic matrix mult routine  
        return MatrixRes<T>(mres);
    }
}

template<class T> template<class U>
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::operator*(const Matrix<U>& m) const
{
    Matrix<typename promote_trait<T,U>::T_promote> m1(*this);
    Matrix<typename promote_trait<T,U>::T_promote> m2(m);
    return m1*m2;
}

template<class U,class T>
MatrixRes<typename promote_trait<T,U>::T_promote> operator*(const vector<U>& v,\
    const Matrix<T>& m)
{
    Matrix<typename promote_trait<T,U>::T_promote> m1(v); 
    Matrix<typename promote_trait<T,U>::T_promote> m2(m); 
    return m1*m2;
}

template<class T>
Scalar<typename promote_trait<T,int>::T_promote> operator*(const int& i,\
     const Scalar<T>& s)
{
  typedef typename promote_trait<T,int>::T_promote T_promote;
  Scalar<T_promote> s1(i);
  Scalar<T_promote> s2(s);
  return Scalar<T_promote>(s1.scalar()*s2.scalar());
}

template<class T>
MatrixRes<typename promote_trait<T,int>::T_promote> operator*(const int& i,\
     const Matrix<T>& m)
{ return mult_scalar(m,i);}

template<class T>
Scalar<typename promote_trait<T,double>::T_promote> operator*(const double& d,\
     const Scalar<T>& s)
{
  typedef typename promote_trait<T,double>::T_promote T_promote;
  Scalar<T_promote> s1(d);
  Scalar<T_promote> s2(s);
  return Scalar<T_promote>(s1.scalar()*s2.scalar());
}

template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> operator*(const double& d,\
    const Matrix<T>& m)
{ return mult_scalar(m,d);}

template<class U, class T>
Scalar<typename promote_trait<T,complex<U> >::T_promote> operator*(const complex<U>& c,\
     const Scalar<T>& s)
{
  typedef typename promote_trait<T,complex<U> >::T_promote T_promote;
  Scalar<T_promote> s1(c);
  Scalar<T_promote> s2(s);
  return Scalar<T_promote>(s1.scalar()*s2.scalar());
}

template<class U, class T>
MatrixRes<typename promote_trait<T,complex<U> >::T_promote> operator*(const complex<U>& c,\
    const Matrix<T>& m)
{ return mult_scalar(m,c);}

template<class T>
Scalar<complex<T> > operator*(const Scalar<complex<T> >& s, const int& i)
{
  return Scalar<complex<T> >(s.scalar()*i);
}

template<class T>
MatrixRes<complex<T> > operator*(const Matrix<complex<T> >& m,const double& d)
{ return mult_scalar(m,d);}

template<class T>
Scalar<complex<T> > operator*(const Scalar<complex<T> >& s, const double& d)
{
  return Scalar<complex<T> >(s.scalar()*d);
}

template<class T>
MatrixRes<complex<T> > operator*(const Matrix<complex<T> >& m,const int& i)
{ return mult_scalar(m,i);}
// ++++++++++++++++++++++++++++++++++++++++++++ class Matrix::mult_element
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ---------------------------------------------------------- mult_element
template<class T>
MatrixRes<T> Matrix<T>::mult_element(const Matrix<T>& m) const
{
// --- Multiply element
// Current Matrix being pointed to by this and Matrix m must both be of same 
// rank and dims must be the same
    if(!samesize(m))
        throw MatrixLengthError();
    // Create a copy of current Matrix
    Matrix<T> mres=duplicate();   
    // uses gslice_array -=
    (*mres.data_ptr())[mres.getgslice()]*=(*m.data_ptr())[m.getgslice()]; 
    return MatrixRes<T>(mres);
}

// +++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator/
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator/
template<class U, class T>
MatrixRes<typename promote_trait<T,U>::T_promote> div_scalar(const Matrix<T>& m,const U& u)
{
// Scalar division
    typedef typename promote_trait<T,U>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());
    *mres.data_ptr()/=(T_promote)u;
    return MatrixRes<T_promote>(mres);
}

template<class T> template<class U>
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::operator/(const Scalar<U>& s) const
{ return div_scalar(*this,s.scalar());}
 
template<class T>
MatrixRes<T> Matrix<T>::operator/(const Matrix<T>& m) const
{
// --- Scalar division
// Create a copy of current matrix
   T t(m.scalar());
   return (*this)/t;
}

template<class T> template<class U>
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::operator/(const Matrix<U>& m) const
{
    Matrix<typename promote_trait<T,U>::T_promote> m1(*this);
    Matrix<typename promote_trait<T,U>::T_promote> m2(m);
    return m1/m2;
}

template<class T>
MatrixRes<typename promote_trait<T,int>::T_promote> operator/(const Matrix<T>& m, const int& i)
{
// --- Scalar division
// Create a copy of current matrix
//  
    typedef typename promote_trait<T,int>::T_promote T_promote;
    T_promote t=i;
    return div_scalar(m,t);
}

template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> operator/(const Matrix<T>& m, const double& d)
{
// --- Scalar division
// Create a copy of current matrix
//  
    typedef typename promote_trait<T,double>::T_promote T_promote;
    T_promote t=d;
    return div_scalar(m,t);
}

template<class T>
Scalar<typename promote_trait<T,int>::T_promote> operator/(const Scalar<T>& s, const int& i)
{
    typedef typename promote_trait<T,int>::T_promote T_promote;
    T_promote t=i;
    return Scalar<T_promote>(s.scalar()/t);
}

template<class T>
Scalar<typename promote_trait<T,double>::T_promote> operator/(const Scalar<T>& s, const double& d)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    T_promote t=d;
    return Scalar<T_promote>(s.scalar()/t);
}

template<class T>
Scalar<complex<T> > operator/(const Scalar<complex<T> >& s, const double& d)
{
  return Scalar<complex<T> >(s.scalar()/d);
}

template<class T>
MatrixRes<complex<T> > operator/(const Matrix<complex<T> >& m,const double& d)
{ return div_scalar(m,d);}

template<class T>
Scalar<complex<T> > operator/(const Scalar<complex<T> >& s, const int& i)
{
  T t=i;
  return Scalar<complex<T> >(s.scalar()/t);
}

template<class T>
MatrixRes<complex<T> > operator/(const Matrix<complex<T> >& m,const int& i)
{ return div_scalar(m,i);}
// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator*=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator*=

template<class T>
Matrix<T>& Matrix<T>::operator*=(const T& t)
{
// --- Matrix*=Scalar
   // uses gslice_array computed assignment...rhs must be a valarray
   valarray<T> v(t,size());
   (*_val)[getgslice()]*=v;
   return *this;
}

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

template<class T>
const Matrix<T>& Matrix<T>::operator*=(const T& t) const
{
// --- Matrix*=Scalar
   // uses gslice_array computed assignment...rhs must be a valarray
   valarray<T> v(t,size());
   (*_val)[getgslice()]*=v;
   return *this;
}

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

template<class T>
Matrix<T>& Matrix<T>::operator*=(const Matrix<T>& m)
{
// --- Matrix*=Matrix
   (*this)=(*this)*m;
   return *this;
}

template<class T> template<class U>
Matrix<T>& Matrix<T>::operator*=(const Matrix<U>& m)
{ Matrix<T> mu(m); return this->operator*=(mu);}

// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator/=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator/=

template<class T>
Matrix<T>& Matrix<T>::operator/=(const T& t)
{
// --- Matrix/=Scalar
   // uses gslice_array computed assignment...rhs must be a valarray
   valarray<T> v(t,size());
   (*_val)[getgslice()]/=v;
   return *this;
}

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

template<class T>
const Matrix<T>& Matrix<T>::operator/=(const T& t) const
{
// --- Matrix/=Scalar
   // uses gslice_array computed assignment...rhs must be a valarray
   valarray<T> v(t,size());
   (*_val)[getgslice()]/=v;
   return *this;
}

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

template<class T>
Matrix<T>& Matrix<T>::operator/=(const Matrix<T>& m)
{
// --- Matrix/=Matrix
   (*this)=(*this)/m;
   return *this;
}

template<class T>
const Matrix<T>& Matrix<T>::operator/=(const Matrix<T>& m) const
{
// --- Matrix/=Matrix
   (*this)=(*this)/m;
   return *this;
}

template<class T> template<class U>
Matrix<T>& Matrix<T>::operator/=(const Matrix<U>& m)
{ Matrix<T> mu(m); return this->operator/=(mu);}

template<class T> template<class U>
const Matrix<T>& Matrix<T>::operator/=(const Matrix<U>& m) const
{ Matrix<T> mu(m); return this->operator/=(mu);}


// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator==
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator==

template<class T> 
bool Matrix<T>::operator==(const Matrix<T>& m) const
{
    if (dimdata()!=m.dimdata())
        return false;
    valarray<bool> msk(false,size());
    transform(begin(),end(),m.begin(),&msk[0],not_equal_to<T>());
    if (msk.max())
        return false;
    else 
        return true;
}

// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::operator!=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator!=

template<class T>
bool Matrix<T>::operator!=(const Matrix<T>& m) const
{
    return !((*this)==m);
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::row
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- row

template<class T>
Matrix<T> Matrix<T>::row(const size_t i) const
// --- Access first dimension of Matrix
{
    return indexdim(1,i);
}

template<class T>
Matrix<T> Matrix<T>::row(const slice& sl) const
// --- Access first dimension of Matrix
{
    return indexdim(1,sl); 
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::col
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- col

template<class T>
Matrix<T> Matrix<T>::col(const size_t j) const
// --- Access second dimension of Matrix
{
     return indexdim(2,j);
}

template<class T>
Matrix<T> Matrix<T>::col(const slice& sl) const
// --- Access second dimension of Matrix
{
     return indexdim(2,sl);
}                  
                 
// +++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::layer
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- layer

template<class T>
Matrix<T> Matrix<T>::layer(const size_t k) const
// --- Access third dimension of Matrix
{
    return indexdim(3,k);
}

template<class T>
Matrix<T> Matrix<T>::layer(const slice& sl) const
// --- Access third dimension of Matrix
{
    return indexdim(3,sl);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::indexdim
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------- indexdim

template<class T>
Matrix<T> Matrix<T>::indexdim(const size_t idim,const size_t i) const
// --- Access nth dimension of Matrix
{
    // check that Matrix is of expected rank
    size_t dm=idim-1;
    if (dm<0 || dm>=rank())
        throw MatrixRankError();
    valarray<size_t> dms(dims());
    // check that index is in correct range
    if (i<0 || i>=dms[dm])
        throw MatrixIndexError();
    valarray<bool> mask(true,rank());
    mask[dm]=false;
    valarray<size_t> str(strides());
    size_t st=start()+i*str[dm];
    Dim_var dim=new Dim(st,dms[mask],str[mask],storageorder());                         
    return Matrix<T>(data_objvar(),dim,cdimno(),true,DEEP,!empty());
}

template<class T>
Matrix<T> Matrix<T>::indexdim(const size_t idim,const slice& sl) const
{
// --- Access nth dimension of Matrix object
    size_t dm=idim-1;
    // Check that matrix is of expected rank
    if (dm<0 || dm>=rank())
        throw MatrixRankError();
    valarray<size_t> dms(dims());
    int sz=dms[dm];
    if (sl.start()<0 || sl.start()>sz || sl.size()<0 || \
        sl.size()>sz)
            throw MatrixIndexError();
    // find starting position
    valarray<size_t> str(strides());
    size_t st=start()+sl.start()*str[dm];
    dms[dm]=sl.size();
    str[dm]=sl.stride()*str[dm];
    Dim_var dim=new Dim(st,dms,str,storageorder());                          
    return Matrix<T>(data_objvar(),dim,cdimno(),true,DEEP,!empty());
}
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::dot
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------------- dot
   
template<class T> template<class U>
inline typename promote_trait<T,U>::T_promote \
     Matrix<T>::dot(const Matrix<U>& m) const
{
// --- Calculate dot product of two rank 1 matrices
     if (rank()!=1 || m.rank()!=1)
         throw MatrixRankError();
     if (dims()[0]!=m.dims()[0])
         throw MatrixLengthError();
     typedef typename promote_trait<T,U>::T_promote T_promote;
     T_promote sum=0;
     sum=inner_product(begin(),end(),m.begin(),sum);
     return sum;
     
}

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)
{
    typedef typename promote_trait<T,U>::T_promote T_promote;
    int sz=ms.size();
    if(sz!=v.size())
        throw MatrixLengthError();
    if (!ms.is_homogeneous_seq())
        throw MatrixDomainError();
    if (ms.elt(0).rank()!=1 )
        throw MatrixRankError();   
    MatrixSeq<T_promote> msres(sz);
    
    for(register int i=0;i<ms.size();i++)
        msres[i]=ms[i].dot(v[i]);
    return MatrixResSeq<T_promote>(msres);   
}

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)
{
    typedef typename promote_trait<T,U>::T_promote T_promote;
    int sz=ms.size();
    if (!ms.is_homogeneous_seq())
        throw MatrixDomainError();
    if (ms.elt(0).rank()!=1 )
        throw MatrixRankError();
    if(m.rank()!=1)
        throw MatrixRankError();   
    MatrixSeq<T_promote> msres(sz);
    
    for(register int i=0;i<ms.size();i++)
        msres[i]=ms[i].dot(m);
    return MatrixResSeq<T_promote>(msres);   
}

template<class U, class T, class ContainerT >
MatrixResSeq<typename promote_trait<T,U>::T_promote> dot(const Matrix<U>& m,\
    MatrixBaseSeq<T,ContainerT >& ms)
{ return dot(ms,m);}

// +++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::cross
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- cross

template<class T> template<class U>
MatrixRes<typename promote_trait<T,U>::T_promote> Matrix<T>::cross(const Matrix<U>& m) const
{
// --- Calculate vector cross product of two rank 1 matrices
    if(rank()!=1 || m.rank()!=1)
       throw MatrixRankError();
    size_t n=dims()[0];
    if (n!=m.dims()[0])
        throw MatrixLengthError();
    if (n<1 || n>3) 
        throw MatrixDomainError();
    typedef typename promote_trait<T,U>::T_promote T_promote;
    Dim d(n);
    Matrix<T_promote> m1(d);
    Matrix<T_promote> m2(d);
    transform(begin(),end(),m.begin()+1,m1.begin()+(n-1),\
        multiplies<T_promote>());
    transform(m.begin(),m.end(),begin()+1,m2.begin()+(n-1),\
        multiplies<T_promote>());
    m1-=m2;
    return MatrixRes<T_promote>(m1);
}

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)
{
    typedef typename promote_trait<T,U>::T_promote T_promote;
    int sz=ms.size();
    if(sz!=v.size())
        throw MatrixLengthError();
    if (!ms.is_homogeneous_seq())
        throw MatrixDomainError();
    if (ms.elt(0).rank()!=1 )
        throw MatrixRankError();   
    MatrixSeq<T_promote> msres(sz);
    
    for(register int i=0;i<ms.size();i++)
        msres[i]=ms[i].cross(v[i]);
    return MatrixResSeq<T_promote>(msres);   
}

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)
{
    typedef typename promote_trait<T,U>::T_promote T_promote;
    int sz=ms.size();
    if (!ms.is_homogeneous_seq())
        throw MatrixDomainError();
    if (ms.elt(0).rank()!=1 )
        throw MatrixRankError();
    if(m.rank()!=1)
        throw MatrixRankError();   
    MatrixSeq<T_promote> msres(sz);
    
    for(register int i=0;i<ms.size();i++)
        msres[i]=ms[i].cross(m);
    return MatrixResSeq<T_promote>(msres);   
}

template<class U, class T, class ContainerT >
MatrixResSeq<typename promote_trait<T,U>::T_promote> cross(const Matrix<U>& m,\
    MatrixBaseSeq<T,ContainerT >& ms)
{
    typedef typename promote_trait<T,U>::T_promote T_promote;
    int sz=ms.size();
    if (!ms.is_homogeneous_seq())
        throw MatrixDomainError();
    if (ms.elt(0).rank()!=1 )
        throw MatrixRankError();
    if(m.rank()!=1)
        throw MatrixRankError();   
    MatrixSeq<T_promote> msres(sz);
    
    for(register int i=0;i<ms.size();i++)
        msres[i]=m.cross(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}



// +++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::matmult
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- matmult

template<class T>    
void Matrix<T>::matmult(const Matrix<T>& m,Matrix<T>& mres) const
{
// --- Generic matrix multiplication routine
// mres is for storing result in.
    // check that there are no rank 0 matrices (i.e scalars)
    if(rank()==0 || m.rank()==0)
        throw MatrixRankError();
    register int i,j;
   
    if (rank()==1)
    {
       // current Matrix object is of rank 1 
        if (m.rank()==1)
           mres=dot(m);  // vector product if both matrices are rank 1
        else 
        {
            int msz=m.dims()[m.rank()-1];
            for(j=0;j<msz;j++)
            {
                Matrix<T> mref=mres.indexdim(mres.rank(),j).duplicate(SHALLOW);
                matmult(m.indexdim(m.rank(),j),mref);
            }
        }
     }
     else
     {
         if (m.rank()==1)
         {
             int sz=dims()[0];
             for(i=0;i<sz;i++)   
             {
                 Matrix<T> mref=mres.row(i).duplicate(SHALLOW);
                 row(i).matmult(m,mref);
             }
          }           
          else
          {
              int sz=dims()[0];
              int msz=m.dims()[m.rank()-1];
              for(i=0;i<sz;i++)
                  for(j=0;j<msz;j++)
                  {
                      Matrix<T> mref=mres.row(i).indexdim(mres.rank()-1,j).\
                          duplicate(SHALLOW);
                      row(i).matmult(m.indexdim(m.rank(),j),\
                          mref);
                  }
          }
      }            
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::solve
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- solve    
template<class T>
inline MatrixRes<typename promote_trait<T,double>::T_promote>\
Matrix<T>::solve(const Matrix<T>& B) const
{
   typedef typename promote_trait<T,double>::T_promote T_promote;
   // scalar case ..trivial
   if (rank()==0)
   {
      if (B.rank()!=0)
          throw MatrixRankError();
      Matrix<T_promote> denom(*this);
      return (B/denom);
   }
   else
   { 
       // ... ensure B is of correct rank
       if (B.rank()!=rank()-1 && B.rank()!=rank())
           throw MatrixRankError();
       // A is Rank 2
       if (rank()==2)
       {    
            MatrixSVD<T> svd(*this);
            return (svd.solve(B));   
       }  
       else
       { 
           Matrix<T> A,B2;
           A.clone(*this);    // create a clone of current matrix
           B2.clone(B); // create a clone of B matrix i.e. pt to same data
           if (B2.rank()==A.rank()) // make B rank 1 or 2
               B2.changerank(2);
           else
               B2.changerank(1);
           A.changerank(2);  // make A rank 2   
           return (A.solve(B2)); 
       }
    }
}
// ++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::pinverse
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------- pinverse 
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> Scalar<T>::pinverse() const  
{
   typedef typename promote_trait<T,double>::T_promote T_promote;
   T_promote I=1;
   return Scalar<T_promote>(I/scalar());
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> Matrix<T>::pinverse() const  
{
// calculate the pseudo-inverse of a rectangular matrix 
   typedef typename promote_trait<T,double>::T_promote T_promote;
   size_t i;
   Matrix<T_promote> I;
   valarray<size_t> dms(dims());
   int sz;
                        
   switch(rank())
   {
       case 0: // trivial if scalars
           I.assign_val(T_promote(1));
           return (I/this->scalar());
           
       case 1: // 1/norm(*this) for rank 1 matrices
           I.assign_val(T_promote(1));
           return (I/norm(*this));
               
       case 3:
       // must be m x m x n or m x m x m
           if(dms[0]!=dms[1])
               throw MatrixDomainError();
       case 2: 
       // m x n or m x m 
           sz=dms[0]; 
           I.resize(Dim(sz,sz));
           I.unit_matrix();
           return solve(I);
                        
       default:
       // cannot invert other matrices...
           throw MatrixDomainError();
               
   }
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::resize
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- resize
template<class T>
void Matrix<T>::resize(const Dim& d)
{
// resize current matrix object (providing it is not
// a subarray)
    if (_subarray)
        throw MatrixDomainError(); // not subarrays
    // resize dimension info
    _diminfo=new Dim(d.dims(),d.storageorder());
    // reinitialise space on stack (Var<T> operator= in initailise
    // will handle freeing of memory
    initialise();
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::reshape
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- reshape     
template<class T>
void Matrix<T>::reshape(const Dim& d)
{
// NOT resize i.e. size() remains unchanged...just dimension
// object changes
   if (size()!=d.size())
       throw MatrixLengthError();
   _diminfo=new Dim(d);

} 
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pow 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- pow
template<class T>
Scalar<typename promote_trait<T,double>::T_promote> pow(const Scalar<T>& s, const T& t)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    return Scalar<T_promote>(std::pow(s.scalar(),t));
}

template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const Matrix<T>& m, \
  const Scalar<T>& s)
{ return pow(m,s.scalar());}
  
template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const Matrix<T>& m, const T& t)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());
    (*mres.data_ptr())=pow((*mres.data_ptr()),(T_promote)t);
    return MatrixRes<T_promote>(mres);
}

template<class T>
MatrixRes<complex<T> > pow(const Matrix<complex<T> >& m, const complex<T>& t)
{
    Matrix<complex<T> > mres(m.duplicate());
    for(register int i=0;i<m.size();i++)
        (*mres.data_ptr())[i]=std::pow((*mres.data_ptr())[i],t);
    return MatrixRes<complex<T> >(mres);
}

template<class T>
Scalar<typename promote_trait<T,double>::T_promote> pow(const T& t, \
  const Scalar<T>& s)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    return Scalar<T_promote>(std::pow(t,s.scalar()));
}

template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const Scalar<T>& s, \
  const Matrix<T>& m)
{ return pow(s.scalar(),m);}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const T& t, \
    const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());
    (*mres.data_ptr())=pow((T_promote)t,(*mres.data_ptr()));
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > pow(const complex<T>& t, \
    const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());
    for(register int i=0;i<m.size();i++)
        (*mres.data_ptr())[i]=std::pow(t,(*mres.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

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)
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;
    T_promote arg1(s1.scalar());
    T_promote arg2(s2.scalar());
    return Scalar<T_promote>(std::pow(arg1,arg2));
}

template<class T>
MatrixRes<typename promote_trait<T,double>::T_promote> pow(const Matrix<T>& m1, \
    const Matrix<T>& m2)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    if(m1.rank()==0)
        return pow(m1.scalar(),m2);
    if(m2.rank()==0)
        return pow(m1,m2.scalar());
    if(m1.dimdata()!=m2.dimdata())
        throw MatrixLengthError();    
    Matrix<T_promote> m_arg1(m1.duplicate());
    Matrix<T_promote> m_arg2(m2.duplicate());
    (*m_arg1.data_ptr())=pow((*m_arg1.data_ptr()),(*m_arg2.data_ptr()));
    return MatrixRes<T_promote>(m_arg1);
}

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)
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;
    Matrix<T_promote> m2(m.duplicate());
    Scalar<T_promote> s2(s.scalar());
    return pow(s2,m2);
}

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)
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;
    Matrix<T_promote> m2(m.duplicate());
    Scalar<T_promote> s2(s.scalar());
    return pow(m2,s2);
}

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)
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;
    Matrix<T_promote> m1_arg(m1.duplicate());
    Matrix<T_promote> m2_arg(m2.duplicate());
    return pow(m1_arg,m2_arg);
}

template<class T>
MatrixRes<complex<T> > pow(const Matrix<complex<T> >& m1, \
     const Matrix<complex<T> >& m2)
{
    if(m1.rank()==0)
        return pow(m1.scalar(),m2);
    if(m2.rank()==0)
        return pow(m1,m2.scalar());
    if(m1.dimdata()!=m2.dimdata())
        throw MatrixLengthError(); 
    Matrix<complex<T> > mres(m1.duplicate());
    const valarray<complex<T> >& v2((*m2.data_ptr())[m2.getgslice()]);
    for(register int i=0;i<m1.size();i++)
        (*mres.data_ptr())[i]=std::pow((*mres.data_ptr())[i],v2[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT >  
MatrixResSeq<T> pow(const MatrixBaseSeq<T, ContainerT >& ms, const T& t) 
{
    MatrixSeq<T> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=pow(ms[i],t);
    return MatrixResSeq<T>(msres);   
}

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) 
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=pow(ms[i],s);
    return MatrixResSeq<T_promote>(msres);   
}

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) 
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=pow(ms[i],m);
    return MatrixResSeq<T_promote>(msres);   
}

template<class T, class ContainerT >   
MatrixResSeq<T> pow(const T& t, const MatrixBaseSeq<T, ContainerT >& ms)
{
    MatrixSeq<T> msres(ms.size());
      for(register int i=0;i<ms.size();i++)
        msres[i]=pow(t,ms[i]);
    return MatrixResSeq<T>(msres);   
}

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)
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;   
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=pow(s,ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}

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)
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;   
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=pow(m,ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}

template<class T, class ContainerT >  
MatrixResSeq<T> pow(const MatrixBaseSeq<T, ContainerT >& ms, \
    const vector<T>& v) 
{
    int sz=ms.size();
    if(sz!=v.size())
        throw MatrixLengthError();
    MatrixSeq<T> msres(sz);
    for(register int i=0;i<sz;i++)
        msres[i]=pow(ms[i],v[i]);
    return MatrixResSeq<T>(msres);   
}

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) 
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;   
    int sz=ms.size();
    if(sz!=v.size())
        throw MatrixLengthError();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<sz;i++)
        msres[i]=pow(ms[i],v[i]);
    return MatrixResSeq<T_promote>(msres);   
}

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) 
{
    typedef typename promote_trait<typename promote_trait<T,U>::T_promote,double>::T_promote \
        T_promote;   
    int sz=ms.size();
    if(sz!=v.size())
        throw MatrixLengthError();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<sz;i++)
        msres[i]=pow(ms[i],v[i]);
    return MatrixResSeq<T_promote>(msres);   
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::power
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- power                 
template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> \
Matrix<T>::power(const double& n) const
{  
     typedef typename promote_trait<T,double>::T_promote T_promote;   
     T_promote p(n);
     T_promote halfp(n/2);
     if (rank()==0)
         return MatrixRes<T_promote>(std::pow((T_promote)scalar(),p));
     if(rank()==1)   
         return MatrixRes<T_promote>(std::pow(dot(*this),halfp));
     // verify that matrix is square
     if(!dimvar()->is_square())
         throw MatrixDomainError();
     Matrix<T_promote> m;
     int ip=(int)n;
     if(ip<0)
        return(power(-ip)).pinverse();
        
     if (ip==0)
     {
         int sz=dims()[0];
         m.resize(Dim(sz,sz));
         m.unit_matrix();
         return MatrixRes<T_promote>(m);
     }
     if (ip % 2 !=0)
     {
         m=(*this)*power(ip-1);
         return MatrixRes<T_promote>(m);
     }
     m=power(ip/2);
     return (m*m);
}

// +++++++++++++++++++++++++++++++++++++++++++++ class Matrix::eigenvalues
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------- eigenvalues
template<class T> 
ScalarSeq<double> eigenvalues(const Matrix<T>& m)
{
    MatrixEig eig(m);
    return eig.eigenvalues(); 
}

// ++++++++++++++++++++++++++++++++++++++++++++ class Matrix::eigenvectors
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ---------------------------------------------------------- eigenvectors     
template<class T> 
MatrixResSeq<double> eigenvectors(const Matrix<T>& m)
{
    MatrixEig eig(m);  
    return MatrixResSeq<double>(eig.eigenvectors()); 
}

// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::plusreduce
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ plusreduce
template<class T> 
MatrixRes<T> Matrix<T>::plusreduce(const size_t idim) const
{
   if (rank()==0)
       return MatrixRes<T>(*this);
   else
   {
       int sz=rank();
       valarray<size_t> d(dims());
       if (idim<=0 || idim>sz)
           throw MatrixIndexError();
       valarray<bool> msk(true,sz);
       msk[idim-1]=false; 
       valarray<size_t> dm(d[msk]);
       Matrix<T> mres(dm);
       Matrix<T> mref;
       if (sz==1)
       {     
           valarray<T> v((*_val)[getgslice()]);
           mres.assign(v.sum());
       }
       else
       {         
           size_t idx=1;
           size_t idim_2=idim-1;
           if(idim==1)
           {
               idx=2;
               idim_2=1;
           }
           register size_t i;
           for(i=0;i<d[idx-1];i++)
           {
                mref=mres.indexdim(1,i).duplicate(SHALLOW);
                indexdim(idx,i).plusreduce(idim_2,mref);
           }
       }        
       return MatrixRes<T>(mres);
   }
}

template<class T>
void Matrix<T>::plusreduce(const size_t idim,Matrix<T>& mres) const
{
   // create local temporary copy of mres (pointing to same data)
   Matrix<T> mref;          
    
         
   if (rank()==1)
   {
       if (idim!=1)
           throw MatrixIndexError();
       else
       {
           valarray<T> v((*_val)[getgslice()]);
           mres.assign(v.sum());
       }
   }
   else
   {
       size_t idx=1;
       size_t idim_2=idim-1;
       if(idim==1)
       {
           idx=2;
           idim_2=1;
       }
       register size_t i;
       int sz=dims()[idx-1];
       for(i=0;i<sz;i++)
       {
           mref=mres.indexdim(1,i).duplicate(SHALLOW);
           indexdim(idx,i).plusreduce(idim_2,mref);
       }    
   }
}
// ++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::changerank
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ changerank     
template<class T>   
void Matrix<T>::changerank(const size_t rk)
{
 // change rank of current matrix
   
    size_t *dm=new size_t[rk];
    register size_t i;
    DimStorageOption order=storageorder();
     
    if (rk!=rank())  // any change?
    {
        if (rk>rank()) // uprank
        {
            for (i=0;i<rk-rank();i++)
                dm[i]=1;
            for (i=rk-rank();i<rk;i++)
                dm[i]=dims()[i-rk+rank()];         
       
        }
        else          // downrank
        {
            dm[0]=1;
            for (i=0;i<rank()-rk+1;i++)
                dm[0]*=dims()[i];
            for (i=rank()-rk+1;i<rank();i++)
                dm[i-rank()+rk]=dims()[i];
        }         
         
        reshape(Dim(start(),setdims(dm,rk),computestrides(dm,rk,order),order));
        delete [] dm;
    }
}

// +++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::display
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- display
     
template<class T>        
Void Matrix<T>::display() const
{
// --- Display function for Matrix class object
// Output each element in matrix. Print a newline character between
// each dimension
// _strides[0] contains no. of elements between each no. in Matrix (not
// necessarily one)
    Void sink;
    cout << "Matrix Object, ";

    dimvar()->display();
    
    iterator itr=begin();
    size_t pos=itr.eltno(); 
    valarray<size_t> str(computestrides(dims()));
   
    register int i;
    
    cout << endl;
       
    while(itr<end())
    {
        if (rank()>0)
        {
         // for each dimension (except lowest) see if newline character reqd
            for (i=rank()-1;i>0;i--)
            {
                if ( pos!=0 && pos % str[i-1]==0)
                    cout << endl;
            }  
        }      
        cout << *itr << "\t";
        itr++;
        pos=itr.eltno();
               
    } 
    cout << endl;
    cout << endl; 
    return sink; 
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++ class Matrix::atan2
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- atan2
template<class T> Scalar<typename promote_trait<T,double>::T_promote> \
atan2(const Scalar<T>& s1, const Scalar<T>& s2)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Scalar<T_promote> s3(s1.scalar());
    Scalar<T_promote> s4(s2.scalar());
    return Scalar<T_promote>(std::atan2(s3,s4));
}

template<class T> MatrixRes<typename promote_trait<T,double>::T_promote> \
atan2(const Matrix<T>& m, const Scalar<T>& s)
{ return atan2(m,s.scalar());}

template<class T> MatrixRes<typename promote_trait<T,double>::T_promote> \
atan2(const Matrix<T>& m, const T& t)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> m_arg1(m.duplicate());
    (*m_arg1.data_ptr())=std::atan2((*m_arg1.data_ptr()),(T_promote)t);
    return MatrixRes<T_promote>(m_arg1);
}

template<class T> MatrixRes<typename promote_trait<T,double>::T_promote> \
atan2(const Scalar<T>& s, const Matrix<T>& m)
{ return atan2(s.scalar(),m);}

template<class T> MatrixRes<typename promote_trait<T,double>::T_promote> \
atan2(const T& t, const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> m_arg1(m.duplicate());
    (*m_arg1.data_ptr())=atan2((T_promote)t,(*m_arg1.data_ptr()));
    return MatrixRes<T_promote>(m_arg1);
}

template<class T> MatrixRes<typename promote_trait<T,double>::T_promote> \
atan2(const Matrix<T>& m1, const Matrix<T>& m2)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    if(m1.rank()==0)
        return atan2(m1.scalar(),m2);
    if(m2.rank()==0)
        return atan2(m1,m2.scalar());
    if(m1.dimdata()!=m2.dimdata())
        throw MatrixLengthError();
    Matrix<T_promote> m_arg1(m1.duplicate());
    Matrix<T_promote> m_arg2(m2.duplicate());
    (*m_arg1.data_ptr())=atan2((*m_arg1.data_ptr()),(*m_arg2.data_ptr()));
    return MatrixRes<T_promote>(m_arg1);
}

template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const \
    MatrixBaseSeq<T, ContainerT>& ms,const T& t)
{   
    typedef typename promote_trait<T,double>::T_promote T_promote;
    int sz=ms.size();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<ms.size();i++)
        msres[i]=atan2(ms[i],t);
    return MatrixResSeq<T_promote>(msres);   
}

template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const \
    MatrixBaseSeq<T, ContainerT>& ms,const Scalar<T>& s)
{   
    typedef typename promote_trait<T,double>::T_promote T_promote;
    int sz=ms.size();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<ms.size();i++)
        msres[i]=atan2(ms[i],s);
    return MatrixResSeq<T_promote>(msres);   
}

template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const \
    MatrixBaseSeq<T, ContainerT>& ms,const Matrix<T>& m)
{   
    typedef typename promote_trait<T,double>::T_promote T_promote;
    int sz=ms.size();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<ms.size();i++)
        msres[i]=atan2(ms[i],m);
    return MatrixResSeq<T_promote>(msres);   
}

template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const T& t,\
    const MatrixBaseSeq<T, ContainerT>& ms)
{   
    typedef typename promote_trait<T,double>::T_promote T_promote;
    int sz=ms.size();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<ms.size();i++)
        msres[i]=atan2(t,ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}

template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const Scalar<T>& s,\
    const MatrixBaseSeq<T, ContainerT>& ms)
{   
    typedef typename promote_trait<T,double>::T_promote T_promote;
    int sz=ms.size();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<ms.size();i++)
        msres[i]=atan2(s,ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}

template<class T, class ContainerT>
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan2(const Matrix<T>& m,\
    const MatrixBaseSeq<T, ContainerT>& ms)
{   
    typedef typename promote_trait<T,double>::T_promote T_promote;
    int sz=ms.size();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<ms.size();i++)
        msres[i]=atan2(m,ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}

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)
{   
    typedef typename promote_trait<T,double>::T_promote T_promote;
    int sz=ms.size();
    if(sz!=v.size())
        throw MatrixLengthError();
    MatrixSeq<T_promote> msres(sz);
    for(register int i=0;i<ms.size();i++)
        msres[i]=atan2(ms[i],v[i]);
    return MatrixResSeq<T_promote>(msres);   
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ hilbert
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- hilbert 
inline MatrixRes<double> hilbert(int n)
{
// return a square hilbert matrix of size n x n where the (i,j)th element
// (0 <= i < n), (0 < =j < n) is 1/(i+j+1)
   if (n<=0)
       throw MatrixDomainError();
   Matrix<double> H(Dim(n,n));
   int i,j;
   for (i=0;i<n;i++)
       for (j=0;j<n;j++)
           H(i,j)=(double)1.0/(i+j+1);
   return MatrixRes<double>(H);
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ cos
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- cos 
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> cos(const Scalar<T>& s)
{
  typedef typename promote_trait<T,double>::T_promote T_promote;
  return Scalar<T_promote>(std::cos(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> cos(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=cos(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > cos(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());   
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=std::cos((*mres.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> cos(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=cos(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
  

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ sin 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- sin  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> sin(const Scalar<T>& s)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    return Scalar<T_promote>(std::sin(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> sin(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=sin(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > sin(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());   
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=std::sin((*mres.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> sin(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=sin(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
  

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ tan 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- tan  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> tan(const Scalar<T>& s)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    return Scalar<T_promote>(std::tan(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> tan(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=tan(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > tan(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());   
    for(register int i=0;i<mres.size();i++)
    {
        (*mres.data_ptr())[i]=std::sin((*mres.data_ptr())[i]);
        (*mres.data_ptr())[i]/=std::cos((*mres.data_ptr())[i]);
    }
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> tan(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=tan(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
  

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ sqrt 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ sqrt  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> sqrt(const Scalar<T>& s)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    return Scalar<T_promote>(std::sqrt(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> sqrt(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=sqrt(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > sqrt(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());   
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=std::sqrt((*mres.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> sqrt(\
    const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote; 
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=sqrt(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ log 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- log  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> log(const Scalar<T>& s)
{
  typedef typename promote_trait<T,double>::T_promote T_promote;
  return Scalar<T_promote>(std::log(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> log(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=log(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > log(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());   
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=std::log((*mres.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> log(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=log(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ log10 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- log10  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> log10(const Scalar<T>& s)
{
  typedef typename promote_trait<T,double>::T_promote T_promote;
  return Scalar<T_promote>(std::log10(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> log10(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=log10(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
Scalar<complex<T> > log10(const Scalar<complex<T> >& s)
{
  T log10=std::log(10.);
  Scalar<complex<T> > res(std::log(s.scalar()));
  return Scalar<complex<T> >(res.scalar()/log10);
}

template<class T> 
MatrixRes<complex<T> > log10(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate()); 
    T log10=std::log(10.);  
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=std::log((*mres.data_ptr())[i]);
    (*mres.data_ptr())/=log10;
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> log10(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=log10(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ exp 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- exp  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> exp(const Scalar<T>& s)
{
  typedef typename promote_trait<T,double>::T_promote T_promote;
  return Scalar<T_promote>(std::exp(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> exp(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=exp(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > exp(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());   
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=std::exp((*mres.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> exp(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=exp(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pinverse
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- pinverse 

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> pinverse(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
         msres[i]=ms[i].pinverse();
    return MatrixResSeq<T_promote>(msres);   
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ acos 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ acos 

template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> acos(const Scalar<T>& s)
{
  typedef typename promote_trait<T,double>::T_promote T_promote;
  return Scalar<T_promote>(std::acos(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> acos(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=acos(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> acos(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=acos(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
  

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ asin 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ asin
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> asin(const Scalar<T>& s)
{ return Scalar<T>(std::asin(s.scalar()));}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> asin(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=asin(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> asin(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=asin(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
  
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ atan 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ atan  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> atan(const Scalar<T>& s)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    return Scalar<T_promote>(std::atan(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> atan(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=std::atan(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> atan(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=atan(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
} 
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ cosh 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ cosh  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> cosh(const Scalar<T>& s)
{
  typedef typename promote_trait<T,double>::T_promote T_promote;
  return Scalar<T_promote>(std::cosh(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> cosh(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=cosh(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > cosh(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());   
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=std::cosh((*mres.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> cosh(const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=cosh(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
  

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ sinh 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ sinh 
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> sinh(const Scalar<T>& s)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    return Scalar<T_promote>(std::sinh(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> sinh(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m.duplicate());   
    (*mres.data_ptr())=sinh(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > sinh(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());   
    for(register int i=0;i<mres.size();i++)
        (*mres.data_ptr())[i]=std::sinh((*mres.data_ptr())[i]);
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> sinh(\
    const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote; 
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=sinh(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ tanh 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ tanh  
template<class T> 
Scalar<typename promote_trait<T,double>::T_promote> tanh(const Scalar<T>& s)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    T_promote res=std::sinh(s.scalar());
    return Scalar<T_promote>(res/std::cosh(s.scalar()));
}

template<class T> 
MatrixRes<typename promote_trait<T,double>::T_promote> tanh(const Matrix<T>& m)
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    Matrix<T_promote> mres(m);
    (*mres.data_ptr())=sinh(*mres.data_ptr());
    (*mres.data_ptr())/=cosh(*mres.data_ptr());
    return MatrixRes<T_promote>(mres);
}

template<class T> 
MatrixRes<complex<T> > tanh(const Matrix<complex<T> >& m)
{
    Matrix<complex<T> > mres(m.duplicate());
    for(register int i=0;i<mres.size();i++)
    {
        (*mres.data_ptr())[i]=std::sinh((*mres.data_ptr())[i]);
        (*mres.data_ptr())[i]/=std::cosh((*mres.data_ptr())[i]);
    }
    return MatrixRes<complex<T> >(mres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> tanh(\
    const MatrixBaseSeq<T, ContainerT >& ms)
{
    typedef typename promote_trait<T,double>::T_promote T_promote; 
    MatrixSeq<T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=tanh(ms[i]);
    return MatrixResSeq<T_promote>(msres);   
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ max 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- max  
template<class T> 
T max(const Matrix<T>& m)
{ return m.data_ptr()->max();}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ min 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- min  
template<class T> 
T min(const Matrix<T>& m)
{ return m.data_ptr()->min();}
// ++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator*
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator*     
template<class T, class ContainerT >   
MatrixResSeq<T> MatrixBaseSeq<T, ContainerT >::operator*(const T& t) const
{
    MatrixSeq<T> msres(*this);
    msres*=t;
    return MatrixResSeq<T>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator*\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<int>& v) 
{
    MatrixSeq<typename promote_trait<T,int>::T_promote> msres(ms);
    msres*=v;
    return MatrixResSeq<typename promote_trait<T,int>::T_promote>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator*\
    (const vector<int>& v, const MatrixBaseSeq<T, ContainerT >& ms) 
{ return ms*v;}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator*\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<double>& v) 
{
    MatrixSeq<typename promote_trait<T,double>::T_promote> msres(ms);
    msres*=v;
    return MatrixResSeq<typename promote_trait<T,double>::T_promote>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator*\
    (const vector<double>& v, const MatrixBaseSeq<T, ContainerT >& ms) 
{ return ms*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) 
{
    MatrixSeq<typename promote_trait<T,complex<U> >::T_promote> msres(ms);
    msres*=v;
    return MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote>(msres);
}

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) 
{ return ms*v;}

template<class T, class ContainerT > template<class U>  
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator*(const Scalar<U>& s) const 
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]*s;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U>  
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator*(const Matrix<U>& m) const 
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]*m;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator*(const vector<Matrix<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]*ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator*(const vector<Scalar<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]*ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator*(const vector<MatrixRes<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]*ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator*\
    (const int& i,const MatrixBaseSeq<T, ContainerT >& ms)
{
    return ms*i;
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator*\
    (const double& d,const MatrixBaseSeq<T, ContainerT >& ms)
{
    return ms*d;
}

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)
{
    return ms*c;
}   
   
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)
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=s*ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

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)
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=u*ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

// +++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator*=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator*=    
template<class T, class ContainerT >
MatrixBaseSeq<T,ContainerT >& operator*=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<int>& v)
{ ms.mul_assign(v); return ms;}

template<class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator*=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<double>& v)
{ ms.mul_assign(v); return ms;}

template<class U, class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator*=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<complex<U> >& v)
{ ms.mul_assign(v); return ms;}

// ++++++++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator^
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------- operator^    
template<class T, class ContainerT >  
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator^(const double& t) const
{
    typedef typename promote_trait<T,double>::T_promote T_promote;   
    MatrixSeq<T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i].power(t);
    
    return MatrixResSeq<T_promote>(msres);
}    
    
template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,double>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator^(const vector<U>& v) const
{
    typedef typename promote_trait<T,double>::T_promote T_promote;
    MatrixSeq<T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i].power((double)v[i]);
    
    return MatrixResSeq<T_promote>(msres);
}
// ++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator+
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator+     
template<class T, class ContainerT >   
MatrixResSeq<T> MatrixBaseSeq<T, ContainerT >::operator+(const T& t) const
{
    MatrixSeq<T> msres(*this);
    msres+=t;
    return MatrixResSeq<T>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator+\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<int>& v) 
{
    MatrixSeq<typename promote_trait<T,int>::T_promote> msres(ms);
    msres+=v;
    return MatrixResSeq<typename promote_trait<T,int>::T_promote>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator+\
    (const vector<int>& v, const MatrixBaseSeq<T, ContainerT >& ms) 
{ return ms+v;}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator+\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<double>& v) 
{
    MatrixSeq<typename promote_trait<T,double>::T_promote> msres(ms);
    msres+=v;
    return MatrixResSeq<typename promote_trait<T,double>::T_promote>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator+\
    (const vector<double>& v, const MatrixBaseSeq<T, ContainerT >& ms) 
{ return ms+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) 
{
    MatrixSeq<typename promote_trait<T,complex<U> >::T_promote> msres(ms);
    msres+=v;
    return MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote>(msres);
}

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) 
{ return ms+v;}

template<class T, class ContainerT > template<class U>  
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator+(const Scalar<U>& s) const 
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]+s;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U>  
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator+(const Matrix<U>& m) const 
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]+m;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator+(const vector<Matrix<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]+ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator+(const vector<Scalar<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]+ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator+(const vector<MatrixRes<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]+ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator+\
    (const int& i,const MatrixBaseSeq<T, ContainerT >& ms)
{
    return ms+i;
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator+\
    (const double& d,const MatrixBaseSeq<T, ContainerT >& ms)
{
    return ms+d;
}

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)
{
    return ms+c;
}
   
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)
{
    return ms+s;
}

   
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)
{
    return ms+u;
}

// +++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator+=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator+=    
template<class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator+=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<int>& v)
{ ms.add_assign(v); return ms;}

template<class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator+=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<double>& v)
{ ms.add_assign(v); return ms;}

template<class U, class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator+=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<complex<U> >& v)
{ ms.add_assign(v); return ms;}

// ++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator-
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator-
template<class T, class ContainerT >   
MatrixResSeq<T> MatrixBaseSeq<T, ContainerT >::operator-() const
{
    int sz=this->size();
    MatrixSeq<T> msres(sz);
    for(register int i=0;i<sz;i++)
        msres[i]=-this->elt(i);
    return MatrixResSeq<T>(msres);  
}

template<class T, class ContainerT >   
MatrixResSeq<T> MatrixBaseSeq<T, ContainerT >::operator-(const T& t) const
{
    MatrixSeq<T> msres(*this);
    msres-=t;
    return MatrixResSeq<T>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator-\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<int>& v) 
{
    MatrixSeq<typename promote_trait<T,int>::T_promote> msres(ms);
    msres-=v;
    return MatrixResSeq<typename promote_trait<T,int>::T_promote>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator-\
    (const vector<int>& v, const MatrixBaseSeq<T, ContainerT >& ms) 
{ 
     MatrixSeq<typename promote_trait<T,int>::T_promote> msres(ms.size());
     for(register int i=0;i<ms.size();i++)
         msres[i]=(-ms[i])+v[i];
     return MatrixResSeq<typename promote_trait<T,int>::T_promote>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator-\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<double>& v) 
{
    MatrixSeq<typename promote_trait<T,double>::T_promote> msres(ms);
    msres-=v;
    return MatrixResSeq<typename promote_trait<T,double>::T_promote>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator-\
    (const vector<double>& v, const MatrixBaseSeq<T, ContainerT >& ms) 
{ 
     MatrixSeq<typename promote_trait<T,double>::T_promote> msres(ms.size());
     for(register int i=0;i<ms.size();i++)
         msres[i]=(-ms[i])+v[i];
     return MatrixResSeq<typename promote_trait<T,double>::T_promote>(msres);
}

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) 
{
    MatrixSeq<typename promote_trait<T,complex<U> >::T_promote> msres(ms);
    msres-=v;
    return MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote>(msres);
}

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) 
{ 
     MatrixSeq<typename promote_trait<T,complex<U> >::T_promote> msres(ms.size());
     for(register int i=0;i<ms.size();i++)
         msres[i]=(-ms[i])+v[i];
     return MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote>(msres);
}

template<class T, class ContainerT > template<class U>  
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator-(const Scalar<U>& s) const 
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]-s;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U>  
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator-(const Matrix<U>& m) const 
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]-m;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator-(const vector<Matrix<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]-ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator-(const vector<Scalar<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=(*this)[i]-ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator-(const vector<MatrixRes<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]-ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator-\
    (const int& iv,const MatrixBaseSeq<T, ContainerT >& ms)
{
     MatrixSeq<typename promote_trait<T,int>::T_promote> msres(ms.size());
     for(register int i=0;i<ms.size();i++)
         msres[i]=(-ms[i])+iv;
         
     return MatrixResSeq<typename promote_trait<T,int>::T_promote>(msres);
}

template<class T, class ContainerT > 
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator-\
    (const double& d,const MatrixBaseSeq<T, ContainerT >& ms)
{
    MatrixSeq<typename promote_trait<T,double>::T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
         msres[i]=(-ms[i])+d;
    return MatrixResSeq<typename promote_trait<T,double>::T_promote>(msres);
}

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)
{
    MatrixSeq<typename promote_trait<T,complex<U> >::T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
         msres[i]=(-ms[i])+c;
    return MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote>(msres);
}   
   
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)
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
         msres[i]=(-ms[i])+s;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);   
}

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)
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
         msres[i]=(-ms[i])+u;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);   
}

// +++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator-=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator-=    
template<class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator-=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<int>& v)
{ ms.sub_assign(v); return ms;}

template<class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator-=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<double>& v)
{ ms.sub_assign(v); return ms;}

template<class U, class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator-=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<complex<U> >& v)
{ ms.sub_assign(v); return ms;}

// ++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator/
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- operator/     

template<class T, class ContainerT >   
MatrixResSeq<T> MatrixBaseSeq<T, ContainerT >::operator/(const T& t) const
{
    MatrixSeq<T> msres(*this);
    msres/=t;
    return MatrixResSeq<T>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,int>::T_promote> operator/\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<int>& v) 
{
    MatrixSeq<typename promote_trait<T,int>::T_promote> msres(ms);
    msres/=v;
    return MatrixResSeq<typename promote_trait<T,int>::T_promote>(msres);
}

template<class T, class ContainerT >   
MatrixResSeq<typename promote_trait<T,double>::T_promote> operator/\
    (const MatrixBaseSeq<T, ContainerT >& ms,const vector<double>& v) 
{
    MatrixSeq<typename promote_trait<T,double>::T_promote> msres(ms);
    msres/=v;
    return MatrixResSeq<typename promote_trait<T,double>::T_promote>(msres);
}

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) 
{
    MatrixSeq<typename promote_trait<T,complex<U> >::T_promote> msres(ms);
    msres/=v;
    return MatrixResSeq<typename promote_trait<T,complex<U> >::T_promote>(msres);
}

template<class T, class ContainerT > template<class U>  
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator/(const Scalar<U>& s) const 
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]/s;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U>  
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator/(const Matrix<U>& m) const 
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]/m;
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator/(const vector<Matrix<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]/ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator/(const vector<Scalar<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]/ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

template<class T, class ContainerT > template<class U> 
MatrixResSeq<typename promote_trait<T,U>::T_promote> \
    MatrixBaseSeq<T, ContainerT >::operator/(const vector<MatrixRes<U> >& ms) const
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(this->size());
    for(register int i=0;i<this->size();i++)
        msres[i]=(*this)[i]/ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

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)
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=s/ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}
   
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)
{
    MatrixSeq<typename promote_trait<T,U>::T_promote> msres(ms.size());
    for(register int i=0;i<ms.size();i++)
        msres[i]=u/ms[i];
    
    return MatrixResSeq<typename promote_trait<T,U>::T_promote>(msres);
}

// +++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::operator/=
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator/=    
template<class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator/=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<int>& v)
{ ms.div_assign(v); return ms;}

template<class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator/=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<double>& v)
{ ms.div_assign(v); return ms;}

template<class U, class T, class ContainerT>
MatrixBaseSeq<T,ContainerT>& operator/=\
    (MatrixBaseSeq<T,ContainerT>& ms,const vector<complex<U> >& v)
{ ms.div_assign(v); return ms;}

// +++++++++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::sum
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------------- sum
     
template<class T, class ContainerT > 
MatrixRes<T> MatrixBaseSeq<T, ContainerT >::sum() const
{
    Matrix<T> mres;
    int i;
    if (this->size()>0)
         mres=(*this)[0];
    for(i=1;i<this->size();i++)
    mres+=(*this)[i];
    return MatrixRes<T>(mres);
}
   
// ++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::plus_scan
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- plus_scan     
template<class T, class ContainerT > 
MatrixResSeq<T> MatrixBaseSeq<T, ContainerT >::plus_scan() const
{
    Matrix<T> msum;
    MatrixSeq<T> mseq;
    int i;
    if (this->size()>0)
    {
        msum=(*this)[0];
        mseq.push_back(msum);
        mseq[0].set_copyoption(SHALLOW);
    }
    for (i=1;i<this->size();i++)
    {
        msum+=(*this)[i];
        mseq.push_back(msum);
        mseq[i].set_copyoption(SHALLOW);
           
    }
    return MatrixResSeq<T>(mseq);
}

// ++++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::display
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- display     
template<class T, class ContainerT > 
void MatrixBaseSeq<T, ContainerT >::display() const 
{
    for_each(this->begin(),this->end(),mem_fun_ref(&Matrix<T>::display));
}
// +++++++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::norm
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ norm     
template<class T, class ContainerT > 
MatrixResSeq<T> norm(const MatrixBaseSeq<T, ContainerT >& ms)
{
   MatrixSeq<T> mseq;
   for(register int i=0;i<ms.size();i++)
   {
       const ContainerT& m=ms[i];
       const T& nrm=norm(m);
       mseq.push_back(nrm);
   }
   return MatrixResSeq<T>(mseq);
}
// +++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::seq_to_matrix
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------- seq_to_matrix
template<class T, class ContainerT >
Matrix<T> MatrixBaseSeq<T, ContainerT >::seq_to_matrix(size_t dimno) const
{
    size_t sz=this->size();
    if(sz==0)
        throw MatrixDomainError();
    size_t rk=this->elt(0).rank()+1;
    if(dimno<=0 || dimno>rk)
        throw MatrixIndexError();
    Dim d;
    DimStorageOption order=ROW_MAJOR_ORDER;
    if(rk>=1)
    {
      const ContainerT& elt0=this->elt(0);
      d=elt0.dimdata();
      order=elt0.storageorder();
    }
    valarray<size_t> dm_elt0(d.dims());
    valarray<size_t> dm(rk);
    for(register int i=0;i<dimno-1;i++)
        dm[i]=dm_elt0[i];
    dm[dimno-1]=sz;
    for(register int i=dimno;i<rk;i++)
        dm[i]=dm_elt0[i-1];
    Matrix<T> mat(Dim(dm,order));
    for(register int i=0;i<sz;i++)
        mat.indexdim(dimno,i)=this->elt(i);
    return MatrixRes<T>(mat);
} 
// ++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::is_homogeneous_seq
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------- is_homogeneous_seq
template<class T, class ContainerT >
bool MatrixBaseSeq<T, ContainerT >::is_homogeneous_seq() const 
{ 
// check that each elt of matrix seq has same rank and shape as first elt i.e.
// that matrix elts are homogeneous
    bool res=true;
    if (this->size()<2) 
        return res;
    const ContainerT& m=this->elt(0);
    int i=1;
    int sz=this->size();
    while(i<sz && res)
    {
        res=res && this->elt(i).samesize(m);
        i++;
    }
    return res;
}
// +++++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::is_scalar_seq
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------- is_scalar_seq
template<class T, class ContainerT >
bool MatrixBaseSeq<T, ContainerT >::is_scalar_seq() const 
{ 
// returns true if each element of seq is a scalar
    int sz=this->size();
    if(!sz)
      return false;
    bool res=true;
    int i=0;
    while(i<sz && res)
    {
        res=res && !this->elt(i).rank();
        i++;
    }
    return res;
}
// +++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::get_element_seq
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ get_element_seq 
template<class T, class ContainerT >
MatrixResSeq<T> \
MatrixBaseSeq<T, ContainerT >::get_element_seq(const vector<int>& v) const
{
    if(!is_homogeneous_seq())
        throw MatrixDomainError(\
            "MatrixBaseSeq<T>::get_element_seq]: sequence not homogeneous");
    const Matrix<T>& elt0=this->elt(0);
    if (v.size()!=elt0.rank())
        throw MatrixRankError();       
    valarray<size_t> sz(elt0.dims());
    valarray<size_t> st(computestrides(sz));
    for(register int j=0;j<sz.size();j++)
      sz[j]=0;
    size_t off=inner_product(v.begin(),v.end(),&st[0],0);
    gslice g(off,sz,st);
    return get_submatrix_seq(g);
}

template<class T, class ContainerT >
MatrixResSeq<T> \
MatrixBaseSeq<T, ContainerT >::get_element_seq(valarray<size_t> v) const
{
    if(!is_homogeneous_seq())
        throw MatrixDomainError(\
            "MatrixBaseSeq<T>::get_element_seq]: sequence not homogeneous");
    int num=v.size();
    const Matrix<T>& elt0=this->elt(0);
    if (num!=elt0.rank())
        throw MatrixRankError();       
    valarray<size_t> sz(elt0.dims());
    valarray<size_t> st(computestrides(sz));
    for(register int j=0;j<sz.size();j++)
      sz[j]=0;
    size_t off=inner_product(&v[0],&v[num],&st[0],0);
    gslice g(off,sz,st);
    return get_submatrix_seq(g);
  
}
// +++++++++++++++++++++++++++++++++++++++ class MatrixBaseSeq::get_submatrix_seq
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ get_submatrix_seq 
template<class T, class ContainerT >
MatrixResSeq<T> \
MatrixBaseSeq<T, ContainerT >::get_submatrix_seq(const gslice& g) const
{
    if(!is_homogeneous_seq())
        throw MatrixDomainError(\
            "MatrixBaseSeq<T>::get_submatrix_seq]: sequence not homogeneous");
    int num=this->size();
    MatrixSeq<T> mseq(num);
    valarray<size_t> sz(g.size());
    valarray<size_t> st(g.stride());
    valarray<size_t> fact(g.size().size());
    int off=g.start();

    for(register int i=0;i<num;i++)
    {
      if(i==0)
        fact=computestrides(this->elt(0).dims())/st;
      gslice elt_g=this->elt(i).getgslice();
      valarray<size_t> elt_st=elt_g.stride();
      for(register int j=0;j<elt_st.size();j++)
        elt_st[j]*=fact[j];
      int elt_off=off+elt_g.start();
      elt_g=gslice(elt_off,sz,elt_st);
      mseq[i]=this->elt(i)(g).duplicate(DEEP);
    }
    return MatrixResSeq<T>(mseq);
} 

// ++++++++++++++++++++++++++++++++++++++++++++++ class MatrixSubSeq::operator[]
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------------ operator[]

template<class T>
MatrixSubSeq<T> \
MatrixSeq<T>::operator[](const DummyArg& arg) const
{ return MatrixSubSeq<T>(*this);}

template<class T>
MatrixSubSeq<T> \
MatrixSubSeq<T>::operator[](size_t idx) const
{
    if(!this->is_homogeneous_seq())
        throw MatrixDomainError(\
            "MatrixSubSeq<T>::operator[]: sequence not homogeneous");
    int sz=this->size();
    MatrixSeq<T> mseq(sz);
    for(register int i=0;i<sz;i++)
      mseq[i]=this->elt(i)[idx].duplicate(SHALLOW);
    return MatrixSubSeq<T>(mseq);
}

template<class T>
MatrixSubSeq<T> \
MatrixSubSeq<T>::operator[](const DummyArg& arg) const
{
    if(!this->is_homogeneous_seq())
        throw MatrixDomainError(\
            "MatrixSubSeq<T>::operator[]: sequence not homogeneous");
    int sz=this->size();
    MatrixSeq<T> mseq(sz);
    for(register int i=0;i<sz;i++)
        mseq[i]=this->elt(i)[arg].duplicate(SHALLOW);
    return MatrixSubSeq<T>(mseq);
}

template<class T>
MatrixSubSeq<T> 
MatrixSubSeq<T>::operator[](const slice& sl) const
{
    if(!this->is_homogeneous_seq())
        throw MatrixDomainError(\
            "MatrixSubSeq<T>::operator[]: sequence not homogeneous");
    int sz=this->size();
    MatrixSeq<T> mseq(sz);
    for(register int i=0;i<sz;i++)
        mseq[i]=this->elt(i)[sl].duplicate(SHALLOW);
    return MatrixSubSeq<T>(mseq);
}

// +++++++++++++++++++++++++++++++++++++++++++++ class ostream::operator<<
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator<<
template<class T>
inline ostream& operator<< (ostream& os,const Scalar<T>& s)
{
  T val=s.scalar();
  int wid=os.width();
  if (wid==0) wid=10;
  //if(!(os.flags() & ios::fixed) && !(os.flags() & ios::scientific))
  //    os.setf(ios::fixed,ios::floatfield);
  os << setw(wid) << val ;
  return os;
}

template<class T>
inline ostream& operator<< (ostream& os,const Matrix<T>& m)
{
    typename Matrix<T>::iterator itr=m.begin();
    valarray<size_t> str(m.rank());
    //size_t ind=itr.index();
    size_t pos=itr.eltno();  // will be zero
    register int i;
    int wid=os.width();
    if (wid==0) wid=10;
    //if(!(os.flags() & ios::fixed) && !(os.flags() & ios::scientific))
    //    os.setf(ios::fixed,ios::floatfield);

                  
    while(itr<m.end())
    {
        if (m.rank()>0)
        {
            //valarray<size_t> str(m.rank());
            if(itr==m.begin())
            {
            // output dimension info
                valarray<size_t> mdm(m.dims());
                os << endl << "[ " << mdm[0];
                for(i=1;i<m.rank();i++)
                    os << " x " << mdm[i] ;
                os << " ]" << endl;
                str=computestrides(mdm);               
            } 
            // for each dimension (except lowest) see if newline character reqd
            for (i=m.rank()-1;i>0;i--)
            {  
                if ( pos!=0 && pos % str[i-1]==0)
                    os << endl;     
            }
            os << " ";
        }      
        T val=*itr;//(*m.data_ptr())[ind];
        os << setw(wid) << val ;//<< " ";
        itr++;
        //ind=itr.index();
        pos=itr.eltno();
               
    } 
    //os << endl;
    return os;
}

template<class T,class ContainerT >
ostream& operator<<(ostream& os, const MatrixBaseSeq<T, ContainerT >& ms)
{
   for(register int i=0;i<ms.size();i++)       
       os << ms[i] << endl;
   return os;
}

// +++++++++++++++++++++++++++++++++++++++++++++ class istream::operator>>
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ operator>> 

template<class T>
inline istream& operator>> (istream& is,const Matrix<T>& m)
{
      typename Matrix<T>::iterator itr;    
      for(itr=m.begin();itr!=m.end();itr++)
      {
          T& elt=*itr;
          is >> elt;
      }
     
      return is;     
}

template<class T,class ContainerT >
istream& operator>> (istream& is, const MatrixBaseSeq<T, ContainerT >& ms)
{
   for(register int i=0;i<ms.size();i++)       
       is >> ms[i];
   return is;
}


}  // end namespace QDOS



#endif  // QDOS_MATRIX_CC

