Matrix Code

迈不过友情╰ 2022-09-25 05:26 222阅读 0赞

## Matrix.h ##

/************************************************************************
    * Copyright(c) 2012  Gavin Liu
    * All rights reserved.
    *
    * File: matrix.h
    * Brief: 
    * Version: 1.0
    * Author: Gavin Liu
    * Email: gavinliuyi@msn.com
    * Date: 2012/09/28
    * History:
    ************************************************************************/
    #ifndef _MATRIX_H
    #define _MATRIX_H
    
    #if defined(DEBUG) || defined(_DEBUG)
    #define RANGE_CHECK_
    #endif
    
    // For MFC projects in VC ++ !!!
    #if defined(max) || defined(min)
    #undef max
    #undef min
    #endif
    
    #if defined(abs)
    #undef abs
    #endif
    
    #include <cmath>
    #include <valarray>
    #include <complex>
    #include <limits>
    #include <stdexcept>
    #include <algorithm>
    #include <iomanip>
    
    
    #ifndef MAX_MATRIX_CACHE
    #define MAX_MATRIX_CACHE    4
    #endif
    
    
    namespace techsoft {
    
    using std::istream;
    using std::ostream;
    using std::valarray;
    using std::complex;
    using std::numeric_limits;
    using std::slice;
    using std::slice_array;
    using std::abs;
    using std::sqrt;
    using std::swap;
    using std::setw;
    
    #ifdef RANGE_CHECK_
    using std::out_of_range;
    using std::invalid_argument;
    using std::runtime_error;
    using std::exception;
    #endif
    
    #if defined(_MSC_VER)
    using ::ceil;
    using ::floor;
    using ::srand;
    using ::rand;
    using ::sqrt;
    using ::sqrtl;
    #else
    using std::ceil;
    using std::floor;
    using std::srand;
    using std::rand;
    using std::size_t;
    #endif
    
    //
    // Temporary helping classes
    //
    class mslice;
    class gmslice;
    template <class T> class Mref;
    template <class T> class Cmat_iter;
    template <class T> class Mat_iter;
    template <class T> class mslice_matrix;
    template <class T> class gmslice_matrix;
    
    enum MatArrayType { C_ARRAY, FORTRAN_ARRAY };
    
    template <class T>
    class matrix
    {
    public:
       typedef T value_type;
       typedef valarray<T> array_type;
    
       // Constructors
       matrix ();
       matrix (size_t nRow, size_t nCol);
       matrix (size_t nRow, size_t nCol, const T& e);
       matrix (size_t nRow, size_t nCol, const T* Array, MatArrayType aType = C_ARRAY);
       matrix (size_t nRow, size_t nCol, const valarray<T>& vec);
       ~matrix ();
    
       // Copy constructors
       matrix (const matrix<T>& m);
       matrix (const mslice_matrix<T>& sm);
       matrix (const gmslice_matrix<T>& sm);
    #ifndef _TS_NO_MEMBER_TEMPLATES
       template <class X> matrix (const matrix<X>& m);
    #endif
    
       // Assignment operators
       matrix<T>& operator= (const matrix<T>& m);
       matrix<T>& operator= (const mslice_matrix<T>& sm);
       matrix<T>& operator= (const gmslice_matrix<T>& sm);
       matrix<T>& operator= (const T& e);
    #ifndef _TS_NO_MEMBER_TEMPLATES
       template <class X> matrix<T>& operator= (const matrix<X>& m);
    #endif
    
       size_t size () const throw() { return mPtr->nrow * mPtr->ncol; }
       size_t typesize () const throw() { return sizeof(T); }
    
       size_t colsize () const throw() { return mPtr->nrow; }   // No. of elements in a column, i.e., it's row number
       size_t rowsize () const throw() { return mPtr->ncol; }   // No. of elements in a row, i.e., it's column number
    
       size_t rowno () const throw() { return mPtr->nrow; }     // No. of rows
       size_t colno () const throw() { return mPtr->ncol; }     // No. of columns
    
       void resize (size_t nRow, size_t nCol, const T& dval = T(0));
       void free ();
    
       Mat_iter<T> row (size_t i);
       Cmat_iter<T> row (size_t i) const;
       Mat_iter<T> column (size_t i);
       Cmat_iter<T> column (size_t i) const;
       Mat_iter<T> diag (int i=0);
       Cmat_iter<T> diag (int i=0) const;
       T trace (int i=0) const;
    
       Mat_iter<T> operator[] (size_t i) { return row(i); }
       Cmat_iter<T> operator[] (size_t i) const { return row(i); }
       Mat_iter<T> operator[] (int i) { return row(i); }
       Cmat_iter<T> operator[] (int i) const { return row(i); }
    
       Mat_iter<T> operator() (size_t i) { return column(i); }
       Cmat_iter<T> operator() (size_t i) const { return column(i); }
       Mat_iter<T> operator() (int i) { return column(i); }
       Cmat_iter<T> operator() (int i) const { return column(i); }
    
       Mref<T> operator() (size_t i, size_t j);
       const T& operator() (size_t i, size_t j) const;
    
       mslice_matrix<T> operator[] (mslice ms);
       gmslice_matrix<T> operator[] (gmslice gm);
    
       const matrix<T> operator[] (mslice ms) const;
       const matrix<T> operator[] (gmslice gm) const;
    
       // Unary operators
       matrix<T> operator+ () const { return *this; }
       matrix<T> operator- () const;
       matrix<T> operator~ () const;
       matrix<T> operator! () const;
    
       // Computed assignment
       matrix<T>& operator+= (const matrix<T>& m);
       matrix<T>& operator-= (const matrix<T>& m);
       matrix<T>& operator*= (const matrix<T>& m);
       matrix<T>& operator/= (const matrix<T>& m);
       matrix<T>& operator*= (const T& x);
       matrix<T>& operator/= (const T& x);
    
       // Miscellaneous methods
       void null ();
       void unit ();
       void rand (int rmin=-1, int rmax=1, int rseed=0);
    
       T sum () const { return mPtr->val.sum(); }
       T min () const { return mPtr->val.min(); }
       T max () const { return mPtr->val.max(); }
    
       matrix<T> apply (T (*fn)(T)) const;
       matrix<T> apply (T (*fn)(const T&)) const;
    
       matrix<T> apply (T (*fn)(size_t,size_t,T)) const;
       matrix<T> apply (T (*fn)(size_t,size_t,const T&)) const;
    
       // Utility methods
       bool solve (const valarray<T>& v, valarray<T>& s) const;
       bool solve (const matrix<T>& v, matrix<T>& s) const;
    
       matrix<T> adj () const;
       T cofact (size_t row, size_t col) const;
    
       T det () const;
       T cond () const;
       size_t rank () const;
    
       T norm1 () const;
       T norm2 () const;
       T normI () const;
       T normF () const;
    
       bool lud (valarray<size_t>& ri, T* pDet = NULL);
       void lubksb (const valarray<size_t>& ri, const valarray<T>& v, valarray<T>& s) const;
       void lumpove (const matrix<T>& ludm, const valarray<size_t>& ri, const valarray<T>& v, valarray<T>& s) const;
       bool solve_lu (const valarray<T>& v, valarray<T>& s) const { return solve( v, s); }
    
       bool inv ();
       bool inv_lu ();
       bool inv_sv ();
       bool inv_qr ();
    
       bool svd (matrix<T>& vc, valarray<T>& w);
       void svbksb (const matrix<T>& vc, const valarray<T>& w, const valarray<T>& b, valarray<T>& x) const;
       bool solve_sv (const valarray<T>& v, valarray<T>& s) const;
    
       void qrd (matrix<T>& r);
       void qrbksb (const matrix<T>& r, const valarray<T>& v, valarray<T>& s) const;
       bool solve_qr (const valarray<T>& v, valarray<T>& s) const;
    
       bool chold ();
       void cholbksb (const valarray<T>& v, valarray<T>& s) const;
       bool solve_chol (const valarray<T>& v, valarray<T>& s) const;
    
       bool eigen (valarray<T>& eival) const;
       bool eigen (valarray<T>& eival, matrix<T>& eivec) const;
    
       bool eigen (valarray<T>& rev, valarray<T>& iev) const;
       bool eigen (valarray<T>& rev, valarray<T>& iev, matrix<T>& eivec) const;
    
       // Type of matrices
       bool isSquare () const throw() { return (mPtr->nrow == mPtr->ncol); }
       bool isSingular () const;
       bool isDiagonal () const;
       bool isScalar () const;
       bool isUnit () const throw();
       bool isNull () const throw();
       bool isSymmetric () const;
       bool isSkewSymmetric () const;
       bool isUpperTriangular () const;
       bool isLowerTriangular () const;
       bool isRowOrthogonal () const;
       bool isColOrthogonal () const;
    
    private:
    
       struct base_mat
       {
          valarray<T> val;
          size_t nrow, ncol;
          int    refcnt;
    
          base_mat (size_t row, size_t col)
                : nrow( row), ncol( col), val( row * col) { refcnt = 1; }
    
          base_mat (const T& v, size_t row, size_t col)
                : nrow( row), ncol( col), val( v, row * col) { refcnt = 1; }
    
          base_mat (const T* v, size_t row, size_t col)
                : nrow( row), ncol( col), val( v, row * col) { refcnt = 1; }
    
          base_mat (size_t row, size_t col, const valarray<T>& v)
                : nrow( row), ncol( col), val( v) { refcnt = 1; }
    
          ~base_mat () {;}
       };
    
       base_mat *mPtr;
       enum AllocType { Allocate, Free, Record };
    
       bool shared () { return (mPtr->refcnt > 1); }
       void clone ();
       base_mat *allocator (AllocType bOpt, size_t nrow = 0, size_t ncol = 0);
    
       void tred2 (valarray<T>& d, valarray<T>& e, bool eivec);
       bool tql2 (valarray<T>& d, valarray<T>& e, bool eivec);
       void balanc (matrix<T>& v, bool eivec);
       bool hqr2 (valarray<T>& d, valarray<T>& e, matrix<T>& v, bool eivec);
    
       friend class Mat_iter<T>;
       friend class Cmat_iter<T>;
       friend class Mref<T>;
    };
    
    
    /
    // Non-member binary operators
    //
    template <class T> matrix<T> operator+ (const matrix<T>& m1, const matrix<T>& m2);
    template <class T> matrix<T> operator- (const matrix<T>& m1, const matrix<T>& m2);
    
    template <class T> matrix<T> operator* (const matrix<T>& m1, const matrix<T>& m2);
    template <class T> matrix<T> operator* (const matrix<T>& m, const T& el);
    template <class T> matrix<T> operator* (const T& el, const matrix<T>& m);
    
    template <class T> matrix<T> operator/ (const matrix<T>& m1, const matrix<T>& m2);
    template <class T> matrix<T> operator/ (const matrix<T>& m, const T& el);
    template <class T> matrix<T> operator/ (const T& el, const matrix<T>& m);
    
    template <class T> valarray<T> operator* (const matrix<T>& m, const valarray<T>& v);
    template <class T> valarray<T> operator* (const valarray<T>& v, const matrix<T>& m);
    template <class T> valarray<T> operator/ (const matrix<T>& m, const valarray<T>& v);
    template <class T> valarray<T> operator/ (const valarray<T>& v, const matrix<T>& m);
    
    template <class T> bool operator== (const matrix<T>& m1, const matrix<T>& m2);
    template <class T> bool operator!= (const matrix<T>& m1, const matrix<T>& m2);
    
    
    /
    // Non-member functions
    //
    template <class T> matrix<T> pow (const matrix<T>& m, size_t n);
    template <class T> matrix<T> abs (const matrix<T>& m);
    template <class T> matrix<T> floor (const matrix<T>& m);
    template <class T> matrix<T> ceil (const matrix<T>& m);
    
    template <class T> void mswap (const Mref<T>& x, const Mref<T>& y);
    template <class T> void mswap (const Mat_iter<T>& x, const Mat_iter<T>& y);
    template <class T> void mswap (matrix<T>& x, matrix<T>& y);
    
    template <class T> istream& operator>> (istream& is, Mref<T> el);
    template <class T> istream& operator>> (istream& is, valarray<T>& v);
    template <class T> istream& operator>> (istream& is, Mat_iter<T> v);
    template <class T> istream& operator>> (istream& is, matrix<T>& m);
    
    template <class T> ostream& operator<< (ostream& os, const Mref<T>& el);
    template <class T> ostream& operator<< (ostream &os, const valarray<T>& v);
    template <class T> ostream& operator<< (ostream& os, const Mat_iter<T>& v);
    template <class T> ostream& operator<< (ostream& os, const Cmat_iter<T>& v);
    template <class T> ostream& operator<< (ostream &os, const matrix<T>& m);
    
    float epsilon (complex<float> v);
    double epsilon (complex<double> v);
    long double epsilon (complex<long double> v);
    
    template <class T> T epsilon (const T& v);
    template <class T> bool isVecEq (const valarray<T>& v1, const valarray<T>& v2);
    
    
    //
    // Temporary helping classes
    //
    
    template <class T>
    class Mref
    {
    public:
       Mref (matrix<T>& mm, size_t ii) : m(mm), i(ii) {;}
    
       operator T () const { return m.mPtr->val[i]; }
       T* operator& () const { if (m.shared()) m.clone(); return &m.mPtr->val[i]; }
    
       T& operator= (const Mref<T>& e) const { if (m.shared()) m.clone(); T te = e; return m.mPtr->val[i] = te; }
       T& operator= (const T& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] = e; }
       T& operator+= (const T& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] += e; }
       T& operator-= (const T& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] -= e; }
       T& operator*= (const T& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] *= e; }
       T& operator/= (const T& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] /= e; }
    
    #ifndef _TS_NO_MEMBER_TEMPLATES
       template <class X> T& operator= (const X& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] = e; }
       template <class X> T& operator+= (const X& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] += e; }
       template <class X> T& operator-= (const X& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] -= e; }
       template <class X> T& operator*= (const X& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] *= e; }
       template <class X> T& operator/= (const X& e) const { if (m.shared()) m.clone(); return m.mPtr->val[i] /= e; }
    
       template <class X> T operator+ (const X& e) { T y = *this; return y += e;  }
       template <class X> T operator- (const X& e) { T y = *this; return y -= e;  }
       template <class X> T operator* (const X& e) { T y = *this; return y *= e;  }
       template <class X> T operator/ (const X& e) { T y = *this; return y /= e;  }
    #endif
    
       T operator++ () { if (m.shared()) m.clone(); return ++m.mPtr->val[i]; }
       T operator++ (int) { if (m.shared()) m.clone(); return m.mPtr->val[i]++; }
       T operator-- () { if (m.shared()) m.clone(); return --m.mPtr->val[i]; }
       T operator-- (int) { if (m.shared()) m.clone(); return m.mPtr->val[i]--; }
    
       T operator+ () { return m.mPtr->val[i]; }
       T operator- () { return -m.mPtr->val[i]; }
    
       Mref (const Mref& mr) : m(mr.m), i(mr.i) {;}
    private:
       Mref ();
       matrix<T>& m;
       size_t i;
    };
    
    
    template <class T>
    class Cmat_iter
    {
    public:
       Cmat_iter (const matrix<T>& mm, const slice& ss) : m(mm), s(ss) {;}
    
       const T& operator[] (size_t i) const;
       operator valarray<T> () const { return m.mPtr->val[s]; }
       size_t size () const { return s.size(); }
    
       Cmat_iter(const Cmat_iter& mi) : m(mi.m), s(mi.s) {;}
    private:
       Cmat_iter();
       Cmat_iter& operator= (const Cmat_iter&);
    
       const matrix<T>& m;
       slice s;
    };
    
    
    template <class T>
    class Mat_iter
    {
    public:
       Mat_iter (matrix<T>& mm, const slice& ss) : m(mm), s(ss) {;}
    
       Mref<T> operator[] (size_t i) const;
       operator valarray<T>() const { return m.mPtr->val[s]; }
    
       void operator= (const valarray<T>& v) const;
       void operator= (const T& e) const;
       void operator*= (const valarray<T>& v) const;
       void operator/= (const valarray<T>& v) const;
       void operator+= (const valarray<T>& v) const;
       void operator-= (const valarray<T>& v) const;
       size_t size () const { return s.size(); }
    
       Mat_iter(const Mat_iter& mi) : m(mi.m), s(mi.s) {;}
    private:
       Mat_iter();
       Mat_iter<T>& operator= (const Mat_iter<T>& mat);
    
       matrix<T>& m;
       slice s;
    };
    
    class mslice
    {
    public:
       mslice (size_t sRow, size_t sCol, size_t nRow, size_t nCol)
          : srow_(sRow), scol_(sCol), nrow_(nRow), ncol_(nCol) {;}
       mslice (const mslice& ms)
          : srow_(ms.srow_), scol_(ms.scol_),
            nrow_(ms.nrow_), ncol_(ms.ncol_) {;}
    
       size_t size () const { return nrow_ * ncol_; }
       size_t start (size_t colno) const { return srow_*colno+scol_; }
       size_t end (size_t colno) const { return (srow_+nrow_-1)*colno+(scol_+ncol_); }
    
       size_t colsize () const { return nrow_; }
       size_t rowsize () const { return ncol_; }
    
       size_t rowno () const { return nrow_; }
       size_t colno () const { return ncol_; }
    
    private:
       mslice ();
       size_t srow_,scol_, nrow_, ncol_;
    };
    
    template <class T>
    class mslice_matrix
    {
    public:
       typedef T value_type;
       mslice_matrix (matrix<T>& mm, const mslice& ss)
          : s(ss), m(mm) {;}
    
       const mslice& get_slice() const { return s; }
       matrix<T>& get_ref_mat () const { return m; }
    
       void operator= (const matrix<T>& m2) const;
       void operator= (const T& e) const;
       void operator+= (const matrix<T>& m2) const;
       void operator-= (const matrix<T>& m2) const;
       void operator*= (const T& e) const;
       void operator/= (const T& e) const;
    
       mslice_matrix (const mslice_matrix& ms) : s(ms.s), m(ms.m) {;}
    private:
       mslice_matrix ();
       mslice_matrix& operator= (const mslice_matrix&);
    
       mslice s;
       matrix<T>& m;
    };
    
    
    enum MATIRX_TYPE { GENERAL, DIAGONAL, TRIDIAGONAL, UTRIANG, LTRIANG };
    
    class gmslice
    {
    public:
       gmslice (MATIRX_TYPE mtype) 
          : MType_(mtype), ncol_(0) {;}
       gmslice (size_t nCol, const valarray<size_t>& sRow, const valarray<size_t>& nRow)
          : MType_(GENERAL), ncol_(nCol), srow_(sRow), rsize_(nRow) {;}
       gmslice (const gmslice& gm)
          : MType_(gm.MType_), ncol_(gm.ncol_),
            srow_(gm.srow_), rsize_(gm.rsize_) {;}
    
       size_t size () const { return srow_.size() * ncol_; }
       size_t start (size_t rowno) const { return srow_[rowno]; }
       size_t end (size_t rowno) const { return srow_[rowno]+rsize_[rowno]; }
    
    
       size_t colsize () const { return srow_.size(); }
       size_t rowsize () const { return ncol_; }
    
       size_t rowno () const { return srow_.size(); }
       size_t colno () const { return ncol_; }
    
       MATIRX_TYPE getype() { return MType_; }
    
       size_t ncol_;
       valarray<size_t> srow_,rsize_;
    
    private:
       gmslice();
       MATIRX_TYPE MType_;
    };
    
    
    
    template <class T>
    class gmslice_matrix
    {
    public:
       typedef T value_type;
       gmslice_matrix (matrix<T>& mm, const gmslice& ss)
          : s(ss), m(mm) {;}
    
       const gmslice& get_slice() const { return s; }
       matrix<T>& get_ref_mat () const { return m; }
    
       void operator= (const matrix<T>& m2) const;
       void operator= (const T& e) const;
       void operator+= (const matrix<T>& m2) const;
       void operator-= (const matrix<T>& m2) const;
       void operator*= (const T& e) const;
       void operator/= (const T& e) const;
    
       gmslice_matrix (const gmslice_matrix& ms) : s(ms.s), m(ms.m) {;}
    private:
       gmslice_matrix ();
       gmslice_matrix& operator= (const gmslice_matrix&);
    
       gmslice s;
       matrix<T>& m;
    };
    
    
    }   // namespace techsoft
    
    //#include "matrix.cpp"
    
    #endif // _MATRIX_H

--------------------

## Matrix.cpp ##

/************************************************************************
    * Copyright(c) 2012  Gavin Liu
    * All rights reserved.
    *
    * File: matrix.cpp
    * Brief: 
    * Version: 1.0
    * Author: Gavin Liu
    * Email: gavinliuyi@msn.com
    * Date: 2012/09/28
    * History:
    ************************************************************************/
    #ifndef _MATRIX_CPP
    #define _MATRIX_CPP
    
    #include"matrix.h"
    
    namespace techsoft {
    
    
    template<class T> inline
    const T& max (const T& x, const T& y)
    {
       return (x < y) ? y : x;
    }
    
    template<class T> inline
    const T& min (const T& x, const T& y)
    {
       return (y < x) ? y : x;
    }
    
    template <class T> inline 
    void mswap (const Mref<T>& x, const Mref<T>& y)
    {
       T z = x;
       x = y;
       y = z;
    }
    
    template <class T> inline 
    void mswap (const Mat_iter<T>& x, const Mat_iter<T>& y)
    {
       size_t n = min( x.size(), y.size());
       for (size_t i=0; i < n; i++)
       {
          T z = x[i];
          x[i] = y[i];
          y[i] = z;
       }
    }
    
    template <class T> inline 
    void mswap (matrix<T>& x, matrix<T>& y)
    {
       matrix<T> z = x;
       x = y;
       y = z;
    }
    
    inline float epsilon (complex<float>) 
    {
       return std::numeric_limits<float>::epsilon() * 100.0f;
    }
    
    inline double epsilon (complex<double>) 
    {
       return std::numeric_limits<double>::epsilon() * 1e4;
    }
    
    inline long double epsilon (complex<long double>) 
    {
       return std::numeric_limits<long double>::epsilon() * 1e4;
    }
    
    inline double epsilon (const double& v) 
    { 
       double ep = std::numeric_limits<double>::epsilon() * 1e4;
       return v > 1.0 ? v * ep : ep;
    }
    
    inline long double epsilon (const long double& v) 
    { 
       long double ep = std::numeric_limits<long double>::epsilon() * 1e4;
       return v > 1.0 ? v * ep : ep;
    }
    
    template <class T> inline 
    T epsilon (const T& v) 
    { 
       T e = std::numeric_limits<T>::epsilon() * 100;
       return v > T(1) ? v * e : e;
    }
    
    //
    // Implementation of Cmat_iter class
    //
    
    template <class T> inline const T&
    Cmat_iter<T>::operator[] (size_t i) const
    {
       #ifdef RANGE_CHECK_
       if (i >= s.size())
          throw out_of_range( "matrix<T>::operator[][]: Subscript out of range!");
       #endif
       return m.mPtr->val[s.start()+i*s.stride()];
    }
    
    
    //
    // Implementation of Mat_iter class
    //
    
    template <class T> inline Mref<T>
    Mat_iter<T>::operator[] (size_t i) const
    {
       #ifdef RANGE_CHECK_
       if (i >= s.size())
          throw out_of_range( "matrix<T>::operator[][]: Subscript out of range!");
       #endif
       return Mref<T>( m, s.start()+i*s.stride());
    }
    
    template <class T> inline void
    Mat_iter<T>::operator= (const valarray<T>& v) const
    {
       #ifdef RANGE_CHECK_
       if (size() != v.size())
          throw invalid_argument( "matrix<T>::operator[]: size mismatch in assignment operation!");
       #endif
    
       if (m.shared())
          m.clone();
       m.mPtr->val[s] = v;
    }
    
    template <class T> inline void
    Mat_iter<T>::operator= (const T& e) const
    {
       if (m.shared())
          m.clone();
       m.mPtr->val[s] = e;
    }
    
    template <class T> inline void
    Mat_iter<T>::operator*= (const valarray<T>& v) const
    {
       #ifdef RANGE_CHECK_
       if (size() != v.size())
          throw invalid_argument( "matrix<T>::operator*=: size mismatch!");
       #endif
    
       if (m.shared())
          m.clone();
       m.mPtr->val[s] *= v;
    }
    
    template <class T> inline void
    Mat_iter<T>::operator/= (const valarray<T>& v) const
    {
       #ifdef RANGE_CHECK_
       if (size() != v.size())
          throw invalid_argument( "matrix<T>::operator/=: size mismatch!");
       #endif
    
       if (m.shared())
          m.clone();
       m.mPtr->val[s] /= v;
    }
    
    template <class T> inline void
    Mat_iter<T>::operator+= (const valarray<T>& v) const
    {
       #ifdef RANGE_CHECK_
       if (size() != v.size())
          throw invalid_argument( "matrix<T>::operator+=: size mismatch!");
       #endif
    
       if (m.shared())
          m.clone();
       m.mPtr->val[s] += v;
    }
    
    template <class T> inline void
    Mat_iter<T>::operator-= (const valarray<T>& v) const
    {
       #ifdef RANGE_CHECK_
       if (size() != v.size())
          throw invalid_argument( "matrix<T>::operator-=: size mismatch!");
       #endif
    
       if (m.shared())
          m.clone();
       m.mPtr->val[s] -= v;
    }
    
    
    //
    // Implementation of mslice_matrix class
    //
    
    template <class T>
    void mslice_matrix<T>::operator= (const matrix<T>& m2) const
    {
       #ifdef RANGE_CHECK_
       if (s.rowno() != m2.rowno() || s.colno() != m2.colno())
          throw out_of_range( "mslice_matrix<T>operator= Non-conforming matrix size!");
       #endif
    
       T *pv = &m(0,0);
       size_t start = s.start( m.colno());
    
       for (size_t i=0; i < s.rowno(); i++)
       {
          size_t ipos = m.colno() * i;
          for (size_t j=0; j < s.colno(); j++)
             pv[start+ipos+j] = m2(i,j);
       }
    }
    
    template <class T>
    void mslice_matrix<T>::operator+= (const matrix<T>& m2) const
    {
       #ifdef RANGE_CHECK_
       if (s.rowno() != m2.rowno() || s.colno() != m2.colno())
          throw out_of_range( "mslice_matrix<T>::operator+= Non-conforming matrix size!");
       #endif
    
       T *pv = &m(0,0);
       size_t start = s.start( m.colno());
    
       for (size_t i=0; i < s.rowno(); i++)
       {
          size_t ipos = m.colno() * i;
          for (size_t j=0; j < s.colno(); j++)
             pv[start+ipos+j] += m2(i,j);
       }
    }
    
    template <class T>
    void mslice_matrix<T>::operator-= (const matrix<T>& m2) const
    {
       #ifdef RANGE_CHECK_
       if (s.rowno() != m2.rowno() || s.colno() != m2.colno())
          throw out_of_range( "mslice_matrix<T>::operator+= Non-conforming matrix size!");
       #endif
    
       T *pv = &m(0,0);
       size_t start = s.start( m.colno());
    
       for (size_t i=0; i < s.rowno(); i++)
       {
          size_t ipos = m.colno() * i;
          for (size_t j=0; j < s.colno(); j++)
             pv[start+ipos+j] -= m2(i,j);
       }
    }
    
    template <class T>
    void mslice_matrix<T>::operator= (const T& e) const
    {
       T *pv = &m(0,0);
       size_t start = s.start( m.colno());
    
       for (size_t i=0; i < s.rowno(); i++)
       {
          size_t ipos = m.colno() * i;
          for (size_t j=0; j < s.colno(); j++)
             pv[start+ipos+j] = e;
       }
    }
    
    template <class T>
    void mslice_matrix<T>::operator*= (const T& e) const
    {
       T *pv = &m(0,0);
       size_t start = s.start( m.colno());
    
       for (size_t i=0; i < s.rowno(); i++)
       {
          size_t ipos = m.colno() * i;
          for (size_t j=0; j < s.colno(); j++)
             pv[start+ipos+j] *= e;
       }
    }
    
    template <class T>
    void mslice_matrix<T>::operator/= (const T& e) const
    {
       T *pv = &m(0,0);
       size_t start = s.start( m.colno());
    
       for (size_t i=0; i < s.rowno(); i++)
       {
          size_t ipos = m.colno() * i;
          for (size_t j=0; j < s.colno(); j++)
             pv[start+ipos+j] /= e;
       }
    }
    
    
    //
    // Implementation of mslice_matrix class
    //
    
    template <class T>
    void gmslice_matrix<T>::operator= (const matrix<T>& m2) const
    {
       #ifdef RANGE_CHECK_
       if (s.rowno() != m2.rowno() || s.colno() != m2.colno())
          throw out_of_range( "mslice_matrix<T>operator= Non-conforming matrix size!");
       #endif
    
       T *pv = &m(0,0);
       const T *pm = &m2(0,0);
    
       for (size_t i=0,ipos=0; i < s.rowno(); ++i,ipos+=s.colno())
          for (size_t j=s.start(i); j < s.end(i); j++)
          {
             size_t jpos = ipos+j;
             pv[jpos] = pm[jpos];
          }
    }
    
    template <class T>
    void gmslice_matrix<T>::operator+= (const matrix<T>& m2) const
    {
       #ifdef RANGE_CHECK_
       if (s.rowno() != m2.rowno() || s.colno() != m2.colno())
          throw out_of_range( "mslice_matrix<T>::operator+= Non-conforming matrix size!");
       #endif
    
       T *pv = &m(0,0);
       const T *pm = &m2(0,0);
    
       for (size_t i=0,ipos=0; i < s.rowno(); ++i,ipos+=s.colno())
          for (size_t j=s.start(i); j < s.end(i); j++)
          {
             size_t jpos = ipos+j;
             pv[jpos] += pm[jpos];
          }
    }
    
    template <class T>
    void gmslice_matrix<T>::operator-= (const matrix<T>& m2) const
    {
       #ifdef RANGE_CHECK_
       if (s.rowno() != m2.rowno() || s.colno() != m2.colno())
          throw out_of_range( "mslice_matrix<T>::operator+= Non-conforming matrix size!");
       #endif
    
       T *pv = &m(0,0);
       const T *pm = &m2(0,0);
    
       for (size_t i=0,ipos=0; i < s.rowno(); ++i,ipos+=s.colno())
          for (size_t j=s.start(i); j < s.end(i); j++)
          {
             size_t jpos = ipos+j;
             pv[jpos] -= pm[jpos];
          }
    }
    
    template <class T>
    void gmslice_matrix<T>::operator= (const T& e) const
    {
       T *pv = &m(0,0);
    
       for (size_t i=0,ipos=0; i < s.rowno(); ++i,ipos+=s.colno())
          for (size_t j=s.start(i); j < s.end(i); j++)
             pv[ipos+j] = e;
    }
    
    template <class T>
    void gmslice_matrix<T>::operator*= (const T& e) const
    {
       T *pv = &m(0,0);
    
       for (size_t i=0,ipos=0; i < s.rowno(); ++i,ipos+=s.colno())
          for (size_t j=s.start(i); j < s.end(i); j++)
             pv[ipos+j] *= e;
    }
    
    template <class T>
    void gmslice_matrix<T>::operator/= (const T& e) const
    {
       T *pv = &m(0,0);
    
       for (size_t i=0,ipos=0; i < s.rowno(); ++i,ipos+=s.colno())
          for (size_t j=s.start(i); j < s.end(i); j++)
             pv[ipos+j] /= e;
    }
    
    
    //
    // Implementation of matrix class
    //
    
    template <class T> inline
    matrix<T>::matrix ()
    {
        mPtr = allocator( Allocate, 0, 0);
    }
    
    template <class T> inline
    matrix<T>::matrix (size_t nRow, size_t nCol)
    {
        mPtr = allocator( Allocate, nRow, nCol);
    }
    
    template <class T>
    matrix<T>::matrix (size_t nRow, size_t nCol, const T& e)
    {
        mPtr = new base_mat( e, nRow, nCol);
        allocator( Record);
    }
    
    template <class T>
    matrix<T>::matrix (size_t nRow, size_t nCol, const valarray<T>& vec)
    {
       mPtr = new base_mat( nRow, nCol, vec);
       allocator( Record);
    }
    
    template <class T>
    matrix<T>::matrix (size_t row, size_t col, const T* val, MatArrayType aType)
    {
        if (aType == C_ARRAY)
        {
            mPtr = new base_mat( val, row, col);
            allocator( Record);
        }
        else
        {
            mPtr = allocator( Allocate, row, col);
            T *pv = &mPtr->val[0];
    
            for (size_t k=0, j=0; j < col; j++)
                for (size_t i=0; i < row; i++)
                    pv[i*col+j] = val[k++];
        }
    }
    
    template <class T> inline
    matrix<T>::matrix (const matrix<T>& m)
    {
       mPtr = m.mPtr;
       mPtr->refcnt++;
    }
    
    template <class T>
    matrix<T>::matrix (const mslice_matrix<T>& sm)
    {
       const mslice& s = sm.get_slice();
       matrix<T>& m =  sm.get_ref_mat();
       mPtr = allocator( Allocate, s.rowno(), s.colno());
    
       for (size_t i=0; i < s.rowno(); i++)
          for (size_t j=0; j < s.colno(); j++)
             mPtr->val[i*s.colno()+j] = m.mPtr->val[s.start(m.colno())+i*m.colno()+j];
    }
    
    template <class T>
    matrix<T>::matrix (const gmslice_matrix<T>& sm)
    {
       const gmslice& s = sm.get_slice();
       const matrix<T>& m =  sm.get_ref_mat();
       mPtr = allocator( Allocate, s.rowno(), s.colno());
       mPtr->val = T(0);
       T *pv = &mPtr->val[0];
       const T *pm = &m.mPtr->val[0];
    
       for (size_t i=0,ipos=0; i < s.rowno(); ++i,ipos+=s.colno())
          for (size_t j=s.start(i); j < s.end(i); j++)
          {
             size_t jpos = ipos+j;
             pv[jpos] = pm[jpos];
          }
    }
    
    #ifndef _TS_NO_MEMBER_TEMPLATES
    template <class T> template <class X>
    matrix<T>::matrix (const matrix<X>& m)
    {
       mPtr = allocator( Allocate, m.rowno(), m.colno());
       for (size_t i=0; i < rowno(); i++)
          for (size_t j=0; j < colno(); j++)
          {
             X x = m(i,j);
             mPtr->val[i*colno()+j] = x;
          }
    }
    #endif
    
    // Internal copy constructor
    template <class T>
    void matrix<T>::clone ()
    {
       if (mPtr->refcnt == 1)
          return;
       base_mat *mp = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       mp->val = mPtr->val;
       mPtr->refcnt--;
       mPtr = mp;
    }
    
    // Internal matrix allocator/deallocator
    template <class T>
    #if defined(__MWERKS__) || defined(_MSC_VER) || defined(__GNUC__)
    typename
    #endif
    matrix<T>::base_mat*
    matrix<T>::allocator (AllocType bOpt, size_t nrow, size_t ncol)
    {
       static base_mat *mlist[MAX_MATRIX_CACHE+1];
       static int instCount = 0;
       static volatile int semCount = 0;
       static bool first = true;
    
       if (first)
       {
          first = false;
          for (int k=0; k <= MAX_MATRIX_CACHE; k++)
             mlist[k] = NULL;   
       }
    
       switch (bOpt)
       {
       case Record:
          instCount++;
          return NULL;
    
       case Allocate:
          instCount++;
          base_mat *mp;
          if (++semCount == 1)
             for (int i=0; i < MAX_MATRIX_CACHE; i++)
                if (mlist[i] != NULL && mlist[i]->nrow == nrow && mlist[i]->ncol == ncol)
                {
                   mp = mlist[i];
                   mlist[i] = NULL;
                   --semCount;
                   mp->refcnt = 1;
                   return mp;
                }
          --semCount;
          mp = new base_mat( nrow, ncol);
          return mp;
    
       case Free:
          if (--instCount)
          {
             if (++semCount == 1)
                for (int i=0; i < MAX_MATRIX_CACHE; i++)
                {
                   swap( mlist[i], mPtr);
                   if (mPtr == NULL)
                   {
                      --semCount;
                      return NULL;
                   }
                }
             --semCount;
          }
          else                               // clear the cached matrices
          {
             ++semCount;
             for (int i=0; i < MAX_MATRIX_CACHE; i++)
                if (mlist[i] != NULL)
                {
                   delete mlist[i];
                   mlist[i] = NULL;
                }
             --semCount;
          }
          delete mPtr;
          mPtr = NULL;
          break;
       }
       return NULL;
    }
    
    // destructor
    template <class T> inline
    matrix<T>::~matrix ()
    {
       if (--mPtr->refcnt == 0)
         allocator( Free);
    }
    
    template <class T> void
    matrix<T>::resize (size_t nRow, size_t nCol, const T& dval)
    {
       if (mPtr->nrow == nRow && mPtr->ncol == nCol)
          return;
    
       base_mat *mp = allocator( Allocate, nRow, nCol);
    
       for (size_t i=0; i < nRow; i++)
          for (size_t j=0; j < nCol; j++)
             mp->val[i*nCol+j] = (i >= mPtr->nrow || j >= mPtr->ncol) ?
                                      dval : mPtr->val[i*mPtr->ncol+j];
       if (--mPtr->refcnt == 0)
          allocator( Free);
       mPtr = mp;
    }
    
    template <class T> void
    matrix<T>::free ()
    {
       if (--mPtr->refcnt == 0)
          allocator( Free);
       mPtr = allocator( Allocate, 0, 0);
    }
    
    template <class T> inline mslice_matrix<T>
    matrix<T>::operator[] (mslice ms)
    {
       #ifdef RANGE_CHECK_
       if (ms.end( colno()-1) > size())
           throw invalid_argument( "matrix<T>::operator[]: sub-matrix is out of range!");
       #endif
       return mslice_matrix<T>( *this, ms);
    }
    
    template <class T> inline gmslice_matrix<T>
    matrix<T>::operator[] (gmslice gm)
    {
       #ifdef RANGE_CHECK_
       if (gm.size() > size())
           throw invalid_argument( "matrix<T>::operator[]: sub-matrix is out of range!");
       #endif
    
       if (gm.getype() != GENERAL)
       {
          gm.ncol_ = mPtr->ncol;
          gm.srow_.resize( mPtr->nrow);
          gm.rsize_.resize( mPtr->nrow);
       }
       size_t i;
       switch(gm.getype())
       {
       case DIAGONAL:
          gm.rsize_ = 1;
          for (i=0; i < gm.ncol_; ++i)
             gm.srow_[i] = i;
          break;
    
       case TRIDIAGONAL:
          gm.rsize_ = 3;
          gm.rsize_[0] = 2;
          gm.rsize_[gm.ncol_-1] = 2;
          gm.srow_[0] = 0;
          for (i=1; i < gm.ncol_; ++i)
             gm.srow_[i] = i-1;
          break;
    
       case UTRIANG:
          for (i=0; i < gm.ncol_; ++i)
          {
             gm.srow_[i] = i;
             gm.rsize_[i] = gm.ncol_-i;
          }
          break;
    
       case LTRIANG:
          gm.srow_ = 0;
          for (i=0; i < gm.ncol_; ++i)
             gm.rsize_[i] = i+1;
          break;
       }
       return gmslice_matrix<T>( *this, gm);
    }
    
    template <class T> inline matrix<T>&
    matrix<T>::operator= (const matrix<T>& m)
    {
       m.mPtr->refcnt++;
       if (--mPtr->refcnt == 0)
          allocator( Free);
       mPtr = m.mPtr;
       return *this;
    }
    
    template <class T> matrix<T>&
    matrix<T>::operator= (const mslice_matrix<T>& sm)
    {
       const mslice& s = sm.get_slice();
       matrix<T>& m =  sm.get_ref_mat();
       base_mat *mp = allocator( Allocate, s.rowno(), s.colno());
       T *pv = &mp->val[0];
       const T *pm = &m.mPtr->val[0];
    
       for (size_t i=0; i < s.rowno(); i++)
          for (size_t j=0; j < s.colno(); j++)
             pv[i*s.colno()+j] = pm[s.start(m.colno())+i*m.colno()+j];
    
       if (--mPtr->refcnt == 0)
          allocator( Free);
       mPtr = mp;
    
       return *this;
    }
    
    template <class T> matrix<T>&
    matrix<T>::operator= (const gmslice_matrix<T>& sm)
    {
       const gmslice& s = sm.get_slice();
       const matrix<T>& m =  sm.get_ref_mat();
       base_mat *mp = allocator( Allocate, s.rowno(), s.colno());
       mp->val = T(0);
       T *pv = &mp->val[0];
       const T *pm = &m(0,0);
    
       for (size_t i=0,ipos=0; i < s.rowno(); ++i,ipos+=s.colno())
          for (size_t j=s.start(i); j < s.end(i); j++)
          {
             size_t jpos = ipos+j;
             pv[jpos] = pm[jpos];
          }
    
       if (--mPtr->refcnt == 0)
          allocator( Free);
       mPtr = mp;
    
       return *this;
    }
    
    template <class T>
    matrix<T>& matrix<T>::operator= (const T& el)
    {
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       }
       mPtr->val = el;
    
       return *this;
    }
    
    #ifndef _TS_NO_MEMBER_TEMPLATES
    template <class T> template <class X>
    matrix<T>& matrix<T>::operator= (const matrix<X>& m)
    {
       if (--mPtr->refcnt == 0)
          allocator( Free);
    
       mPtr = allocator( Allocate, m.rowno(), m.colno());
       const X *pm = &m(0,0);
       T *pv = &mPtr->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
          pv[i] = pm[i];
    
       return *this;
    }
    #endif
    
    template <class T> const matrix<T>
    matrix<T>::operator[] (mslice ms) const
    {
       #ifdef RANGE_CHECK_
       if (ms.end( colno()-1) > size())
           throw invalid_argument( "matrix<T>::operator[]: sub-matrix is out of range!");
       #endif
    
       matrix<T> tm(ms.rowno(),ms.colno());
       for (size_t i=0; i < tm.rowno(); i++)
          for (size_t j=0; j < tm.colno(); j++)
             tm.mPtr->val[i*tm.colno()+j] = mPtr->val[ms.start(colno())+i*colno()+j];
       return tm;
    }
    
    template <class T> const matrix<T>
    matrix<T>::operator[] (gmslice gm) const
    {
       #ifdef RANGE_CHECK_
       if (gm.size() > size())
           throw invalid_argument( "matrix<T>::operator[]: sub-matrix is out of range!");
       #endif
    
       if (gm.getype() != GENERAL)
       {
          gm.ncol_ = mPtr->ncol;
          gm.srow_.resize( mPtr->nrow);
          gm.rsize_.resize( mPtr->nrow);
       }
       size_t k;
       switch(gm.getype())
       {
       case DIAGONAL:
          gm.rsize_ = 1;
          for (k=0; k < gm.ncol_; ++k)
             gm.srow_[k] = k;
          break;
    
       case TRIDIAGONAL:
          gm.rsize_ = 3;
          gm.rsize_[0] = 2;
          gm.rsize_[gm.ncol_-1] = 2;
          gm.srow_[0] = 0;
          for (k=1; k < gm.ncol_; ++k)
             gm.srow_[k] = k-1;
          break;
    
       case UTRIANG:
          for (k=0; k < gm.ncol_; ++k)
          {
             gm.srow_[k] = k;
             gm.rsize_[k] = gm.ncol_-k;
          }
          break;
    
       case LTRIANG:
          gm.srow_ = 0;
          for (k=0; k < gm.ncol_; ++k)
             gm.rsize_[k] = k+1;
          break;
       }
    
       matrix<T> tm(gm.rowno(), gm.colno(), T(0));
       const T *pv = &mPtr->val[0];
       T *pm = &tm.mPtr->val[0];
    
       for (size_t i=0,ipos=0; i < gm.rowno(); ++i,ipos+=gm.colno())
          for (size_t j=gm.start(i); j < gm.end(i); j++)
          {
             size_t jpos = ipos+j;
             pm[jpos] = pv[jpos];
          }
       return tm;
    }
    
    template <class T> inline Mat_iter<T>
    matrix<T>::row (size_t i)
    {
       #ifdef RANGE_CHECK_
       if (i >= mPtr->nrow)
          throw out_of_range( "matrix<T>::row: subscript is out of range!");
       #endif
       return Mat_iter<T>( *this, slice( i * mPtr->ncol, mPtr->ncol, 1));
    }
    
    template <class T> inline Cmat_iter<T>
    matrix<T>::row (size_t i) const
    {
       #ifdef RANGE_CHECK_
       if (i >= mPtr->nrow)
          throw out_of_range( "matrix<T>::row: subscript is out of range!");
       #endif
       return Cmat_iter<T>( *this, slice( i * mPtr->ncol, mPtr->ncol, 1));
    }
    
    template <class T> inline Mat_iter<T>
    matrix<T>::column (size_t i)
    {
       #ifdef RANGE_CHECK_
       if (i >= mPtr->ncol)
          throw out_of_range( "matrix<T>::column: subscript is out of range!");
       #endif
       return Mat_iter<T>( *this, slice( i, mPtr->nrow, mPtr->ncol));
    }
    
    template <class T> inline Cmat_iter<T>
    matrix<T>::column (size_t i) const
    {
       #ifdef RANGE_CHECK_
       if (i >= mPtr->ncol)
          throw out_of_range( "matrix<T>::column: subscript is out of range!");
       #endif
       return Cmat_iter<T>( *this, slice( i, mPtr->nrow, mPtr->ncol));
    }
    
    template <class T> inline Mat_iter<T>
    matrix<T>::diag (int i)
    {
       #ifdef RANGE_CHECK_
       if (i >= int(mPtr->ncol) || int(mPtr->nrow)+i < 1)
          throw out_of_range( "matrix<T>::diag: subscript is out of range!");
       #endif
       size_t st,sz;
       if (i >= 0)
       {
          st = i;
          sz = mPtr->ncol - i;
          if (sz > mPtr->nrow)
            sz = mPtr->nrow;
       }
       else
       {
          st = -i * int(mPtr->ncol);
          sz = int(mPtr->nrow) + i;
          if (sz > mPtr->ncol)
             sz = mPtr->ncol;
       }
       return Mat_iter<T>( *this, slice( st, sz, mPtr->ncol+1));
    }
    
    template <class T> inline Cmat_iter<T>
    matrix<T>::diag (int i) const
    {
       #ifdef RANGE_CHECK_
       if (i >= int(mPtr->ncol) || int(mPtr->nrow)+i < 1)
          throw out_of_range( "matrix<T>::diag: subscript is out of range!");
       #endif
    
       size_t st,sz;
       if (i >= 0)
       {
          st = i;
          sz = mPtr->ncol - i;
          if (sz > mPtr->nrow)
            sz = mPtr->nrow;
       }
       else
       {
          st = -i * int(mPtr->ncol);
          sz = int(mPtr->nrow) + i;
          if (sz > mPtr->ncol)
             sz = mPtr->ncol;
       }
       return Cmat_iter<T>( *this, slice( st, sz, mPtr->ncol+1));
    }
    
    template <class T> inline Mref<T>
    matrix<T>::operator() (size_t i, size_t j)
    {
       #ifdef RANGE_CHECK_
       if (i >= mPtr->nrow || j >= mPtr->ncol)
          throw out_of_range( "matrix<T>::operator(): subscript is out of range!");
       #endif
       return Mref<T>( *this, i * mPtr->ncol + j);
    }
    
    template <class T> inline const T&
    matrix<T>::operator() (size_t i, size_t j) const
    {
       #ifdef RANGE_CHECK_
       if (i >= mPtr->nrow || j >= mPtr->ncol)
          throw out_of_range( "matrix<T>::operator(): subscript is out of range!");
       #endif
       return mPtr->val[i*mPtr->ncol+j];
    }
    
    template <class T> matrix<T>&
    matrix<T>::operator+= (const matrix<T>& m)
    {
       #ifdef RANGE_CHECK_
       if (rowno() != m.rowno() || colno() != m.colno())
          throw invalid_argument( "matrix<T>::operator+=: Inconsistent matrix size in addition!");
       #endif
    
       base_mat *mp;
       if (shared())
          mp = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       else
          mp = mPtr;
    
       const T *pm = &m(0,0);
       const T *pv = &mPtr->val[0];
       T *pv2 = &mp->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
          pv2[i] = pv[i] + pm[i];
    
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = mp;
       }
       return *this;
    }
    
    template <class T> matrix<T>&
    matrix<T>::operator-= (const matrix<T>& m)
    {
       #ifdef RANGE_CHECK_
       if (rowno() != m.rowno() || colno() != m.colno())
          throw invalid_argument( "matrix<T>::operator-=: Inconsistent matrix size in addition!");
       #endif
    
       base_mat *mp;
       if (shared())
          mp = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       else
          mp = mPtr;
    
       const T *pm = &m(0,0);
       const T *pv = &mPtr->val[0];
       T *pv2 = &mp->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
          pv2[i] = pv[i] - pm[i];
    
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = mp;
       }
       return *this;
    }
    
    template <class T> matrix<T>&
    matrix<T>::operator*= (const T& x)
    {
       base_mat *mp;
    
       if (shared())
          mp = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       else
          mp = mPtr;
    
       const T *pv = &mPtr->val[0];
       T *pv2 = &mp->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
          pv2[i] = pv[i] * x;
    
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = mp;
       }
       return *this;
    }
    
    template <class T> matrix<T>&
    matrix<T>::operator*= (const matrix<T>& m)
    {
       size_t nr1 = mPtr->nrow;
       size_t nc1 = mPtr->ncol;
       size_t nc2 = m.mPtr->ncol;
    
       #ifdef RANGE_CHECK_
       size_t nr2 = m.mPtr->nrow;
       if (nc1 != nr2)
          throw invalid_argument( "matrix<T>::operator*=: Inconsistent matrix size in multiplication!");
       #endif
    
       matrix<T> tm( nr1, nc2);
       const T *pm1 = &mPtr->val[0];
       const T *pm2 = &m(0,0);
       T *pm = &tm.mPtr->val[0];
    
       for (size_t i=0; i < nr1; i++)
          for (size_t j=0; j < nc2; j++)
          {
             pm[i*nc2+j] = T(0);
             for (size_t k=0; k < nc1; k++)
                pm[i*nc2+j] += pm1[i*nc1+k] * pm2[k*nc2+j];
          }
       *this = tm;
       return *this;
    }
    
    template <class T> inline matrix<T>&
    matrix<T>::operator/= (const T& x)
    {
       *this *= 1/x;
       return *this;
    }
    
    template <class T> matrix<T>&
    matrix<T>::operator/= (const matrix<T>& m)
    {
       *this *= !m;
       return *this;
    }
    
    template <class T> matrix<T> 
    matrix<T>::operator- () const
    {
       matrix<T> m( mPtr->nrow, mPtr->ncol);
       T *pv2 = &m.mPtr->val[0];
       const T *pv = &mPtr->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
          pv2[i] = -pv[i];
    
       return m;
    }
    
    template <class T>
    matrix<T> matrix<T>::operator~ () const
    {
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       matrix<T> tm( n, m);
    
       const T *pm = &mPtr->val[0];
       T *ptm = &tm.mPtr->val[0];
    
       for (size_t i=0; i < m; i++)
          for (size_t j=0; j < n; j++)
             ptm[j*m+i] = pm[i*n+j];
    
       return tm;
    }
    
    template <class T>
    inline matrix<T> matrix<T>::operator! () const
    {
       matrix<T> m = *this;
       if (!m.inv())
          m.null();
       return m;
    }
    
    template <class T> inline matrix<T>
    operator+ (const matrix<T>& m1, const matrix<T>& m2)
    {
       matrix<T> m = m1;
       return m += m2;
    }
    
    template <class T> inline matrix<T>
    operator- (const matrix<T>& m1, const matrix<T>& m2)
    {
       matrix<T> m = m1;
       return m -= m2;
    }
    
    template <class T> inline matrix<T>
    operator* (const matrix<T>& m1, const matrix<T>& m2)
    {
       matrix<T> m = m1;
       return m *= m2;
    }
    
    template <class T> inline matrix<T>
    operator* (const matrix<T>& m, const T& el)
    {
       matrix<T> m2 = m;
       return m2 *= el;
    }
    
    template <class T> inline matrix<T>
    operator* (const T& el, const matrix<T>& m)
    {
       return m * el;
    }
    
    template <class T> valarray<T>
    operator* (const matrix<T>& m, const valarray<T>& v)
    {
       #ifdef RANGE_CHECK_
       if (m.rowsize() != v.size())
          throw invalid_argument( "matrix<T>::operator*: Inconsistent matrix size!");
       #endif
    
       valarray<T> vt( m.rowno());
       for (size_t i=0; i < m.rowno(); i++)
       {
          vt[i] = T(0);
          for (size_t j=0; j < m.colno(); j++)
          {
             T el = m(i,j);
             vt[i] += el * v[j];
          }
       }
       return vt;
    }
    
    template <class T> valarray<T>
    operator* (const valarray<T>& v, const matrix<T>& m)
    {
       #ifdef RANGE_CHECK_
       if (v.size() != m.rowno())
          throw invalid_argument( "matrix<T>::operator*: Inconsistent matrix size!");
       #endif
    
       valarray<T> vt( m.colno());
       for (size_t i=0; i < m.colno(); i++)
       {
          vt[i] = T(0);
          for (size_t j=0; j < v.size(); ++j)
          {
             T el = m(j,i);
             vt[i] += el * v[j];
          }
       }
       return vt;
    }
    
    template <class T> inline matrix<T>
    operator/ (const matrix<T>& m1, const matrix<T>& m2)
    {
       matrix<T> m = m1;
       return m *= !m2;
    }
    
    template <class T> inline matrix<T>
    operator/ (const matrix<T>& m, const T& el)
    {
       matrix<T> m2 = m;
       return m2 *= 1/el;
    }
    
    template <class T> inline matrix<T>
    operator/ (const T& el, const matrix<T>& m)
    {
       matrix<T> m2 = !m;
       return m2 *= el;
    }
    
    template <class T> inline valarray<T>
    operator/ (const matrix<T>& m, const valarray<T>& v)
    {
       valarray<T> v2 = T(1) / v;
       return m * v2;
    }
    
    template <class T> inline valarray<T>
    operator/ (const valarray<T>& v, const matrix<T>& m)
    {
       return v * !m;
    }
    
    template <class T> bool
    operator== (const matrix<T>& m1, const matrix<T>& m2)
    {
       if (m1.rowno() != m2.rowno() || m1.colno() != m2.colno())
          return false;
    
       size_t n = m1.size();
       const T *pm1 = &m1(0,0);
       const T *pm2 = &m2(0,0);
    
       for (size_t i=0; i < n; i++)
          if (abs( pm1[i] - pm2[i]) > epsilon( pm1[i]))
             return false;
    
       return true;
    }
    
    template <class T> inline bool
    operator!= (const matrix<T>& m1, const matrix<T>& m2)
    {
       return !(m1 == m2);
    }
    
    template <class T> bool 
    isVecEq (const valarray<T>& v1, const valarray<T>& v2)
    {
       if (v1.size() != v2.size())
          return false;
    
       for (size_t i=0; i < v1.size(); i++)
       {
          T d = v1[i] - v2[i];
          if (abs( d) > epsilon( v1[i]))
             return false;
       }
       return true;
    }
    
    template <class T> inline void 
    matrix<T>::unit ()
    {
       null();
       diag(0) = T(1);
    }
    
    template <class T> inline void 
    matrix<T>::null ()
    {
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       }
       mPtr->val = T(0);
    }
    
    inline int rand_ () { return rand(); }
    
    template <class T> void
    matrix<T>::rand (int rmin, int rmax, int rseed)
    {
       if (!rmin && !rmax)
       {
          rmin = -1;
          rmax = 1;
       }
       if (rmin > rmax)
       {
          int tmp = rmin;
          rmin = rmax;
          rmax = tmp;
       }
       int rdiff = rmax - rmin;
    
       if (rseed)
          srand( rseed);
    
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       }
    
       T *pv = &mPtr->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
       {
          int rn = rand_();
          if (rdiff > 1)
          {
             rn %= rdiff;
             pv[i] = T(rn) - T(rdiff) / T(2);
             if (rdiff > 1)
             {
                T f = ((T)rand_()) / T(RAND_MAX);
                if (rn > rdiff)
                   pv[i] -= f;
                else
                   pv[i] += f;
             }
          }
          else
             pv[i] = T(rn) / T(RAND_MAX);
       }
    }
    
    template <class T> matrix<T>
    pow (const matrix<T>& m, size_t n)
    {
       #ifdef RANGE_CHECK_
       if (m.colno() != m.rowno() || m.rowno() == 0)
          throw invalid_argument( "pow(): Non-square matrix!");
       #endif
    
       matrix<T> x(m.rowno(), m.colno());
       if (n == 0)
       {
          x.unit();
          return x;
       }
       matrix<T> m1 = m;
       bool first = true;
    
       while (1)
       {
         if (n & 1)
         {
           if (first)
           {
             x = m1;
             first = false;
           }
           else
             x *= m1;
         }
         n >>= 1;
         if (n == 0)
            break;
    
         m1 *= m1;
       }
       return x;
    }
    
    template <class T> matrix<T>
    abs (const matrix<T>& m)
    {
       matrix<T> a(m.rowno(), m.colno());
       T *pa = &a(0,0);
       const T *pm = &m(0,0);
       size_t n = m.size();
       for (size_t i=0; i < n; i++)
          pa[i] = T(abs( pm[i]));
       return a;
    }
    
    template <class T> matrix<T>
    floor (const matrix<T>& m)
    {
       matrix<T> a(m.rowno(), m.colno());
       T *pa = &a(0,0);
       const T *pm = &m(0,0);
       size_t n = m.size();
       for (size_t i=0; i < n; i++)
          pa[i] = T(floor( pm[i]));
       return a;
    }
    
    template <class T> matrix<T>
    ceil (const matrix<T>& m)
    {
       matrix<T> a(m.rowno(), m.colno());
       T *pa = &a(0,0);
       const T *pm = &m(0,0);
       size_t n = m.size();
       for (size_t i=0; i < n; i++)
          pa[i] = T(ceil( pm[i]));
       return a;
    }
    
    template <class T> inline T
    matrix<T>::trace (int i) const
    {
       valarray<T> d = diag( i);
       return d.sum();
    }
    
    ///
    // Computes (a^2 + b^2 )^1/2 without destructive underflow or overflow.
    //
    template <class T> inline 
    T hypot (T a, T b)
    {
       a = abs( a);
       b = abs( b);
       if (a > b)
       {
          T c = b/a;
          return (a * sqrt( T(1) + c * c));
       }
       if (b < epsilon(b))
          return T(0);
       T c = a/b;
       return (b * sqrt( T(1) + c * c));
    }
    
    template <class T> inline
    T sign (T a, T b)
    {
       return (b >= T(0) ? abs(a) : -abs(a));
    }
    
    template <class T>
    bool matrix<T>::lud (valarray<size_t>& ri, T* pDet)
    {
       size_t i,j,k;
       double ta,tb;
    
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          throw invalid_argument( "matrix<T>::lud: Non-square matrix!");
       #endif
    
       if (mPtr->nrow != ri.size())
          ri.resize( mPtr->nrow);
    
       if (shared())
          clone();
    
       if (pDet != NULL)
          *pDet = T(1);
    
       size_t n = mPtr->nrow;
       T *pv = &mPtr->val[0];
    
       for (i=0; i < n; i++)
          ri[i] = i;
    
       for (k=0; k < n-1; k++)
       {
          j = k;
          ta = abs( pv[ri[k]*n+k]);
          for (i=k+1; i < n; i++)
             if ((tb = abs( pv[ri[i]*n+k])) > ta)
             {
                ta = tb;
                j = i;
             }
          if (j != k)
          {
             swap( ri[j], ri[k]);
             if (pDet != NULL)
                *pDet = - *pDet;
          }
          size_t kpos = ri[k] * n;
    
          if (abs( pv[kpos+k]) < epsilon( pv[kpos+k]))
             return false;
    
          if (pDet != NULL)
             *pDet *= pv[kpos+k];
    
          for (i=k+1; i < n; i++)
          {
             size_t ipos = ri[i] * n;
             T a = pv[ipos+k] /= pv[kpos+k];
    
             for (j=k+1; j < n; j++)
                pv[ipos+j] -= a * pv[kpos+j];
          }
       }
       if (pDet != NULL)
          *pDet *= pv[ri[k]*n+k];
    
       return true;
    }
    
    template <class T> void
    matrix<T>::lubksb (const valarray<size_t>& ri, const valarray<T>& v, valarray<T>& s) const
    {
       size_t i,j,k,ip;
       bool nonzero = false;
       size_t n = mPtr->nrow;
       const T *pv = &mPtr->val[0];
    
       #ifdef RANGE_CHECK_
       if (n != v.size())
          throw invalid_argument( "matrix<T>::lubksb: Invalid vector size!");
       #endif
    
       if (n != s.size())
          s.resize( n);
    
       for (k=0,i=0; i < n; i++)
       {
          ip = ri[i];
          s[i] = v[ip];
    
          if (nonzero)
          {
             ip *= n;
             for (j=k; j < i; j++)
                s[i] -= pv[ip+j] * s[j];
          }
          else if (s[i] != T(0))
          {
             k = i;
             nonzero = true;
          }
       }
       for (i=n-1; ; i--)
       {
          ip = ri[i] * n;
          for (j=i+1; j < n; j++)
             s[i] -= pv[ip+j] * s[j];
          s[i] /= pv[ip+i];
          if (i == 0) break;
       }
    }
    
    template <class T> void 
    matrix<T>::lumpove (const matrix<T>& ludm, const valarray<size_t>& ri, const valarray<T>& v, valarray<T>& s) const
    {
       size_t i,j;
       size_t n = mPtr->ncol;
       const T *pa = &mPtr->val[0];
    
       valarray<T> rv(n), rs(n);
       for (i=0; i < n; i++)
       {
          T tmp = -v[i];
          for (j=0; j < n; j++)
             tmp += pa[i*n+j] * s[j];
          rv[i] = tmp;
       }
       ludm.lubksb( ri, rv, rs);
       for (i=0; i < n; i++)
          s[i] -= rs[i];
    }
    
    template <class T> bool 
    matrix<T>::solve (const valarray<T>& v, valarray<T>& s) const
    {
       matrix<T> tm = *this;
       valarray<size_t> ri( mPtr->nrow);
    
       if (tm.lud( ri))
       {
          tm.lubksb( ri, v, s);
          return true;
       }
       return false;
    }
    
    template <class T> bool 
    matrix<T>::solve_chol (const valarray<T>& v, valarray<T>& s) const
    {
       matrix<T> tm = *this;
    
       if (tm.chold())
       {
          tm.cholbksb( v, s);
          return true;
       }
       return false;
    }
    
    template <class T> bool 
    matrix<T>::solve_sv (const valarray<T>& v, valarray<T>& s) const
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow < mPtr->ncol)
          throw invalid_argument( "matrix<T>::solve_sv(): No. of rows must be equal or greater than column!");
       #endif
    
       matrix<T> vc( mPtr->ncol, mPtr->ncol);
       matrix<T> tm = *this;
       valarray<T> w( mPtr->ncol);
    
       if (tm.svd( vc, w))
       {
          tm.svbksb( vc, w, v, s);
          return true;
       }
       return false;
    }
    
    template <class T> bool 
    matrix<T>::solve_qr (const valarray<T>& v, valarray<T>& s) const
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow < mPtr->ncol)
          throw invalid_argument( "matrix<T>::solve_qr(): No. of rows must be equal or greater than column!");
       #endif
    
       matrix<T> r( mPtr->ncol, mPtr->ncol);
       matrix<T> tm = *this;
    
       tm.qrd( r);
       tm.qrbksb( r, v, s);
    
       return true;
    }
    
    template <class T> bool
    matrix<T>::solve (const matrix<T>& v, matrix<T>& s) const
    {
       size_t n = mPtr->nrow;
       #ifdef RANGE_CHECK_
       if (n != v.rowno())
          throw invalid_argument( "matrix<T>::operator*: Inconsistent number of vectors in solve!");
       #endif
    
       if (v.size() != s.size())
          s.resize( v.rowno(), v.colno());
    
       matrix<T> tm = *this;
       valarray<size_t> ri(n);
       valarray<T> vt(n), st(n);
    
       if (!tm.lud( ri))
          return false;
    
       for (size_t i=0; i < v.colno(); i++)
       {
          vt = v.column(i);
          tm.lubksb( ri, vt, st);
          s.column(i) = st;
       }
       return true;
    }
    
    template <class T> bool 
    matrix<T>::chold ()
    {
       #ifdef RANGE_CHECK_
       if (!isSymmetric())
          throw invalid_argument( "matrix<T>::chold: Non-symmetric matrix!");
       #endif
    
       int n = int(mPtr->nrow);
       if (shared())
          clone();
       T *pa = &mPtr->val[0];
       int i,j;
    
       for (i=0; i < n; i++)
          for (j=i; j < n; j++)
          {
             T s = pa[i*n+j];
             for (int k=i-1; k >= 0; k--)
                s -= pa[i*n+k] * pa[j*n+k];
    
             if (i == j)
             {
                if (s  <= epsilon(s))
                   return false;
                pa[i*n+i] = sqrt( s);
             }
             else
                pa[j*n+i] = s / pa[i*n+i];
          }
        for (j=1; j < n; j++)
            for (i=0; i < j; i++)
                pa[i*n+j] = T(0);
    
       return true;
    }
    
    template <class T> void
    matrix<T>::cholbksb (const valarray<T>& v, valarray<T>& s) const
    {
       int i,k,n;
       T a;
    
       const T *pm = &mPtr->val[0];
       n = int(mPtr->nrow);
       i = int(s.size());
    
       if (n != i)
          s.resize( n);
    
       for (i=0; i < n; i++)
       {
          a = v[i];
          for (k=i-1; k >= 0; k--)
             a -= pm[i*n+k] * s[k];
          s[i] = a / pm[i*n+i];
       }
       for (i=n-1; i >= 0; i--)
       {
          a = s[i];
          for (k=i+1; k < n; k++)
             a -= pm[k*n+i] * s[k];
          s[i] = a / pm[i*n+i];
       }
    }
    
    template <class T> void        
    matrix<T>::qrd (matrix<T>& r)
    {
       size_t i,j,k;
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
    
       #ifdef RANGE_CHECK_
       if (m < n)
          throw invalid_argument( "matrix<T>::qrd: Inconsistent matrix size!");
       #endif
    
       if (shared())
          clone();
    
       if (n != r.rowno() || r.rowno() != r.colno())
          r.resize( n, n);
    
       T *pa = &mPtr->val[0];
       T *pr = &r(0,0);
    
       for (k=0; k < n; k++)
       {
          T nrm(0);
          for (i = k; i < m; i++)
             nrm = hypot( nrm, pa[i*n+k]);
    
          if (nrm > epsilon(nrm))
          {
             if (pa[k*n+k] < T(0))
                nrm = -nrm;
    
             for (i = k; i < m; i++)
                pa[i*n+k] /= nrm;
             pa[k*n+k] += T(1);
    
             for (j=k+1; j < n; j++)
             {
                T s(0);
                for (i=k; i < m; i++)
                   s += pa[i*n+k] * pa[i*n+j];
    
                s /= -pa[k*n+k];
                for (i=k; i < m; i++)
                   pa[i*n+j] += s * pa[i*n+k];
             }
          }
          pr[k*n+k] = -nrm;
       }
    
        for (j=1; j < n; j++)
            for (size_t i=0; i < j; i++)
          {
             pr[i*n+j] = pa[i*n+j];
             pr[j*n+i] = T(0);
          }
    
       matrix<T> q(m,n);
       T *pq = &q(0,0);
    
       for (k=n-1; ; k--)
       {
          for (i = 0; i < m; i++)
             pq[i*n+k] = T(0);
    
          pq[k*n+k] = T(1);
          for (j=k; j < n; j++)
          {
             if (pa[k*n+k] != T(0))
             {
                T s(0);
                for (i=k; i < m; i++)
                   s += pa[i*n+k] * pq[i*n+j];
                s /= -pa[k*n+k];
                for (i = k; i < m; i++)
                   pq[i*n+j] += s * pa[i*n+k];
             }
          }
          if (k == 0)
             break;
       }
       *this = q;
    }
    
    template <class T> void
    matrix<T>::qrbksb (const matrix<T>& r, const valarray<T>& v, valarray<T>& s) const
    {
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       size_t i,j;
    
       #ifdef RANGE_CHECK_
       if (m != v.size())
          throw std::invalid_argument( "matrix<T>::qrsolve: Inconsistent vector size!");
       #endif
    
       if (n != s.size())
          s.resize( n);
    
       for (i=0; i < n; i++)
       {
          s[i] = T(0);
          for (j=0; j < m; j++)
             s[i] += mPtr->val[j*n+i] * v[j];
       }
    
       const T *pr = &r(0,0);
       for (i=n-1; ; i--)
       {
          size_t ip = i * n;
          for (j=i+1; j < n; j++)
             s[i] -= pr[ip+j] * s[j];
          s[i] /= pr[ip+i];
          if (i == 0) break;
       }
    }
    
    template <class T> bool 
    matrix<T>::svd (matrix<T>& vc, valarray<T>& w)
    {
       size_t flag,i,its,j,jj,k,l,nm;
       T c,f,h,s,x,y,z,tmp;
    
       if (shared())
          clone();
    
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
    
       if (vc.rowno() != n || vc.colno() != n)
          vc.resize( n, n);
       if (w.size() != n)
          w.resize( n);
    
       T *a = &mPtr->val[0];
       T *v = &vc(0,0);
       valarray<T> rv1(n);
    
       T g(0), scale(0), anorm(0);
       for (i=0; i < n; i++)
       {
          l = i + 1;
          rv1[i] = scale * g;
          g = s = scale = T(0);
          if (i < m)
          {
             for (k=i; k < m; k++)
                scale += abs( a[k*n+i]);
    
             if (scale > epsilon( scale))
             {
                for (k=i; k < m; k++)
                {
                   tmp = a[k*n+i] /= scale;
                   s += tmp * tmp;
                }
                f = a[i*n+i];
                g = -sign( sqrt(s), f);
                h = f * g - s;
                a[i*n+i] = f - g;
    
                for (j=l; j < n; j++)
                {
                   for (s=T(0),k=i; k < m; k++)
                      s += a[k*n+i] * a[k*n+j];
                   f = s / h;
                   for (k=i; k < m; k++)
                      a[k*n+j] += f * a[k*n+i];
                }
                for (k=i; k < m; k++)
                   a[k*n+i] *= scale;
             }
          }
          w[i] = scale * g;
          g = s = scale = T(0);
          if (i < m && i != n-1)
          {
             for (k=l; k < n; k++)
                scale += abs(a[i*n+k]);
    
             if (scale > epsilon(scale))
             {
                for (k=l; k < n; k++)
                {
                   tmp = a[i*n+k] /= scale;
                   s += tmp * tmp;
                }
                f = a[i*n+l];
                g = -sign( sqrt(s), f);
                h = f * g - s;
                a[i*n+l] = f - g;
    
                for (k=l; k < n; k++)
                   rv1[k] = a[i*n+k] / h;
    
                for (j=l; j < m; j++)
                {
                   for (s=T(0),k=l; k < n; k++)
                      s += a[j*n+k] * a[i*n+k];
                   for (k=l; k < n; k++)
                      a[j*n+k] += s * rv1[k];
                }
                for (k=l; k < n; k++)
                   a[i*n+k] *= scale;
             }
          }
          anorm = techsoft::max( anorm, abs( w[i]) + abs( rv1[i]) );
       }
    
       for (i=n-1; ; i--)
       {
          if (i < n-1)
          {
             if (abs(g) > epsilon(g))
             {
                for (j=l; j < n; j++)
                    v[j*n+i] = (a[i*n+j] / a[i*n+l]) / g;
                for (j=l; j < n; j++)
                {
                   for (s=T(0),k=l; k < n; k++)
                      s += a[i*n+k] * v[k*n+j];
    
                   for (k=l; k < n; k++)
                      v[k*n+j] += s * v[k*n+i];
                }
             }
             for (j=l; j < n; j++)
                v[i*n+j] = v[j*n+i] = T(0);
          }
          v[i*n+i] = T(1);
          g = rv1[i];
          l = i;
          if (i == 0)
             break;
       }
    
       for (i=techsoft::min(m,n)-1; ; i--)
       {
          l = i + 1;
          g = w[i];
          for (j=l; j < n; j++)
             a[i*n+j] = T(0);
          if (abs(g) > epsilon(g))
          {
             g = T(1) / g;
             for (j=l; j < n; j++)
             {
                for (s=T(0),k=l; k < m; k++)
                   s += a[k*n+i] * a[k*n+j];
    
                f = (s / a[i*n+i]) * g;
    
                for (k=i; k < m; k++)
                   a[k*n+j] += f * a[k*n+i];
             }
             for (j=i; j < m; j++)
                a[j*n+i] *= g;
    
          }
          else
          {
              for (j=i; j < m; j++)
                  a[j*n+i] = T(0);
          }
          ++a[i*n+i];
          if (i == 0)
             break;
       }
    
       for (k=n-1; ; k--)
       {
          for (its=1; its <= 30; its++)
          {
             flag = 1;
             for (l=k; ; l--)
             {
                nm = l-1;
                if (abs( rv1[l]) < epsilon(rv1[l]))
                {
                   flag = 0;
                   break;
                }
                if (abs( w[nm]) < epsilon(w[nm]))
                   break;
                if (l == 0)
                   break;
             }
             if (flag)
             {
                c = T(0);
                s = T(1);
                for (i=l; i <= k; i++)
                {
                   f = s * rv1[i];
                   rv1[i] = c * rv1[i];
                   if (abs(f) < epsilon(f))
                      break;
                   g = w[i];
                   h = hypot( f, g);
                   w[i] = h;
                   h = T(1) / h;
                   c = g * h;
                   s = -f * h;
                   for (j=0; j < m; j++)
                   {
                      y = a[j*n+nm];
                      z = a[j*n+i];
                      a[j*n+nm] = y * c + z * s;
                      a[j*n+i] = z * c - y * s;
                   }
                }
             }
             z = w[k];
             if (l == k)                 
             {
                if (z < T(0))
                {
                   w[k] = -z;
                   for (j=0; j < n; j++)
                      v[j*n+k] = -v[j*n+k];
                }
                break;
             }
    
             if (its == 30)
                return false;
    
             x = w[l];
             nm = k - 1;
             y = w[nm];
             g = rv1[nm];
             h = rv1[k];
             f = ((y-z) * (y+z) + (g-h) * (g+h)) / (T(2) * h * y);
             g = hypot( f, T(1));
             f = ((x-z) * (x+z) + h * ((y / (f+sign( g, f))) - h)) / x;
             c = s = T(1);                    
             for (j=l; j <= nm; j++)
             {
                i = j + 1;
                g = rv1[i];
                y = w[i];
                h = s * g;
                g = c * g;
                z = hypot( f, h);
                rv1[j] = z;
                c = f / z;
                s = h / z;
                f = x * c + g * s;
                g = g*c - x*s;
                h = y * s;
                y *= c;
                for (jj=0; jj < n; jj++)
                {
                   x = v[jj*n+j];
                   z = v[jj*n+i];
                   v[jj*n+j] = x * c + z * s;
                   v[jj*n+i] = z * c - x * s;
                }
                z = hypot( f, h);
                w[j] = z;
                if (abs(z) > epsilon(z))
                {
                   z = 1.0 / z;
                   c = f * z;
                   s = h * z;
                }
                f = c * g + s * y;
                x = c * y - s * g;
                for (jj=0; jj < m; jj++)
                {
                   y = a[jj*n+j];
                   z = a[jj*n+i];
                   a[jj*n+j] = y * c + z * s;
                   a[jj*n+i] = z * c - y * s;
                }
             }
             rv1[l] = 0.0;
             rv1[k] = f;
             w[k] = x;
          }
          if (k == 0)
             break;
       }
       return true;
    }
    
    template <class T> void
    matrix<T>::svbksb (const matrix<T>& vc, const valarray<T>& w, const valarray<T>& b, valarray<T>& s) const
    {
       size_t i,j,k;
       const T *pm;
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       valarray<T> tmp(n);
    
       pm = &mPtr->val[0];
       for (j=0; j < n; j++)                  
       { 
          T ta(0);
          if (abs(w[j]) > epsilon(w[j]))
          { 
             for (i=0; i < m; i++) 
                ta += pm[i*n+j] * b[i];
             ta /= w[j];                       
          }
          tmp[j] = ta;
       }
    
       pm = &vc.mPtr->val[0];
       for (j=0; j < n; j++)                  
       { 
          T ta(0);
          for (k=0; k < n; k++) 
             ta += pm[j*n+k] * tmp[k];
          s[j] = ta;
       }
    }
    
    template <class T> T
    matrix<T>::norm1 () const
    {
       T nr(0);
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       const T *pm = &mPtr->val[0];
    
       for (size_t j=0; j < n; j++)
       {
          T s(0);
          for (size_t i=0; i < m; i++)
             s += abs( pm[i*n+j]);
          nr = techsoft::max( nr, s);
       }
       return nr;
    }
    
    template <class T> T
    matrix<T>::norm2 () const
    {
       matrix<T> m(*this);
       matrix<T> v(mPtr->ncol,mPtr->ncol);
       valarray<T> w(mPtr->ncol);
       if (m.svd( v, w))
          return w.max();
       return T(0);
    }
    
    template <class T> T
    matrix<T>::normI () const
    {
       T nr(0);
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       const T *pm = &mPtr->val[0];
    
       for (size_t i=0; i < m; i++)
       {
          T s(0);
          for (size_t j=0; j < n; j++)
             s += abs( pm[i*n+j]);
          nr = techsoft::max( nr, s);
       }
       return nr;
    }
    
    template <class T> T
    matrix<T>::normF () const
    {
       T nr(0);
       size_t n = size();
       const T *pm = &mPtr->val[0];
    
       for (size_t i=0; i < n; i++)
          nr += pm[i] * pm[i];
       nr = sqrt( nr);
    
       return nr;
    }
    
    template <class T> T
    matrix<T>::cond () const
    {
       matrix<T> m(*this);
       matrix<T> v(mPtr->ncol,mPtr->ncol);
       valarray<T> w(mPtr->ncol);
       if (m.svd( v, w))
          return w.max()/w.min();
       return T(0);
    }
    
    template <class T> size_t 
    matrix<T>::rank () const
    {
       matrix<T> m(*this);
       matrix<T> v(mPtr->ncol,mPtr->ncol);
       valarray<T> w(mPtr->ncol);
    
       if (!m.svd( v, w))
          return 0;
    
       size_t r = 0;
       for (size_t i=0; i < mPtr->ncol; i++)
          if (abs(w[i]) > epsilon( w[i]))
             ++r;
       return r;
    }
    
    template <class T> T
    matrix<T>::det () const
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          throw invalid_argument( "matrix<T>::det: Determinant of a non-square matrix!");
       #endif
    
       T d;
    
       matrix<T> m(*this);
       valarray<size_t> ri( mPtr->nrow);
       if (!m.lud( ri, &d))
          d = T(0);
    
       return d;
    }
    
    template <class T> T
    matrix<T>::cofact (size_t row, size_t col) const
    {
       size_t i,i1,j,j1;
    
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          throw invalid_argument( "matrix<T>::Cofact(): Cofactor of a non-square matrix!");
    
       if (row >= mPtr->nrow || col >= mPtr->ncol)
          throw out_of_range( "matrix<T>::Cofact(): Index out of range!");
       #endif
    
       size_t n = mPtr->nrow;
       matrix<T> tm (n-1,n-1);
       const T *pm = &mPtr->val[0];
       T *ptm = &tm.mPtr->val[0];
    
       for (i=0,i1=0; i < n; i++)
       {
          if (i == row)
            continue;
          for (j=0,j1=0; j < n; j++)
          {
             if (j == col)
                continue;
             ptm[i1*(n-1)+j1] = pm[i*n+j];
             j1++;
          }
          i1++;
       }
       T cof = tm.det();
       if ((row+col)%2 == 1)
          cof = -cof;
    
       return cof;
    }
    
    template <class T> matrix<T>
    matrix<T>::adj () const
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          throw invalid_argument( "matrix<T>::adj(): Adjoin of a non-square matrix!");
       #endif
    
       size_t n = mPtr->nrow;
       matrix<T> tm(n,n);
       T *pm = &tm.mPtr->val[0];
    
       for (size_t i=0; i < n; i++)
          for (size_t j=0; j < n; j++)
             pm[j*n+i] = cofact(i,j);
       return tm;
    }
    
    template <class T>
    bool matrix<T>::inv ()
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          throw invalid_argument( "matrix<T>::inv(): Inversion of non-square matrix!");
       #endif
    
       if (shared())
          clone();
    
       size_t n = mPtr->nrow;
       T *pv = &mPtr->val[0];
       size_t i,j,k,ipos,kpos;
    
       valarray<size_t> ri(n);
       for (i=0; i < n; i++)
          ri[i] = i;
    
       for (k=0; k < n; k++)
       {
          double ta,tb;
          T a(0);
    
          kpos = k * n;
          i = k;
          ta = abs( pv[kpos+k]);
          for (j=k+1; j < n; j++)
             if ((tb = abs(pv[j*n+k])) > ta)
             {
                ta = tb;
                i = j;
             }
    
          if (ta < epsilon( a))
             return false;
    
          if (i != k)
          {
             swap( ri[k], ri[i]);
             for (ipos=i*n,j=0; j < n; j++)
                swap( pv[kpos+j], pv[ipos+j]);
          }
    
          a = T(1) / pv[kpos+k];
          pv[kpos+k] = T(1);
    
          for (j=0; j < n; j++)
             pv[kpos+j] *= a;
    
          for (i=0; i < n; i++)
          {
             if (i != k)
             {
                ipos = i * n;
                a = pv[ipos+k];
                pv[ipos+k] = T(0);
                for (j=0; j < n; j++)
                   pv[ipos+j] -= a * pv[kpos+j];
             }
          }
       }
       for (j=0; j < n; j++)
       {
          if (j != ri[j])         // Column is out of order
          {
             k = j+1;
             while (j != ri[k])
                 k++;
             for (i=0; i < n; i++)
                swap( pv[i*n+j], pv[i*n+k]);
             swap( ri[j], ri[k]);
          }
       }
       return true;
    }
    
    template <class T>
    bool matrix<T>::inv_lu ()
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          throw invalid_argument( "matrix<T>::inv_lu(): Inversion of non-square matrix!");
       #endif
    
       size_t n = mPtr->nrow;
       size_t i,j,k;
    
       valarray<size_t> ri(n);
       valarray<T> v(n), s(n);
    
       if (!lud( ri))
          return false;
    
       matrix<T> tm(n,n);
       T *pv = &tm.mPtr->val[0];
       for (j=0; j < n; j++)
       {
          for (i=0; i < n; i++)
             v[i] = T(0);
          v[j] = T(1);
          lubksb( ri, v, s);
          for (k=0,i=0; i < n; i++,k+=n)
             pv[k+j] = s[i];
       }
       *this = tm;
       return true;
    }
    
    template <class T>
    bool matrix<T>::inv_qr ()
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          throw invalid_argument( "matrix<T>::inv_qr(): Inversion of non-square matrix!");
       #endif
    
       size_t n = mPtr->nrow;
       size_t i,j,k;
    
       valarray<T> v(n), s(n);
       matrix<T> tm(n,n), r(n,n);
    
       qrd( r);
       T *pv = &tm.mPtr->val[0];
       for (j=0; j < n; j++)
       {
          for (i=0; i < n; i++)
             v[i] = T(0);
          v[j] = T(1);
          qrbksb( r, v, s);
          for (k=0,i=0; i < n; i++,k+=n)
             pv[k+j] = s[i];
       }
       *this = tm;
       return true;
    }
    
    template <class T>
    bool matrix<T>::inv_sv ()
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          throw invalid_argument( "matrix<T>::inv2(): Inversion of non-square matrix!");
       #endif
    
       size_t n = mPtr->ncol;
       size_t i,j,k;
    
       valarray<T> w(n), v(n), s(n);
       matrix<T> vc(n,n);
    
       if (!svd( vc, w))
          return false;
    
       matrix<T> tm(n,n);
       T *pv = &tm.mPtr->val[0];
       for (j=0; j < n; j++)
       {
          for (i=0; i < n; i++)
             v[i] = T(0);
          v[j] = T(1);
          svbksb( vc, w, v, s);
          for (k=0,i=0; i < n; i++,k+=n)
             pv[k+j] = s[i];
       }
       *this = tm;
       return true;
    }
    
    template <class T> inline bool
    matrix<T>::isSingular () const
    {
       if (mPtr->nrow != mPtr->ncol)
          return false;
       return (abs( det()) < epsilon(T(0)));
    }
    
    template <class T> bool
    matrix<T>::isDiagonal () const
    {
       if (mPtr->nrow != mPtr->ncol)
          return false;
    
       size_t n = mPtr->nrow;
       const T *pm = &mPtr->val[0];
    
       for (size_t i=0; i < n; i++)
          for (size_t j=0; j < n; j++)
             if (i != j && abs( pm[i*n+j]) > epsilon( pm[i]))
                return false;
       return true;
    }
    
    template <class T> bool
    matrix<T>::isScalar () const
    {
       if (!isDiagonal())
         return false;
    
       T v = mPtr->val[0];
       size_t n = mPtr->nrow;
       const T *pm = &mPtr->val[0];
    
       for (size_t i=1; i < n; i++)
         if (abs( pm[i*n+i] - v) > epsilon( v))
            return false;
    
       return true;
    }
    
    template <class T> bool
    matrix<T>::isUnit () const throw()
    {
       if (mPtr->nrow != mPtr->ncol)
          return false;
    
       size_t n = mPtr->nrow;
       const T *pm = &mPtr->val[0];
    
       for (size_t k=0; k < n; k++)
         if (abs( pm[k*n+k] - T(1)) > epsilon( T(1)))
           return false;
    
        for (size_t j=1; j < n; j++)
          for (size_t i=0; i < j; i++)
           if (abs( pm[i*n+j]) > epsilon( T(1)) || abs( pm[j*n+i]) > epsilon( T(1)))
              return false;
    
       return true;
    }
    
    template <class T> bool
    matrix<T>::isNull () const throw()
    {
       size_t n = size();
       const T *pm = &mPtr->val[0];
    
       for (size_t i=0; i < n; i++)
          if (abs( pm[i]) > epsilon( pm[i]))
             return false;
    
       return true;
    }
    
    template <class T> bool
    matrix<T>::isSymmetric () const
    {
       if (mPtr->nrow != mPtr->ncol)
          return false;
    
       size_t n = mPtr->nrow;
       const T *pm = &mPtr->val[0];
    
        for (size_t j=1; j < n; j++)
            for (size_t i=0; i < j; i++)
                if (abs( pm[i*n+j] - pm[j*n+i]) > epsilon( pm[0]))
                return false;
       return true;
    }
    
    template <class T> bool
    matrix<T>::isSkewSymmetric () const
    {
       if (mPtr->nrow != mPtr->ncol)
          return false;
    
       size_t n = mPtr->nrow;
       const T *pm = &mPtr->val[0];
    
        for (size_t j=1; j < n; j++)
            for (size_t i=0; i < j; i++)
                if (abs( pm[i*n+j] + pm[j*n+i]) > epsilon( pm[0]))
                return false;
       return true;
    }
    
    template <class T> bool
    matrix<T>::isUpperTriangular () const
    {
       if (mPtr->nrow != mPtr->ncol)
          return false;
    
       size_t n = mPtr->nrow;
       const T *pm = &mPtr->val[0];
    
       for (size_t i=1; i < n; i++)
          for (size_t j=0; j < i-1; j++)
             if (abs( pm[i*n+j]) > epsilon( pm[0]))
                return false;
       return true;
    }
    
    template <class T> bool
    matrix<T>::isLowerTriangular () const
    {
       if (mPtr->nrow != mPtr->ncol)
          return false;
    
       size_t n = mPtr->nrow;
       const T *pm = &mPtr->val[0];
    
       for (size_t j=1; j < n; j++)
          for (size_t i=0; i < j-1; i++)
             if (abs( pm[i*n+j]) > epsilon( pm[0]))
                return false;
       return true;
    }
    
    template <class T> inline bool
    matrix<T>::isRowOrthogonal () const
    {
       matrix<T> tm = *this * ~(*this);
       return tm.isUnit();
    }
    
    template <class T> inline bool
    matrix<T>::isColOrthogonal () const
    {
       matrix<T> tm = ~(*this) * (*this);
       return tm.isUnit();
    }
    
    template <class T> matrix<T> 
    matrix<T>::apply (T (*fn)(T)) const
    {
       matrix<T> a( mPtr->nrow, mPtr->ncol);
       size_t n = a.size();
    
       const T *pm = &mPtr->val[0];
       T *pa = &a.mPtr->val[0];
    
       for (size_t i=0; i < n; i++)
          pa[i] = fn( pm[i]);
    
       return a;
    }
    
    template <class T> matrix<T> 
    matrix<T>::apply (T (*fn)(size_t,size_t,T)) const
    {
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       matrix<T> a( m, n);
    
       const T *pm = &mPtr->val[0];
       T *pa = &a.mPtr->val[0];
       size_t i,j,im,jn;
    
       for (i=0,im=0; i < m; im+=n,i++)
          for (j=0,jn=im; j < n; jn++,j++)
             pa[jn] = fn( i, j, pm[jn]);
    
       return a;
    }
    
    template <class T> matrix<T> 
    matrix<T>::apply (T (*fn)(size_t,size_t,const T&)) const
    {
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       matrix<T> a( m, n);
    
       const T *pm = &mPtr->val[0];
       T *pa = &a.mPtr->val[0];
       size_t i,j,im,jn;
    
       for (i=0,im=0; i < m; im+=n,i++)
          for (j=0,jn=im; j < n; jn++,j++)
             pa[jn] = fn( i, j, pm[jn]);
    
       return a;
    }
    
    template <class T> inline matrix<T> 
    matrix<T>::apply (T (*fn)(const T&)) const
    {
       matrix<T> a( mPtr->nrow, mPtr->ncol);
       size_t n = a.size();
    
       const T *pm = &mPtr->val[0];
       T *pa = &a.mPtr->val[0];
    
       for (size_t i=0; i < n; i++)
          pa[i] = fn( pm[i]);
    
       return a;
    }
    
    template <class T> bool 
    matrix<T>::eigen (valarray<T>& eival) const
    {
       #ifdef RANGE_CHECK_
       if (!isSymmetric())
          return false;
       #endif
    
       if (eival.size() != mPtr->nrow)
          eival.resize( mPtr->nrow);
    
       matrix<T> a(*this);
       valarray<T> e(mPtr->nrow);
    
       a.tred2( eival, e, false);
       return a.tql2( eival, e, false);
    }
    
    template <class T> bool 
    matrix<T>::eigen (valarray<T>& eival, matrix<T>& eivec) const
    {
       #ifdef RANGE_CHECK_
       if (!isSymmetric())
          return false;
       #endif
    
       valarray<T> e(mPtr->nrow);
    
       if (eival.size() != mPtr->nrow)
          eival.resize( mPtr->nrow);
    
       eivec = *this;
       eivec.tred2( eival, e, true);
    
       return eivec.tql2( eival, e, true);
    }
    
    template <class T> bool 
    matrix<T>::eigen (valarray<T>& rev, valarray<T>& iev) const
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          return false;
       #endif
    
       if (rev.size() != mPtr->nrow)
          rev.resize( mPtr->nrow);
       if (iev.size() != mPtr->nrow)
          iev.resize( mPtr->nrow);
    
       matrix<T> eivec, m(*this);
       m.balanc( eivec, false);
    
       return m.hqr2( rev, iev, eivec, false);
    }
    
    template <class T> bool 
    matrix<T>::eigen (valarray<T>& rev, valarray<T>& iev, matrix<T>& eivec) const
    {
       #ifdef RANGE_CHECK_
       if (mPtr->nrow != mPtr->ncol)
          return false;
       #endif
    
       if (rev.size() != mPtr->nrow)
          rev.resize( mPtr->nrow);
       if (iev.size() != mPtr->nrow)
          iev.resize( mPtr->nrow);
       if (eivec.rowno() != mPtr->nrow || eivec.colno() != mPtr->nrow)
          eivec.resize( mPtr->nrow, mPtr->nrow);
    
       matrix<T> m(*this);
       m.balanc( eivec, true);
    
       return m.hqr2( rev, iev, eivec, true);
    }
    
    template <class T> inline
    void cdiv (const T& xr, const T& xi, const T& yr, const T& yi, T& cdivr, T& cdivi) 
    {
       T r,d;
    
       if (abs( yr) > abs( yi)) 
       {
          r = yi/yr;
          d = yr + r*yi;
          cdivr = (xr + r*xi)/d;
          cdivi = (xi - r*xr)/d;
       } 
       else 
       {
          r = yr/yi;
          d = yi + r*yr;
          cdivr = (r*xr + xi)/d;
          cdivi = (r*xi - xr)/d;
       }
    }
    
    
    //  This is derived from the Algol procedure hqr2, by Martin 
    //  and Wilkinson, Handbook for Auto. Comp.,Vol.ii-Linear Algebra, 
    //  and the corresponding Fortran subroutine in EISPACK.
    //
    template <class T> bool
    matrix<T>::hqr2 (valarray<T>& d, valarray<T>& e, matrix<T>& v, bool eivec)
    {
       int i,j,k,l;
    
       int nn = int(mPtr->ncol);
       int n = nn-1;
       int low = 0;
       int high = nn-1;
       T exshift(0);
       T p(0),q(0),r(0),s(0),z(0),t,w,x,y;
       T cdivr, cdivi;
    
       if (shared())
          clone();
    
       T *pv, *pm = &mPtr->val[0];
       if (eivec)
          pv = &v(0,0);
    
       T norm(0);
       for (i=0; i < nn; i++) 
       {
          if (i < low || i > high) 
          {
             d[i] = pm[i*nn+i];
             e[i] = T(0);
          }
          for (j = techsoft::max(i-1,0); j < nn; j++)
             norm += abs( pm[i*nn+j]);
       }
    
       size_t iter = 0;
       while (n >= low) 
       {
          l = n;
          while (l > low) 
          {
             s = abs( pm[(l-1)*nn+l-1]) + abs( pm[l*nn+l]);
             if (s < epsilon(s))
                s = norm;
    
             if (abs( pm[l*nn+l-1]) < epsilon(pm[l*nn+l-1]) * s)
                break;
    
             l--;     
          }
    
          if (l == n) 
          {
             pm[n*nn+n] += exshift;
             d[n] = pm[n*nn+n];
             e[n] = T(0);
             n--;
             iter = 0;
          } 
          else if (l == n-1) 
          {
             w = pm[n*nn+n-1] * pm[(n-1)*nn+n];
             p = (pm[(n-1)*nn+n-1] - pm[n*nn+n]) / T(2);
             q = p * p + w;
             z = sqrt( abs( q));
             pm[n*nn+n] += exshift;
             pm[(n-1)*nn+n-1] += exshift;
             x = pm[n*nn+n];
    
             if (q >= T(0)) 
             {
                if (p >= T(0))
                   z = p + z;
                else
                   z = p - z;
    
                d[n-1] = x + z;
                d[n] = d[n-1];
                if (z != T(0))
                   d[n] = x - w / z;
    
                e[n-1] = T(0);
                e[n] = T(0);
    
                x = pm[n*nn+n-1];
                s = abs( x) + abs( z);
                p = x / s;
                q = z / s;
                r = sqrt( p*p + q*q);
                p = p / r;
                q = q / r;
    
                for (j = n-1; j < nn; j++) 
                {
                   z = pm[(n-1)*nn+j];
                   pm[(n-1)*nn+j] = q * z + p * pm[n*nn+j];
                   pm[n*nn+j] = q * pm[n*nn+j] - p * z;
                }
    
                for (i = 0; i <= n; i++) 
                {
                   z = pm[i*nn+n-1];
                   pm[i*nn+n-1] = q * z + p * pm[i*nn+n];
                   pm[i*nn+n] = q * pm[i*nn+n] - p * z;
                }
    
                if (eivec)
                   for (i = low; i <= high; i++) 
                   {
                      z = pv[i*nn+n-1];
                      pv[i*nn+n-1] = q * z + p * pv[i*nn+n];
                      pv[i*nn+n] = q * pv[i*nn+n] - p * z;
                   }
             } 
             else 
             {
                d[n-1] = x + p;
                d[n] = x + p;
                e[n-1] = z;
                e[n] = -z;
             }
             n = n - 2;
             iter = 0;
          } 
          else 
          {
             x = pm[n*nn+n];
             y = T(0);
             w = T(0);
             if (l < n) 
             {
                y = pm[(n-1)*nn+n-1];
                w = pm[n*nn+n-1] * pm[(n-1)*nn+n];
             }
    
             if (iter == 10) 
             {
                exshift += x;
                for (i = low; i <= n; i++)
                   pm[i*nn+i] -= x;
    
                s = abs( pm[n*nn+n-1]) + abs( pm[(n-1)*nn+n-2]);
                x = y = 0.75 * s;
                w = -0.4375 * s * s;
             }
    
             if (iter == 30) 
             {
                 s = (y - x) / T(2);
                 s = s * s + w;
                 if (s > T(0)) 
                 {
                     s = sqrt( s);
                     if (y < x)
                        s = -s;
    
                     s = x - w / ((y - x) / T(2) + s);
                     for (i = low; i <= n; i++)
                        pm[i*nn+i] -= s;
    
    
                     exshift += s;
                     x = y = w = 0.964;
                 }
             }
    
             if (++iter > 250)
                return false;
    
             int m = n-2;
             while (m >= l) 
             {
                z = pm[m*nn+m];
                r = x - z;
                s = y - z;
                p = (r * s - w) / pm[(m+1)*nn+m] + pm[m*nn+m+1];
                q = pm[(m+1)*nn+m+1] - z - r - s;
                r = pm[(m+2)*nn+m+1];
                s = abs( p) + abs( q) + abs( r);
                p = p / s;
                q = q / s;
                r = r / s;
                if (m == l)
                   break;
    
                if (abs( pm[m*nn+m-1]) * (abs( q) + abs( r)) <
                   epsilon(r) * (abs( p) * (abs( pm[(m-1)*nn+m-1]) + abs( z) +
                   abs( pm[(m+1)*nn+m+1])))) {
                      break;
                }
                m--;
             }
    
             for (i = m+2; i <= n; i++) 
             {
                pm[i*nn+i-2] = T(0);
                if (i > m+2)
                   pm[i*nn+i-3] = T(0);
             }
    
             for (k = m; k <= n-1; k++) 
             {
                bool notlast = (k != n-1);
                if (k != m) 
                {
                   p = pm[k*nn+k-1];
                   q = pm[(k+1)*nn+k-1];
                   r = notlast ? pm[(k+2)*nn+k-1] : T(0);
                   x = abs( p) + abs( q) + abs( r);
                   if (x > epsilon(x)) 
                   {
                      p = p / x;
                      q = q / x;
                      r = r / x;
                   }
                }
                if (x < epsilon(x))
                   break;
    
                s = sqrt( p * p + q * q + r * r);
                if (p < T(0))
                   s = -s;
    
                if (s != T(0)) 
                {
                   if (k != m)
                      pm[k*nn+k-1] = -s * x;
                   else if (l != m)
                      pm[k*nn+k-1] = -pm[k*nn+k-1];
    
                   p = p + s;
                   x = p / s;
                   y = q / s;
                   z = r / s;
                   q = q / p;
                   r = r / p;
    
                   for (j = k; j < nn; j++) 
                   {
                      p = pm[k*nn+j] + q * pm[(k+1)*nn+j];
                      if (notlast) 
                      {
                         p = p + r * pm[(k+2)*nn+j];
                         pm[(k+2)*nn+j] = pm[(k+2)*nn+j] - p * z;
                      }
                      pm[k*nn+j] = pm[k*nn+j] - p * x;
                      pm[(k+1)*nn+j] = pm[(k+1)*nn+j] - p * y;
                   }
    
                   for (i = 0; i <= techsoft::min( n, k+3); i++) 
                   {
                      p = x * pm[i*nn+k] + y * pm[i*nn+k+1];
                      if (notlast) 
                      {
                         p = p + z * pm[i*nn+k+2];
                         pm[i*nn+k+2] = pm[i*nn+k+2] - p * r;
                      }
                      pm[i*nn+k] = pm[i*nn+k] - p;
                      pm[i*nn+k+1] = pm[i*nn+k+1] - p * q;
                   }
    
                   if (eivec)
                      for (i = low; i <= high; i++) 
                      {
                         p = x * pv[i*nn+k] + y * pv[i*nn+k+1];
                         if (notlast) 
                         {
                            p += z * pv[i*nn+k+2];
                            pv[i*nn+k+2] -= p * r;
                         }
                         pv[i*nn+k] -= p;
                         pv[i*nn+k+1] -= p * q;
                      }
                }
             }
          }
       }
    
       if (norm < epsilon(norm))
          return true;
    
       for (n = nn-1; n >= 0; n--) 
       {
          p = d[n];
          q = e[n];
    
          if (q == T(0)) 
          {
             int l = n;
             pm[n*nn+n] = T(1);
             for (i = n-1; i >= 0; i--) 
             {
                w = pm[i*nn+i] - p;
                r = T(0);
                for (j = l; j <= n; j++)
                   r += pm[i*nn+j] * pm[j*nn+n];
    
                if (e[i] < T(0)) 
                {
                   z = w;
                   s = r;
                } 
                else 
                {
                   l = i;
                   if (e[i] == T(0)) 
                   {
                      if (w != T(0))
                         pm[i*nn+n] = -r / w;
                      else
                         pm[i*nn+n] = -r / (epsilon(norm) * norm);
                   } 
                   else 
                   {
                      x = pm[i*nn+i+1];
                      y = pm[(i+1)*nn+i];
                      q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
                      t = (x * s - z * r) / q;
                      pm[i*nn+n] = t;
                      if (abs( x) > abs( z))
                         pm[(i+1)*nn+n] = (-r - w * t) / x;
                      else
                         pm[(i+1)*nn+n] = (-s - y * t) / z;
                   }
    
                   t = abs( pm[i*nn+n]);
                   if ((epsilon(t) * t) * t > T(1)) 
                      for (j = i; j <= n; j++)
                         pm[j*nn+n] /= t;
                }
             }
          } 
          else if (q < T(0)) 
          {
             int l = n-1;
    
             if (abs( pm[n*nn+n-1]) > abs( pm[(n-1)*nn+n])) 
             {
                pm[(n-1)*nn+n-1] = q / pm[n*nn+n-1];
                pm[(n-1)*nn+n] = -(pm[n*nn+n] - p) / pm[n*nn+n-1];
             } 
             else 
             {
                cdiv( 0.0, -pm[(n-1)*nn+n], pm[(n-1)*nn+n-1]-p, q, cdivr, cdivi);
                pm[(n-1)*nn+n-1] = cdivr;
                pm[(n-1)*nn+n] = cdivi;
             }
             pm[n*nn+n-1] = T(0);
             pm[n*nn+n] = T(1);
             for (i = n-2; i >= 0; i--) 
             {
                T ra(0),sa(0),vr,vi;
    
                for (j = l; j <= n; j++) 
                {
                   ra += pm[i*nn+j] * pm[j*nn+n-1];
                   sa += pm[i*nn+j] * pm[j*nn+n];
                }
                w = pm[i*nn+i] - p;
    
                if (e[i] < T(0)) 
                {
                   z = w;
                   r = ra;
                   s = sa;
                } 
                else 
                {
                   l = i;
                   if (e[i] == T(0))
                   {
                      cdiv( -ra, -sa, w, q, pm[i*nn+n-1], pm[i*nn+n]);
                   } 
                   else 
                   {
                      x = pm[i*nn+i+1];
                      y = pm[(i+1)*nn+i];
                      vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
                      vi = (d[i] - p) * 2.0 * q;
                      if (vr == T(0) && vi == T(0)) 
                         vr = epsilon(norm) * norm * (abs( w) + abs(q) + abs(x) + abs(y) + abs(z));
    
                      cdiv( x*r - z*ra + q*sa, x*s - z*sa - q*ra, vr, vi, cdivr, cdivi);
                      pm[i*nn+n-1] = cdivr;
                      pm[i*nn+n] = cdivi;
                      if (abs(x) > (abs(z) + abs(q))) 
                      {
                         pm[(i+1)*nn+n-1] = (-ra - w * pm[i*nn+n-1] + q * pm[i*nn+n]) / x;
                         pm[(i+1)*nn+n] = (-sa - w * pm[i*nn+n] - q * pm[i*nn+n-1]) / x;
                      } 
                      else 
                      {
                         cdiv( -r - y * pm[i*nn+n-1], -s - y * pm[i*nn+n], z, q, cdivr, cdivi);
                         pm[(i+1)*nn+n-1] = cdivr;
                         pm[(i+1)*nn+n] = cdivi;
                      }
                   }
    
                   t = techsoft::max( abs( pm[i*nn+n-1]), abs( pm[i*nn+n]));
                   if ((epsilon(t) * t) * t > T(1)) 
                      for (j = i; j <= n; j++) 
                      {
                         pm[j*nn+n-1] /= t;
                         pm[j*nn+n] /= t;
                      }
                }
             }
          }
       }
    
       if (eivec)
       {
          for (i = 0; i < nn; i++) 
             if (i < low || i > high) 
                for (j = i; j < nn; j++)
                   pv[i*nn+j] = pm[i*nn+j];
    
          for (j = nn-1; j >= low; j--) 
             for (i = low; i <= high; i++) 
             {
                z = T(0);
                for (k = low; k <= techsoft::min( j, high); k++)
                   z += pv[i*nn+k] * pm[k*nn+j];
    
                pv[i*nn+j] = z;
             }
       }
       return true;
    }
    
    ///
    //  This is derived from the Algol procedures orthes and ortran,
    //  by Martin and Wilkinson, Handbook for Auto. Comp.,Vol.ii-Line Algebra, and the corresponding Fortran subroutines in EISPACK.
    //
    template <class T> void 
    matrix<T>::balanc (matrix<T>& v, bool eivec) 
    {
       size_t i,j,lo,hi,m,n;
    
       if (shared())
          clone();
    
       n = mPtr->ncol;
       lo = 0;
       hi = n-1;
       valarray<T> ort(n);
       T *pv, *pm = &mPtr->val[0];
       if (eivec)
          pv = &v(0,0);
    
       for (m=lo+1; m <= hi-1; m++) 
       {
          T scale(0);
          for (i=m; i <= hi; i++)
             scale += abs( pm[i*n+m-1]);
    
          if (scale > epsilon(scale)) 
          {
             T h(0);
             for (i=hi; i >= m; i--) 
             {
                ort[i] = pm[i*n+m-1] / scale;
                h += ort[i] * ort[i];
             }
             T g = sqrt( h);
             if (ort[m] > T(0)) 
                g = -g;
    
             h -= ort[m] * g;
             ort[m] -= g;
    
             for (j=m; j < n; j++) 
             {
                T f(0);
                for (i = hi; i >= m; i--) 
                   f += ort[i] * pm[i*n+j];
                f /= h;
                for (i=m; i <= hi; i++)
                   pm[i*n+j] -= f * ort[i];
            }
    
            for (i=0; i <= hi; i++) 
            {
                T f(0);
                for (j=hi; j >= m; j--)
                   f += ort[j] * pm[i*n+j];
                f /= h;
                for (j=m; j <= hi; j++)
                   pm[i*n+j] -= f * ort[j];
             }
             ort[m] = scale * ort[m];
             pm[m*n+m-1] = scale * g;
          }
       }
    
       if (eivec)
          for (i=0; i < n; i++)
             for (j = 0; j < n; j++)
                pv[i*n+j] = (i == j ? T(1) : T(0));
    
       for (m=hi-1; m >= lo+1; m--) 
       {
          if (abs(pm[m*n+m-1]) > T(0)) 
          {
             for (i=m+1; i <= hi; i++)
                ort[i] = pm[i*n+m-1];
    
             if (eivec)
                for (j=m; j <= hi; j++) 
                {
                   T g(0);
                   for (i=m; i <= hi; i++)
                      g += ort[i] * pv[i*n+j];
    
                   g = (g / ort[m]) / pm[m*n+m-1];
                   for (i=m; i <= hi; i++)
                      pv[i*n+j] += g * ort[i];
                }
          }
       }
    }
    
    template <class T> bool
    matrix<T>::tql2 (valarray<T>& d, valarray<T>& e, bool eivec)
    {
       size_t m,l,iter,i,k;
       T s,r,p,g,f,c,b;
    
       size_t n = mPtr->nrow;
       T *pm = &mPtr->val[0];
    
       for (i=1; i < n; i++)
          e[i-1] = e[i]; 
       e[n-1] = T(0);
    
       for (l=0; l < n; l++) 
       {
          iter=0;
          do 
          {
             for (m=l; m < n-1; m++)
             { 
                T dd = abs( d[m]) + abs( d[m+1]);
                if (abs( e[m]) <= epsilon(e[m])*dd) 
                   break;
             }
             if (m != l) 
             {
                if (iter++ == 30) 
                   return false;
    
                g = (d[l+1] - d[l]) / (T(2) * e[l]);
                r = hypot( g, T(1));
                g = d[m] - d[l] + e[l] / (g + sign( r, g));
                s = c = T(1);
                p = T(0);
                for (i=m-1; ; i--)
                {
                   f = s * e[i];
                   b = c * e[i];
                   r = hypot( f, g);
                   e[i+1] = r;
                   if (r < epsilon(r))
                   { 
                      d[i+1] -= p;
                      e[m] = T(0);
                      break;
                   }
                   s = f / r;
                   c = g / r;
                   g = d[i+1] - p;
                   r = (d[i] - g) * s + T(2) * c * b;
                   p = s * r;
                   d[i+1] = g + p;
                   g = c * r - b;
    
                   if (eivec)
                      for (k=0; k < n; k++)
                      { 
                         f = pm[k*n+i+1];
                         pm[k*n+i+1] = s * pm[k*n+i] + c * f;
                         pm[k*n+i] = c * pm[k*n+i] - s * f;
                      }
                   if (i == l)
                      break;
                }
                if (abs(r) < epsilon(r) && i >= l) 
                   continue;
                d[l] -= p;
                e[l] = g;
                e[m] = T(0);
             }
          } 
          while (m != l);
    
          if (abs(d[l]) < epsilon(d[l]))
             d[l] = T(0);
       }
       return true;
    }
    
    template <class T> void
    matrix<T>::tred2 (valarray<T>& d, valarray<T>& e, bool eivec)
    {
       size_t i,j,k;
       size_t n = mPtr->nrow;
    
       #ifdef RANGE_CHECK_
       if (n != mPtr->ncol)
          throw std::invalid_argument( "matrix<T>::tred2: non-square matrix!");
       #endif
    
       if (shared())
          clone();
       T *pm = &mPtr->val[0];
    
       for (i=n-1; i > 0; i--) 
       {
          T h(0), scale(0);
    
          if (i > 1) 
          {
             for (k=0; k < i; k++)
                scale += abs( pm[i*n+k]);
    
             if (scale < epsilon(scale))
                e[i] = pm[i*n+i-1];
             else 
             {
                for (k=0; k < i; k++) 
                {
                   pm[i*n+k] /= scale;
                   T tmp = pm[i*n+k]; 
                   h += tmp * tmp;
                }
                T f = pm[i*n+i-1];
                T g = sqrt( h);
                if (f >= T(0))
                   g = -g;
    
                e[i] = scale * g;
                h -= f * g;
                pm[i*n+i-1] = f - g;
                f = T(0);
                for (j=0; j < i; j++) 
                {
                   if (eivec)
                      pm[j*n+i] = pm[i*n+j] / h;
                   g = T(0);
                   for (k=0; k <= j; k++)
                      g += pm[j*n+k] * pm[i*n+k];
                   for (k=j+1; k < i; k++)
                      g += pm[k*n+j] * pm[i*n+k];
                   e[j] = g / h;
                   f += e[j] * pm[i*n+j];
                }
                T hh = f / (h+h);
                for (j=0; j < i; j++)
                { 
                   f = pm[i*n+j];
                   g = e[j] - hh * f;
                   e[j] = g;
                   for (k=0; k <= j; k++)
                      pm[j*n+k] -= (f * e[k] + g * pm[i*n+k]);
                }
             }
          } 
          else
             e[i] = pm[i*n+i-1];
          d[i] = h;
       }
    
       if (eivec)
          d[0] = T(0);
       e[0] = T(0);
    
       for (i=0; i < n; i++)
       { 
          if (eivec)
          {
             if (d[i] != T(0))
             { 
                for (j=0; j < i; j++) 
                {
                   T g(0);
                   for (k=0; k < i; k++)
                      g += pm[i*n+k] * pm[k*n+j];
                   for (k=0; k < i; k++)
                      pm[k*n+j] -= g * pm[k*n+i];
                }
             }
             d[i] = pm[i*n+i];
             pm[i*n+i] = T(1);
             for (j=0; j < i; j++) 
                pm[j*n+i] = pm[i*n+j] = T(0);
          }
          else
             d[i] = pm[i*n+i];
       }
    }
    
    
    /
    //  I/O functions
    //
    
    #ifndef MATRIX_ROW_SEP
    #define MATRIX_ROW_SEP   ('\n')
    #endif
    
    #ifndef MATRIX_COL_SEP
    #define MATRIX_COL_SEP   ('\t')
    #endif
    
    // Input stream function for a single element
    template <class T> inline istream&
    operator>> (istream& is, Mref<T> el)
    {
        T x; is >> x; el = x; return is;
    }
    
    // Output stream function for a single element
    template <class T> inline ostream&
    operator<< (ostream& os, const Mref<T>& el)
    {
        const T x = el;
        os << x;
        return os;
    }
    
    // Output stream function for matrix row/column
    template <class T> ostream&
    operator<< (ostream& os, const Mat_iter<T>& v)
    {
        int w = os.width();
        for (size_t i=0; i < v.size(); i++)
            os << setw( w) << v[i] << MATRIX_COL_SEP;
        return os;
    }
    
    // Iutput stream function for matrix row/column
    template <class T> istream&
    operator>> (istream& is, Mat_iter<T> v)
    {
        for (size_t i=0; i < v.size(); i++)
            is >> v[i];
        return is;
    }
    
    
    // Output stream function for const matrix row/column
    template <class T> ostream&
    operator<< (ostream& os, const Cmat_iter<T>& v)
    {
        int w = os.width();
        for (size_t i=0; i < v.size(); i++)
            os << setw( w) << v[i] << MATRIX_COL_SEP;
        return os;
    }
    
    // Input stream function
    template <class T> istream&
    operator>> (istream& is, matrix<T>& m)
    {
        T *pm = &m(0,0);
        size_t n = m.size();
        for (size_t i=0; i < n; i++)
           is >> pm[i];
        return is;
    }
    
    // Output stream function for matrix
    template <class T> ostream&
    operator<< (ostream &os, const matrix<T>& m)
    {
        int w = os.width();
        for (size_t i=0; i < m.rowno(); i++)
           for (size_t j=0; j < m.colno(); j++)
              os << setw( w) << m(i,j) << (j == m.colno()-1 ? MATRIX_ROW_SEP : MATRIX_COL_SEP);
       return os;
    }
    
    
    // Iutput stream function for valarray
    template <class T> istream&
    operator>> (istream& is, valarray<T>& v)
    {
        for (size_t i=0; i < v.size(); i++)
            is >> v[i];
       return is;
    }
    
    // Output stream function for valarray
    template <class T> ostream&
    operator<< (ostream& os, const valarray<T>& v)
    {
        int w = os.width();
        for (size_t i=0; i < v.size(); i++)
            os << setw( w) << v[i] << MATRIX_COL_SEP;
        return os;
    }
    
    
    //
    // Specializations
    //
    
    template <> inline 
    void matrix<complex<float> >::rand (int rmin, int rmax, int rseed)
    {
       typedef value_type T;
       typedef float VT;
    
       if (!rmin && !rmax)
       {
          rmin = -1;
          rmax = 1;
       }
       if (rmin > rmax)
       {
          int tmp = rmin;
          rmin = rmax;
          rmax = tmp;
       }
       int rdiff = rmax - rmin;
    
       if (rseed)
          srand( rseed);
    
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       }
    
       T *pv = &mPtr->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; ++i)
       {
          VT rd[2];
    
          for (size_t j=0; j < 2; ++j)
          {
             int rn = rand_();
             if (rdiff > 1)
             {
                rn %= rdiff;
                rd[j] = VT(rn) - VT(rdiff) / VT(2);
                if (rdiff > 1)
                {
                   VT f = ((VT)rand_()) / VT(RAND_MAX);
                   if (rn > rdiff)
                      rd[j] -= f;
                   else
                      rd[j] += f;
                }
             }
             else
                rd[j] = VT(rn) / VT(RAND_MAX);
          }
          pv[i] = T( rd[0], rd[1]);
       }
    }
    
    template <> inline 
    void matrix<complex<double> >::rand (int rmin, int rmax, int rseed)
    {
       typedef value_type T;
       typedef double VT;
    
       value_type y = value_type(1.0);
    
       if (!rmin && !rmax)
       {
          rmin = -1;
          rmax = 1;
       }
       if (rmin > rmax)
       {
          int tmp = rmin;
          rmin = rmax;
          rmax = tmp;
       }
       int rdiff = rmax - rmin;
    
       if (rseed)
          srand( rseed);
    
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       }
    
       T *pv = &mPtr->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; ++i)
       {
          VT rd[2];
    
          for (size_t j=0; j < 2; ++j)
          {
             int rn = rand_();
             if (rdiff > 1)
             {
                rn %= rdiff;
                rd[j] = VT(rn) - VT(rdiff) / VT(2);
                if (rdiff > 1)
                {
                   VT f = ((VT)rand_()) / VT(RAND_MAX);
                   if (rn > rdiff)
                      rd[j] -= f;
                   else
                      rd[j] += f;
                }
             }
             else
                rd[j] = VT(rn) / VT(RAND_MAX);
          }
          pv[i] = T( rd[0], rd[1]);
       }
    }
    
    template <> inline 
    void matrix<complex<long double> >::rand (int rmin, int rmax, int rseed)
    {
       typedef value_type T;
       typedef long double VT;
    
       if (!rmin && !rmax)
       {
          rmin = -1;
          rmax = 1;
       }
       if (rmin > rmax)
       {
          int tmp = rmin;
          rmin = rmax;
          rmax = tmp;
       }
       int rdiff = rmax - rmin;
    
       if (rseed)
          srand( rseed);
    
       if (shared())
       {
          --mPtr->refcnt;
          mPtr = allocator( Allocate, mPtr->nrow, mPtr->ncol);
       }
    
       T *pv = &mPtr->val[0];
       size_t n = size();
    
       for (size_t i=0; i < n; ++i)
       {
          VT rd[2];
    
          for (size_t j=0; j < 2; ++j)
          {
             int rn = rand_();
             if (rdiff > 1)
             {
                rn %= rdiff;
                rd[j] = VT(rn) - VT(rdiff) / VT(2);
                if (rdiff > 1)
                {
                   VT f = ((VT)rand_()) / VT(RAND_MAX);
                   if (rn > rdiff)
                      rd[j] -= f;
                   else
                      rd[j] += f;
                }
             }
             else
                rd[j] = VT(rn) / VT(RAND_MAX);
          }
          pv[i] = T( rd[0], rd[1]);
       }
    }
    
    template <> inline
    complex<float> sign (complex<float> a, complex<float> b)
    {
       if (abs( a-b) < epsilon( a))
          return a;
    
       return (b.real() < 0.0f || b.imag() < 0.0f)? -a : a;
    }
    
    template <> inline
    complex<double> sign (complex<double> a, complex<double> b)
    {
       if (abs( a-b) < epsilon( a))
          return a;
    
       return (b.real() < double(0) || b.imag() < double(0))? -a : a;
    }
    
    template <> inline
    complex<long double> sign (complex<long double> a, complex<long double> b)
    {
       if (abs( a-b) < epsilon( a))
          return a;
    
       return (b.real() < 0.0 || b.imag() < 0.0)? -a : a;
    }
    
    template <> inline 
    bool matrix<complex<float> >::svd (matrix<complex<float> >& vc, valarray<complex<float> >& w)
    {
       typedef value_type T;
       typedef float VT;
    
       size_t flag,i,its,j,jj,k,l,nm;
       T c,f,g(0),h,s,x,y,z,tmp;
    
       if (shared())
          clone();
    
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
    
       if (vc.rowno() != n || vc.colno() != n)
          vc.resize( n, n);
       if (w.size() != n)
          w.resize( n);
    
       T *a = &mPtr->val[0];
       T *v = &vc(0,0);
       valarray<T> rv1(n);
    
       VT scale(0), anorm(0);
       for (i=0; i < n; i++)
       {
          l = i + 1;
          rv1[i] = scale * g;
          scale = VT( 0);
          s = T( 0);
          g = T( 0);
    
          if (i < m)
          {
             for (k=i; k < m; k++)
                scale += abs( a[k*n+i]);
    
             if (scale > epsilon( scale))
             {
                for (k=i; k < m; k++)
                {
                   a[k*n+i] /= scale;
                   tmp = a[k*n+i];
                   s += tmp * tmp;
                }
                f = a[i*n+i];
                g = -sign( sqrt(s), f);
                h = f * g - s;
                a[i*n+i] = f - g;
    
                for (j=l; j < n; j++)
                {
                   for (s=T(0),k=i; k < m; k++)
                      s += a[k*n+i] * a[k*n+j];
                   f = s / h;
                   for (k=i; k < m; k++)
                      a[k*n+j] += f * a[k*n+i];
                }
                for (k=i; k < m; k++)
                   a[k*n+i] *= scale;
             }
          }
    
          w[i] = scale * g;
          scale = VT( 0);
          s = T(0);
          g = T(0);
    
          if (i < m && i != n-1)
          {
             for (k=l; k < n; k++)
                scale += abs(a[i*n+k]);
    
             if (scale > epsilon(scale))
             {
                for (k=l; k < n; k++)
                {
                   a[i*n+k] /= scale;
                   tmp = a[i*n+k];
                   s += tmp * tmp;
                }
                f = a[i*n+l];
                g = -sign( sqrt(s), f);
                h = f * g - s;
                a[i*n+l] = f - g;
    
                for (k=l; k < n; k++)
                   rv1[k] = a[i*n+k] / h;
    
                for (j=l; j < m; j++)
                {
                   for (s=T(0),k=l; k < n; k++)
                      s += a[j*n+k] * a[i*n+k];
                   for (k=l; k < n; k++)
                      a[j*n+k] += s * rv1[k];
                }
                for (k=l; k < n; k++)
                   a[i*n+k] *= scale;
             }
          }
          anorm = techsoft::max( anorm, abs( w[i]) + abs( rv1[i]) );
       }
    
       for (i=n-1; ; i--)
       {
          if (i < n-1)
          {
             if (abs(g) > epsilon(g))
             {
                for (j=l; j < n; j++)
                    v[j*n+i] = (a[i*n+j] / a[i*n+l]) / g;
                for (j=l; j < n; j++)
                {
                   for (s=T(0),k=l; k < n; k++)
                      s += a[i*n+k] * v[k*n+j];
    
                   for (k=l; k < n; k++)
                      v[k*n+j] += s * v[k*n+i];
                }
             }
             for (j=l; j < n; j++)
                v[i*n+j] = v[j*n+i] = T(0);
          }
    
          v[i*n+i] = T(1);
          g = rv1[i];
          l = i;
          if (i == 0)
             break;
       }
    
       for (i=techsoft::min(m,n)-1; ; i--)
       {
          l = i + 1;
          g = w[i];
          for (j=l; j < n; j++)
             a[i*n+j] = T(0);
          if (abs(g) > epsilon(g))
          {
             g = T(1) / g;
             for (j=l; j < n; j++)
             {
                for (s=T(0),k=l; k < m; k++)
                   s += a[k*n+i] * a[k*n+j];
    
                f = (s / a[i*n+i]) * g;
    
                for (k=i; k < m; k++)
                   a[k*n+j] += f * a[k*n+i];
             }
             for (j=i; j < m; j++)
                a[j*n+i] *= g;
    
          }
          else
          {
              for (j=i; j < m; j++)
                  a[j*n+i] = T(0);
          }
          a[i*n+i] = a[i*n+i] + T(1);
          if (i == 0)
             break;
       }
    
       for (k=n-1; ; k--)
       {
          for (its=1; its <= 30; its++)
          {
             flag = 1;
             for (l=k; ; l--)
             {
                nm = l-1;
                if (abs( rv1[l]) < epsilon(rv1[l]))
                {
                   flag = 0;
                   break;
                }
                if (abs( w[nm]) < epsilon(w[nm]))
                   break;
                if (l == 0)
                   break;
             }
             if (flag)
             {
                c = T(0);
                s = T(1);
                for (i=l; i <= k; i++)
                {
                   f = s * rv1[i];
                   rv1[i] = c * rv1[i];
                   if (abs(f) < epsilon(f))
                      break;
                   g = w[i];
                   h = sqrt( f * f + g * g);
                   w[i] = h;
                   h = VT(1) / h;
                   c = g * h;
                   s = -f * h;
                   for (j=0; j < m; j++)
                   {
                      y = a[j*n+nm];
                      z = a[j*n+i];
                      a[j*n+nm] = y * c + z * s;
                      a[j*n+i] = z * c - y * s;
                   }
                }
             }
             z = w[k];
             if (l == k)                 
             {
                if (z.real() < VT(0))
                {
                   w[k] = -z;
                   for (j=0; j < n; j++)
                      v[j*n+k] = -v[j*n+k];
                }
                break;
             }
    
             if (its == 30)
                return false;
    
             x = w[l];
             nm = k - 1;
             y = w[nm];
             g = rv1[nm];
             h = rv1[k];
             f = ((y-z) * (y+z) + (g-h) * (g+h)) / (VT(2) * h * y);
             g = sqrt( f * f + T(1) * T(1));
             f = ((x-z) * (x+z) + h * ((y / (f+sign( g, f))) - h)) / x;
             c = s = T(1);                    
             for (j=l; j <= nm; j++)
             {
                i = j + 1;
                g = rv1[i];
                y = w[i];
                h = s * g;
                g = c * g;
                z = sqrt( f * f + h * h);
                rv1[j] = z;
                c = f / z;
                s = h / z;
                f = x * c + g * s;
                g = g*c - x*s;
                h = y * s;
                y *= c;
                for (jj=0; jj < n; jj++)
                {
                   x = v[jj*n+j];
                   z = v[jj*n+i];
                   v[jj*n+j] = x * c + z * s;
                   v[jj*n+i] = z * c - x * s;
                }
                z = sqrt( f * f + h * h);
                w[j] = z;
                if (abs(z) > epsilon(z))
                {
                   z = VT(1) / z;
                   c = f * z;
                   s = h * z;
                }
                f = c * g + s * y;
                x = c * y - s * g;
                for (jj=0; jj < m; jj++)
                {
                   y = a[jj*n+j];
                   z = a[jj*n+i];
                   a[jj*n+j] = y * c + z * s;
                   a[jj*n+i] = z * c - y * s;
                }
             }
             rv1[l] = 0.0;
             rv1[k] = f;
             w[k] = x;
          }
          if (k == 0)
             break;
       }
       return true;
    }
    
    template <> inline 
    bool matrix<complex<double> >::svd (matrix<complex<double> >& vc, valarray<complex<double> >& w)
    {
       typedef value_type T;
       typedef double VT;
    
       size_t flag,i,its,j,jj,k,l,nm;
       T c,f,g(0),h,s,x,y,z,tmp;
    
       if (shared())
          clone();
    
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
    
       if (vc.rowno() != n || vc.colno() != n)
          vc.resize( n, n);
       if (w.size() != n)
          w.resize( n);
    
       T *a = &mPtr->val[0];
       T *v = &vc(0,0);
       valarray<T> rv1(n);
    
       VT scale(0), anorm(0);
       for (i=0; i < n; i++)
       {
          l = i + 1;
          rv1[i] = scale * g;
          scale = VT( 0);
          s = T( 0);
          g = T( 0);
    
          if (i < m)
          {
             for (k=i; k < m; k++)
                scale += abs( a[k*n+i]);
    
             if (scale > epsilon( scale))
             {
                for (k=i; k < m; k++)
                {
                   a[k*n+i] /= scale;
                   tmp = a[k*n+i];
                   s += tmp * tmp;
                }
                f = a[i*n+i];
                g = -sign( sqrt(s), f);
                h = f * g - s;
                a[i*n+i] = f - g;
    
                for (j=l; j < n; j++)
                {
                   for (s=T(0),k=i; k < m; k++)
                      s += a[k*n+i] * a[k*n+j];
                   f = s / h;
                   for (k=i; k < m; k++)
                      a[k*n+j] += f * a[k*n+i];
                }
                for (k=i; k < m; k++)
                   a[k*n+i] *= scale;
             }
          }
    
          w[i] = scale * g;
          scale = VT( 0);
          s = T(0);
          g = T(0);
    
          if (i < m && i != n-1)
          {
             for (k=l; k < n; k++)
                scale += abs(a[i*n+k]);
    
             if (scale > epsilon(scale))
             {
                for (k=l; k < n; k++)
                {
                   a[i*n+k] /= scale;
                   tmp = a[i*n+k];
                   s += tmp * tmp;
                }
                f = a[i*n+l];
                g = -sign( sqrt(s), f);
                h = f * g - s;
                a[i*n+l] = f - g;
    
                for (k=l; k < n; k++)
                   rv1[k] = a[i*n+k] / h;
    
                for (j=l; j < m; j++)
                {
                   for (s=T(0),k=l; k < n; k++)
                      s += a[j*n+k] * a[i*n+k];
                   for (k=l; k < n; k++)
                      a[j*n+k] += s * rv1[k];
                }
                for (k=l; k < n; k++)
                   a[i*n+k] *= scale;
             }
          }
          anorm = techsoft::max( anorm, abs( w[i]) + abs( rv1[i]) );
       }
    
       for (i=n-1; ; i--)
       {
          if (i < n-1)
          {
             if (abs(g) > epsilon(g))
             {
                for (j=l; j < n; j++)
                    v[j*n+i] = (a[i*n+j] / a[i*n+l]) / g;
                for (j=l; j < n; j++)
                {
                   for (s=T(0),k=l; k < n; k++)
                      s += a[i*n+k] * v[k*n+j];
    
                   for (k=l; k < n; k++)
                      v[k*n+j] += s * v[k*n+i];
                }
             }
             for (j=l; j < n; j++)
                v[i*n+j] = v[j*n+i] = T(0);
          }
    
          v[i*n+i] = T(1);
          g = rv1[i];
          l = i;
          if (i == 0)
             break;
       }
    
       for (i=techsoft::min(m,n)-1; ; i--)
       {
          l = i + 1;
          g = w[i];
          for (j=l; j < n; j++)
             a[i*n+j] = T(0);
          if (abs(g) > epsilon(g))
          {
             g = T(1) / g;
             for (j=l; j < n; j++)
             {
                for (s=T(0),k=l; k < m; k++)
                   s += a[k*n+i] * a[k*n+j];
    
                f = (s / a[i*n+i]) * g;
    
                for (k=i; k < m; k++)
                   a[k*n+j] += f * a[k*n+i];
             }
             for (j=i; j < m; j++)
                a[j*n+i] *= g;
    
          }
          else
          {
              for (j=i; j < m; j++)
                  a[j*n+i] = T(0);
          }
          a[i*n+i] = a[i*n+i] + T(1);
          if (i == 0)
             break;
       }
    
       for (k=n-1; ; k--)
       {
          for (its=1; its <= 30; its++)
          {
             flag = 1;
             for (l=k; ; l--)
             {
                nm = l-1;
                if (abs( rv1[l]) < epsilon(rv1[l]))
                {
    
                   flag = 0;
                   break;
                }
                if (abs( w[nm]) < epsilon(w[nm]))
                   break;
                if (l == 0)
                   break;
             }
             if (flag)
             {
                c = T(0);
                s = T(1);
                for (i=l; i <= k; i++)
                {
                   f = s * rv1[i];
                   rv1[i] = c * rv1[i];
                   if (abs(f) < epsilon(f))
                      break;
                   g = w[i];
                   h = sqrt( f * f + g * g);
                   w[i] = h;
                   h = VT(1) / h;
                   c = g * h;
                   s = -f * h;
                   for (j=0; j < m; j++)
                   {
                      y = a[j*n+nm];
                      z = a[j*n+i];
                      a[j*n+nm] = y * c + z * s;
                      a[j*n+i] = z * c - y * s;
                   }
                }
             }
             z = w[k];
             if (l == k)                 
             {
                if (z.real() < VT(0))
                {
                   w[k] = -z;
                   for (j=0; j < n; j++)
                      v[j*n+k] = -v[j*n+k];
                }
                break;
             }
    
             if (its == 30)
                return false;
    
             x = w[l];
             nm = k - 1;
             y = w[nm];
             g = rv1[nm];
             h = rv1[k];
             f = ((y-z) * (y+z) + (g-h) * (g+h)) / (VT(2) * h * y);
             g = sqrt( f * f + T(1) * T(1));
             f = ((x-z) * (x+z) + h * ((y / (f+sign( g, f))) - h)) / x;
             c = s = T(1);                    
             for (j=l; j <= nm; j++)
             {
                i = j + 1;
                g = rv1[i];
                y = w[i];
                h = s * g;
                g = c * g;
                z = sqrt( f * f + h * h);
                rv1[j] = z;
                c = f / z;
                s = h / z;
                f = x * c + g * s;
                g = g*c - x*s;
                h = y * s;
                y *= c;
                for (jj=0; jj < n; jj++)
                {
                   x = v[jj*n+j];
                   z = v[jj*n+i];
                   v[jj*n+j] = x * c + z * s;
                   v[jj*n+i] = z * c - x * s;
                }
                z = sqrt( f * f + h * h);
                w[j] = z;
                if (abs(z) > epsilon(z))
                {
                   z = VT(1) / z;
                   c = f * z;
                   s = h * z;
                }
                f = c * g + s * y;
                x = c * y - s * g;
                for (jj=0; jj < m; jj++)
                {
                   y = a[jj*n+j];
                   z = a[jj*n+i];
                   a[jj*n+j] = y * c + z * s;
                   a[jj*n+i] = z * c - y * s;
                }
             }
             rv1[l] = 0.0;
             rv1[k] = f;
             w[k] = x;
          }
          if (k == 0)
             break;
       }
       return true;
    }
    
    template <> inline 
    bool matrix<complex<long double> >::svd (matrix<complex<long double> >& vc, valarray<complex<long double> >& w)
    {
       typedef value_type T;
       typedef long double VT;
    
       size_t flag,i,its,j,jj,k,l,nm;
       T c,f,g(0),h,s,x,y,z,tmp;
    
       if (shared())
          clone();
    
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
    
       if (vc.rowno() != n || vc.colno() != n)
          vc.resize( n, n);
       if (w.size() != n)
          w.resize( n);
    
       T *a = &mPtr->val[0];
       T *v = &vc(0,0);
       valarray<T> rv1(n);
    
       VT scale(0), anorm(0);
       for (i=0; i < n; i++)
       {
          l = i + 1;
          rv1[i] = scale * g;
          scale = VT( 0);
          s = T( 0);
          g = T( 0);
    
          if (i < m)
          {
             for (k=i; k < m; k++)
                scale += abs( a[k*n+i]);
    
             if (scale > epsilon( scale))
             {
                for (k=i; k < m; k++)
                {
                   a[k*n+i] /= scale;
                   tmp = a[k*n+i];
                   s += tmp * tmp;
                }
                f = a[i*n+i];
                g = -sign( sqrt(s), f);
                h = f * g - s;
                a[i*n+i] = f - g;
    
                for (j=l; j < n; j++)
                {
                   for (s=T(0),k=i; k < m; k++)
                      s += a[k*n+i] * a[k*n+j];
                   f = s / h;
                   for (k=i; k < m; k++)
                      a[k*n+j] += f * a[k*n+i];
                }
                for (k=i; k < m; k++)
                   a[k*n+i] *= scale;
             }
          }
    
          w[i] = scale * g;
          scale = VT( 0);
          s = T(0);
          g = T(0);
    
          if (i < m && i != n-1)
          {
             for (k=l; k < n; k++)
                scale += abs(a[i*n+k]);
    
             if (scale > epsilon(scale))
             {
                for (k=l; k < n; k++)
                {
                   a[i*n+k] /= scale;
                   tmp = a[i*n+k];
                   s += tmp * tmp;
                }
                f = a[i*n+l];
                g = -sign( sqrt(s), f);
                h = f * g - s;
                a[i*n+l] = f - g;
    
                for (k=l; k < n; k++)
                   rv1[k] = a[i*n+k] / h;
    
                for (j=l; j < m; j++)
                {
                   for (s=T(0),k=l; k < n; k++)
                      s += a[j*n+k] * a[i*n+k];
                   for (k=l; k < n; k++)
                      a[j*n+k] += s * rv1[k];
                }
                for (k=l; k < n; k++)
                   a[i*n+k] *= scale;
             }
          }
          anorm = techsoft::max( anorm, abs( w[i]) + abs( rv1[i]) );
       }
    
       for (i=n-1; ; i--)
       {
          if (i < n-1)
          {
             if (abs(g) > epsilon(g))
             {
                for (j=l; j < n; j++)
                    v[j*n+i] = (a[i*n+j] / a[i*n+l]) / g;
                for (j=l; j < n; j++)
                {
                   for (s=T(0),k=l; k < n; k++)
                      s += a[i*n+k] * v[k*n+j];
    
                   for (k=l; k < n; k++)
                      v[k*n+j] += s * v[k*n+i];
                }
             }
             for (j=l; j < n; j++)
                v[i*n+j] = v[j*n+i] = T(0);
          }
    
          v[i*n+i] = T(1);
          g = rv1[i];
          l = i;
          if (i == 0)
             break;
       }
    
       for (i=techsoft::min(m,n)-1; ; i--)
       {
          l = i + 1;
          g = w[i];
          for (j=l; j < n; j++)
             a[i*n+j] = T(0);
          if (abs(g) > epsilon(g))
          {
             g = T(1) / g;
             for (j=l; j < n; j++)
             {
                for (s=T(0),k=l; k < m; k++)
                   s += a[k*n+i] * a[k*n+j];
    
                f = (s / a[i*n+i]) * g;
    
                for (k=i; k < m; k++)
                   a[k*n+j] += f * a[k*n+i];
             }
             for (j=i; j < m; j++)
                a[j*n+i] *= g;
    
          }
          else
          {
              for (j=i; j < m; j++)
                  a[j*n+i] = T(0);
          }
          a[i*n+i] = a[i*n+i] + T(1);
          if (i == 0)
             break;
       }
    
       for (k=n-1; ; k--)
       {
          for (its=1; its <= 30; its++)
          {
             flag = 1;
             for (l=k; ; l--)
             {
                nm = l-1;
                if (abs( rv1[l]) < epsilon(rv1[l]))
                {
                   flag = 0;
                   break;
                }
                if (abs( w[nm]) < epsilon(w[nm]))
                   break;
                if (l == 0)
                   break;
             }
             if (flag)
             {
                c = T(0);
                s = T(1);
                for (i=l; i <= k; i++)
                {
                   f = s * rv1[i];
                   rv1[i] = c * rv1[i];
                   if (abs(f) < epsilon(f))
                      break;
                   g = w[i];
                   h = sqrt( f * f + g * g);
                   w[i] = h;
                   h = VT(1) / h;
                   c = g * h;
                   s = -f * h;
                   for (j=0; j < m; j++)
                   {
                      y = a[j*n+nm];
                      z = a[j*n+i];
                      a[j*n+nm] = y * c + z * s;
                      a[j*n+i] = z * c - y * s;
                   }
                }
             }
             z = w[k];
             if (l == k)                 
             {
                if (z.real() < VT(0))
                {
                   w[k] = -z;
                   for (j=0; j < n; j++)
                      v[j*n+k] = -v[j*n+k];
                }
                break;
             }
    
             if (its == 30)
                return false;
    
             x = w[l];
             nm = k - 1;
             y = w[nm];
             g = rv1[nm];
             h = rv1[k];
             f = ((y-z) * (y+z) + (g-h) * (g+h)) / (VT(2) * h * y);
             g = sqrt( f * f + T(1) * T(1));
             f = ((x-z) * (x+z) + h * ((y / (f+sign( g, f))) - h)) / x;
             c = s = T(1);                    
             for (j=l; j <= nm; j++)
             {
                i = j + 1;
                g = rv1[i];
                y = w[i];
                h = s * g;
                g = c * g;
                z = sqrt( f * f + h * h);
                rv1[j] = z;
                c = f / z;
                s = h / z;
                f = x * c + g * s;
                g = g*c - x*s;
                h = y * s;
                y *= c;
                for (jj=0; jj < n; jj++)
                {
                   x = v[jj*n+j];
                   z = v[jj*n+i];
                   v[jj*n+j] = x * c + z * s;
                   v[jj*n+i] = z * c - x * s;
                }
                z = sqrt( f * f + h * h);
                w[j] = z;
                if (abs(z) > epsilon(z))
                {
                   z = VT(1) / z;
                   c = f * z;
                   s = h * z;
                }
                f = c * g + s * y;
                x = c * y - s * g;
                for (jj=0; jj < m; jj++)
                {
                   y = a[jj*n+j];
                   z = a[jj*n+i];
                   a[jj*n+j] = y * c + z * s;
                   a[jj*n+i] = z * c - y * s;
                }
             }
             rv1[l] = 0.0;
             rv1[k] = f;
             w[k] = x;
          }
          if (k == 0)
             break;
       }
       return true;
    }
    
    template <> inline
    complex<float> matrix<complex<float> >::normF () const
    {
       float nr = 0.0f;
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
       {
          float e = abs( mPtr->val[i]);
          nr += e * e;
       }
       nr = sqrt( nr);
    
       return complex<float>( nr, 0.0f);
    }
    
    template <> inline
    complex<double> matrix<complex<double> >::normF () const
    {
       double nr = 0.0;
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
       {
          double e = abs( mPtr->val[i]);
          nr += e * e;
       }
       nr = sqrt( nr);
    
       return complex<double>( nr, 0.0);
    }
    
    template <> inline
    complex<long double> matrix<complex<long double> >::normF () const
    {
       long double nr = 0.0;
       size_t n = size();
    
       for (size_t i=0; i < n; i++)
       {
          long double e = abs( mPtr->val[i]);
          nr += e * e;
       }
       nr = sqrt( nr);
    
       return complex<long double>( nr, 0.0);
    }
    
    template <> inline
    complex<float> matrix<complex<float> >::norm1 () const
    {
       typedef value_type T;
       float nr = 0.0;
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       const T *pm = &mPtr->val[0];
    
       for (size_t j=0; j < n; j++)
       {
          float s = 0.0f;
          for (size_t i=0; i < m; i++)
             s += abs( pm[i*n+j]);
          nr = techsoft::max( nr, s);
       }
       return T( nr, 0.0f);
    }
    
    template <> inline
    complex<double> matrix<complex<double> >::norm1 () const
    {
       typedef value_type T;
       double nr = 0.0;
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       const T *pm = &mPtr->val[0];
    
       for (size_t j=0; j < n; j++)
       {
          double s = 0.0;
          for (size_t i=0; i < m; i++)
             s += abs( pm[i*n+j]);
          nr = techsoft::max( nr, s);
       }
       return T( nr, 0.0);
    }
    
    template <> inline
    complex<long double> matrix<complex<long double> >::norm1 () const
    {
       typedef value_type T;
       long double nr = 0.0;
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       const T *pm = &mPtr->val[0];
    
       for (size_t j=0; j < n; j++)
       {
          long double s = 0.0;
          for (size_t i=0; i < m; i++)
             s += abs( pm[i*n+j]);
          nr = techsoft::max( nr, s);
       }
       return T( nr, 0.0);
    }
    
    template <>  inline
    complex<float> matrix<complex<float> >::norm2 () const
    {
       matrix<complex<float> > m(*this);
       matrix<complex<float> > v(mPtr->ncol,mPtr->ncol);
    
       size_t i, n = mPtr->ncol;
       float nr = 0.0f;
       valarray<complex<float> > w( n);
    
       if (m.svd( v, w))
          for (i=0; i < n; i++)
             nr = techsoft::max( nr, abs( w[i]));
    
       return complex<float>( nr, 0.0f);
    }
    
    template <>  inline
    complex<double> matrix<complex<double> >::norm2 () const
    {
       matrix<complex<double> > m(*this);
       matrix<complex<double> > v(mPtr->ncol,mPtr->ncol);
    
       size_t i, n = mPtr->ncol;
       double nr = 0.0;
       valarray<complex<double> > w( n);
    
       if (m.svd( v, w))
          for (i=0; i < n; i++)
             nr = techsoft::max( nr, abs( w[i]));
    
       return complex<double>( nr, 0.0);
    }
    
    template <>  inline
    complex<long double> matrix<complex<long double> >::norm2 () const
    {
       matrix<complex<long double> > m(*this);
       matrix<complex<long double> > v(mPtr->ncol,mPtr->ncol);
    
       size_t i, n = mPtr->ncol;
       long double nr = 0.0;
       valarray<complex<long double> > w( n);
    
       if (m.svd( v, w))
          for (i=0; i < n; i++)
             nr = techsoft::max( nr, abs( w[i]));
    
       return complex<long double>( nr, 0.0);
    }
    
    template <> inline
    complex<float> matrix<complex<float> >::normI () const
    {
       typedef value_type T;
       float nr = 0.0f;
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       const T *pm = &mPtr->val[0];
    
       for (size_t i=0; i < m; i++)
       {
          float s = 0.0f;
          for (size_t j=0; j < n; j++)
             s += abs( pm[i*n+j]);
          nr = techsoft::max( nr, s);
       }
       return T( nr, 0.0f);
    }
    
    template <> inline
    complex<double> matrix<complex<double> >::normI () const
    {
       typedef value_type Ty;
       double nr = 0.0;
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       const Ty *pm = &mPtr->val[0];
    
       for (size_t i=0; i < m; i++)
       {
          double s = 0.0;
          for (size_t j=0; j < n; j++)
             s += abs( pm[i*n+j]);
          nr = techsoft::max( nr, s);
       }
       return Ty( nr, 0.0);
    }
    
    template <> inline
    complex<long double> matrix<complex<long double> >::normI () const
    {
       typedef value_type Ty;
       long double nr = 0.0;
       size_t m = mPtr->nrow;
       size_t n = mPtr->ncol;
       const Ty *pm = &mPtr->val[0];
    
       for (size_t i=0; i < m; i++)
       {
          long double s = 0.0;
          for (size_t j=0; j < n; j++)
             s += abs( pm[i*n+j]);
          nr = techsoft::max( nr, s);
       }
       return Ty( nr, 0.0);
    }
    
    
    } // namespace techsoft
    
    #endif // _MATRIX_CC

--------------------

## Test.cpp ##

/************************************************************************
    * Copyright(c) 2011  Gavin Liu
    * All rights reserved.
    *
    * File: test.cpp
    * Brief: 
    * Version: 1.0
    * Author: Gavin Liu
    * Email: gavinliuyi@msn.com
    * Date: 2011/08/28
    * History:
    ************************************************************************/
    
    #include <cstdio>
    #include <cstdlib>
    #include <iostream>
    #include <fstream>
    #include <iomanip>
    #include"matrix.cpp"
    
    #if defined(__BORLANDC__)
    #pragma option -w-aus
    #endif
    
    //
    // Using declarations. 
    //
    using std::cout;
    using std::endl;
    using techsoft::epsilon;
    using techsoft::isVecEq;
    using std::exception;
    
    #ifdef _MSC_VER
    using ::rand;
    #else
    using std::rand;
    using std::size_t;
    #endif
    
    
    //
    // typedef for different matrix types
    //
    typedef techsoft::matrix<float>        fMatrix;
    typedef std::valarray<float>           fVector;
    
    typedef techsoft::matrix<double>       dMatrix;
    typedef std::valarray<double>          dVector;
    
    typedef techsoft::matrix<long double>  ldMatrix;
    typedef std::valarray<long double>     ldVector;
    
    typedef std::complex<float>            fComplex;
    typedef techsoft::matrix<fComplex>     cfMatrix;
    typedef std::valarray<fComplex>        cfVector;
    
    typedef std::complex<double>           dComplex;
    typedef techsoft::matrix<dComplex>     cdMatrix;
    typedef std::valarray<dComplex>        cdVector;
    
    typedef std::complex<long double>      ldComplex;
    typedef techsoft::matrix<ldComplex>    cldMatrix;
    typedef std::valarray<ldComplex>       cldVector;
    
    const char *TestFile = "test.txt";
    
    void test_const ();     // Test constructor/destructor
    void test_submat ();    // Test sub-matrix operations
    void test_op ();        // Test unary operators
    void test_util ();      // Test utility methods
    void test_dcomp ();     // Test decomposition methods
    void test_eigen ();     // Test eigen value and eigen vector
    void test_io ();           // Test I/O
    
    template <class T> T 
    mfn (size_t i, size_t j, const T& v)
    {
       int k = int(i+j);
       return pow( v, k);
    }
    
    
    int main ()
    {
       #if !defined(_MSC_VER)
       using std::remove;
       using std::srand;
       #endif
    
       try
       {
          srand( 0x23657876);
          test_const();
          test_submat();
          test_op();
          test_util();
          test_dcomp();
          test_eigen();
          test_io();
          remove( TestFile);
       }
       catch (const exception& e)
       {
          cout << "Error: " << e.what();
       }
       return 0;
    }
    
    void test_const ()     // Test constructor/destructor
    {
       try
       {
          cout << "Testing constructors.....";
          dMatrix m(3,3);
          dMatrix m2(3,3,1.0);
    
          m = 1.0;
          if (m != m2)
          {
             cout << "failed!\n";
             return;
          }
    
          float fv[] = { 1.34f, 2.54f, 8.23f, 
                         7.34f, -2.3f, 2.45f,
                         -1.2f, 4.50f, 7.34f };
    
          fMatrix mf( 3, 3, fv);
          fMatrix mf2( 3, 3, fv, techsoft::FORTRAN_ARRAY);
          fMatrix mf3 = mf;
          fVector vf( fv,9);
          fMatrix mf4( 3,3,vf);
    
          if (mf != ~mf2 || mf != mf3 || mf != mf4)
          {
             cout << "failed!\n";
             return;
          }
    
          dMatrix ma[5];
          for (size_t i=0; i < 5; i++)
          {
             ma[i].resize( 3, 3);
             if (!ma[i].isNull())
             {
                cout << "failed!\n";
                return;
             }
          }
    
          cfMatrix cm(4,4);
          cm.rand();
          cfMatrix cm2 = cm;
    
    #ifndef _TS_NO_MEMBER_TEMPLATES
          dMatrix md = mf;
          cdMatrix mcd = md;
          mcd = 1.0;
          mcd = cm;
    #endif
       }
       catch (const exception& e)
       {
          cout << "failed!\nError: " << e.what() << endl;
          return;
       }
       cout << "Ok\n";
    }
    
    
    void test_submat ()
    {
       using techsoft::mslice;
       using techsoft::gmslice;
       using techsoft::mswap;
       using techsoft::abs;
       using techsoft::DIAGONAL;
       using techsoft::TRIDIAGONAL;
       using techsoft::UTRIANG;
       using techsoft::LTRIANG;
    
       try
       {
          cout << "Testing subscript/submatrix.....";
    
          dMatrix m(3,3);
          m.rand();
          const dMatrix m2 = m;
          double x,y;
          size_t i=0,j=0;
    
          x = m2[i][j];
          m[i][j] = x;
          y = m[i][j];
          if (abs(x-y) > epsilon(x))
          {
             cout << "failed!\n";
             return;
          }
          x = m(i,j)++;
          y = ++m(i,j);
          if (abs(y-x-2.0) > epsilon(y))
          {
             cout << "failed!\n";
             return;
          }
          x = m(i,j)--;
          y = --m(i,j);
          if (abs(x-y-2.0) > epsilon(x))
          {
             cout << "failed!\n";
             return;
          }
          x = m[i][j]  * m[i+1][j+1];
          x = m2[i][j] / m2[i+1][j+1];
          x = m[i][j]  + m2[i+1][j+1];
          x = m2[i][j] - m[i+1][j+1];
    
          y = x * m[i][j];
          y = m[i][j] * x;
          y = x / m[i][j];
          y = m[i][j] / x;
          y = x + m[i][j];
          y = m[i][j] + x;
          y = x - m[i][j];
          y = m[i][j] - x;
    
    #if defined(_MSC_VER) && _MSC_VER > 1100     // MS VC++ 5.0 is generating wrong code!!!
          x = m[i][j];
          m[i][j] += m[i+1][j+1];
          m[i][j] -= m[i+1][j+1];
          m[i][j] *= m[i+1][j+1];
          m[i][j] /= m[i+1][j+1];
          y = m[i][j];
          if (abs(y-x) > epsilon(y))
          {
             cout << "failed!\n";
             return;
          }
    #endif
    
          x = -m[i][j];
          y = +m[i][j];
          if (abs(y+x) > epsilon(x+y))
          {
             cout << "failed!\n";
             return;
          }
    
          dMatrix a(6,6);
          a.rand();
          const dMatrix b = a;
    
          dMatrix sa = b[mslice(1,1,3,3)];
          a[mslice(1,1,3,3)] = sa;
          if (a != b)
          {
             cout << "failed!\n";
             return;
          }
    
          a[mslice(2,2,1,3)] = b[mslice(2,2,1,3)];
          if (a != b)
          {
             cout << "failed!\n";
             return;
          }
    
          a[mslice(1,0,2,3)] *= 10.0;
          a[mslice(1,0,2,3)] /= 10.0;
          if (a != b)
          {
             cout << "failed!\n";
             return;
          }
    
          a[mslice(1,0,2,3)] += b[mslice(1,0,2,3)];
          a[mslice(1,0,2,3)] -= b[mslice(1,0,2,3)];
          if (a != b)
          {
             cout << "failed!\n";
             return;
          }
    
          a[mslice(0,0,6,6)] = 1.0;
          dMatrix c(6,6,1.0);
          if (a != c)
          {
             cout << "failed!\n";
             return;
          }
    
          a.resize(4,4);
          a.rand();
          c = a[mslice(1,1,2,2)];
          c = a[gmslice(DIAGONAL)];
          c = a[gmslice(TRIDIAGONAL)];
          c = a[gmslice(UTRIANG)];
          c = a[gmslice(LTRIANG)];
    
          c = 1.0;
          c[gmslice(DIAGONAL)] += a;   
          c[gmslice(TRIDIAGONAL)] -= a;      
          c[gmslice(UTRIANG)] *= 12.0;
          c[gmslice(LTRIANG)] /= 12.0;
    
          dMatrix cm = a[gmslice(UTRIANG)];
          const dMatrix ca = a;
          cm = ca[gmslice(UTRIANG)];
    
          dMatrix d(4,4), d2(4,4);
          d.rand();
          d2.rand();
          mswap( d[0][0], d[1][1]);
          mswap( d[2], d[3]);
          mswap( d(1), d(2));
    
          mswap( d[0][0], d2[0][0]);
          mswap( d[2], d2[2]);
          mswap( d(1), d2(1));
          mswap( d, d2);      
    
          a = b;
          dVector v = b[5];
          a[5] = v;
          v = a(3);
          a(3) = v;
          a[2] += v;
          a[2] -= v;
          a[4] *= v;
          a[4] /= v;
          if (a != b)
          {
             cout << "failed!\n";
             return;
          }
    
          v = b.diag();
          a.diag() = v;
    
          dVector v2 = b.diag(-2);
          a.diag(-2) = v2;
          v2 = b.diag(2);
          a.diag(2) = v2;
          if (a != b)
          {
             cout << "failed!\n";
             return;
          }
       }
       catch (const exception& e)
       {
          cout << "failed!\nError: " << e.what() << endl;
          return;
       }
       cout << "Ok\n";
    }
    
    void test_op ()
    {
       try
       {
          cout << "Testing operators.....";
          dMatrix m(6,4);
    
          m.rand();
          dMatrix m2 = m;
          m2 += m;
          m2 -= m;
          if (m2 != m)
          {
             cout << "failed!\n";
             return;
          }
          double te(1.234567f);
          m2 *= te;
          m2 /= te;
          if (m2 != m)
          {
             cout << "failed!\n";
             return;
          }
    
          dMatrix a(4,4);
          a.rand();
          dMatrix b = !a;
          a *= b;
          if (!a.isUnit())
          {
             cout << "failed!\n";
             return;
          }
    
          a.rand();
          a /= a;
          if (!a.isUnit())
          {
             cout << "failed!\n";
             return;
          }
    
          a.resize( 5,3);
          a.rand();
          b = ~a;
          dMatrix c = a * b;
          b = ~b;
          if (a != b)
          {
             cout << "failed!\n";
             return;
          }
          c = a + b;
          b = c - a;
          if (a != b)
          {
             cout << "failed!\n";
             return;
          }
    
          a.resize(4,4);
          a.rand();
          b = !a;
          c = a * b;
          if (!c.isUnit())
          {
             cout << "failed!\n";
             return;
          }
    
          double x = 2.045f;
          c = 1/x * a * x;
          if (a != c)
          {
             cout << "failed!\n";
             return;
          }
    
          c = x * a / x;
          if (a != c)
          {
             cout << "failed!\n";
             return;
          }
    
          c = x / a;
          b = x * !a;
          if (b != c)
          {
             cout << "failed!\n";
             return;
          }
    
          dVector v(4), v2(4);
          v = b[0];
          a.unit();
          v2 = v * a;
          v = a * v2;
          if (!isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
          v2 = v / a;
          v = a / v;
       }
       catch (const exception& e)
       {
          cout << "failed!\nError: " << e.what() << endl;
          return;
       }
       cout << "Ok\n";
    }
    
    void test_util ()    // Test utility methods
    {
       using techsoft::pow;
       using techsoft::abs;
       using techsoft::floor;
       using techsoft::ceil;
       using techsoft::gmslice;
       using techsoft::UTRIANG;
       using techsoft::LTRIANG;
    
       try
       {
          cout << "Testing utility methods.....";
          dMatrix m(6,4);
    
          m.rand();
          dMatrix m2 = m;
          m.unit();
          m.resize( 2, 3);
          m.unit();
          m.resize( 4, 6);
          m.unit();
          m.free();
    
          m.resize( 4, 4);
          m.rand();
          m2 = m * m * m * m * m * m * m;
          if (m2 != pow( m, 7))
          {
             cout << "failed!\n";
             return;
          }
    
          m2 = -m;
          m2 = abs( m2);
    
          m2 = floor( m2);
          if (!m2.isNull() || m2.trace() != 0.0)
          {
             cout << "failed!\n";
             return;
          }
          m2 = ceil( m);
    
          double val;
          size_t rn;
    
          m.rand();
          val = m.norm1();
          val = m.norm2();
          val = m.normI();
          val = m.normF();
          val = m.cond();
          val = m.det();
          rn = m.rank();
          val = m.cofact(0,0);
    
          m2 = m.adj()/m.det();
          if (m2 != !m)
          {
             cout << "failed!\n";
             return;
          }
    
          m2 = m;
          m2.inv_lu();
          if (m2 != !m)
          {
             cout << "failed!\n";
             return;
          }
    
          m2 = m;
          m2.inv_sv();
          if (m2 != !m)
          {
             cout << "failed!\n";
             return;
          }
    
          m2 = m;
          m2.inv_qr();
          if (m2 != !m)
          {
             cout << "failed!\n";
             return;
          }
    
          cdMatrix cm(4,4), cm2(4,4);
          cm.rand();
          cm2 = cm;
          cm2.inv_sv();
          if (cm2 != !cm)
          {
             cout << "failed!\n";
             return;
          }
    
          bool bt = m.isSingular();
          m = double(0);
          bt = m.isNull();
          m.diag() = 1.0;
          bt = m.isUnit();
          m.diag() = 3.0;
          bt = m.isScalar();
          m.diag()[0] = 7.0;
          bt = m.isDiagonal();
    
          m.rand();
          m2 = m[gmslice(UTRIANG)];
          if (!m2.isUpperTriangular())
          {
             cout << "failed!\n";
             return;
          }
          m2 = m[gmslice(LTRIANG)];
          if (!m2.isLowerTriangular())
          {
             cout << "failed!\n";
             return;
          }
    
          dMatrix m3 = ~m2;
          m2[gmslice(UTRIANG)] = m3;
          if (!m2.isSymmetric())
          {
             cout << "failed!\n";
             return;
          }
          m2[gmslice(UTRIANG)] = -m3;
          if (!m2.isSkewSymmetric())
          {
             cout << "failed!\n";
             return;
          }
    #if defined(_MSC_VER) && _MSC_VER > 1100
          m2 = m.apply( log);
          m2 = m.apply( mfn);
    #endif
       }
       catch (const exception& e)
       {
          cout << "failed!\nError: " << e.what() << endl;
       }
       cout << "Ok\n";
    
    }
    
    void test_dcomp ()     // Test constructor/destructor
    {
       using techsoft::gmslice;
       using techsoft::LTRIANG;
    
       try
       {
          cout << "Testing matrix decomposition.....";
          size_t i, n = 4;
          dMatrix a(n,n);
          dVector v(n),s(n),v2(n), w(n);
          std::valarray<size_t> indx(n);
    
          a.rand();
          dMatrix b = a;
          v = a[0];
          v *= double(1.2394f);
    
          b.lud( indx);
          b.lubksb( indx, v, s);
          v2 = a * s;
          if (!isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          a.lumpove( b, indx, v, s);
          v2 = a * s;
          if (!isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          v /= double(-1.7913f);
          b.lubksb( indx, v, s);
          v2 = a * s;
          if (!isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          dMatrix mv(n,3), ms(n,3), mv2(n,3);
          mv.rand();
          a.solve( mv, ms);
          for (i=0; i < 3; i++)
          {
             dVector vx = ms(i);
             mv2(i) = a * vx;
          }
          if (mv != mv2)
          {
             cout << "failed!\n";
             return;
          }
    
          b = a[gmslice(LTRIANG)];
          b *= ~b;
          a = b;
          b.chold();
          b.cholbksb( v, s);
          v2 = a * s;
          if (a != (b * ~b) || !isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          a.solve_chol( v, s);
          v2 = a * s;
          if (!isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          a.rand();
          b = a;
          dMatrix c(n,n);
          b.qrd( c);
          b.qrbksb( c, v, s);
          v2 = a * s;
          if (!b.isRowOrthogonal() || a != (b * c) || !isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          a.solve_qr( v, s);
          v2 = a * s;
          if (!isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          b = a;
          b.svd( c, w);
          b.svbksb( c, w, v, s);
          v2 = a * s;
          dMatrix d(n,n,0.0);
          d.diag() = w;
          if (!b.isColOrthogonal() || !c.isRowOrthogonal()
              || a != (b * d * ~c) || !isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          a.solve_sv( v, s);
          v2 = a * s;
          if (!isVecEq( v, v2))
          {
             cout << "failed!\n";
             return;
          }
    
          a.resize( n+2,n);
          a.rand();
          b = a;
          b.qrd( c);
          v.resize( n+2);
          for (i=0; i < n+2; i++)
             v[i] = ((double)rand())/RAND_MAX;
          b.qrbksb( c, v, s);
          v2 = a * s;
          if (a != b * c)
          {
             cout << "failed!\n";
             return;
          }
    
          b = a;
          b.svd( c, w);
          b.svbksb( c, w, v, s);
          d.diag() = w;
    
          if (a != (b * d * ~c)
              || !b.isColOrthogonal() || !c.isRowOrthogonal())
          {
             cout << "failed!\n";
             return;
          }
    
          cdMatrix ca(n,n);
          cdVector cv(n),cs(n),cv2(n);
    
          ca.rand();
          for (i=0; i < n; i++)
          {
             double ta = ((double)rand())/RAND_MAX;
             double tb = ((double)rand())/RAND_MAX;
             dComplex z(ta,tb);
             cv[i] = z;
          }
    
          cdMatrix cb = ca;
          cdMatrix cd = !ca;
          cs = cd * cv;
          cv2 = ca * cs;
          if (!isVecEq( cv, cv2))
          {
             cout << "failed!\n";
             return;
          }
    
          cb.lud( indx);
          cb.lubksb( indx, cv, cs);
          cv2 = ca * cs;
          if (!isVecEq( cv, cv2))
          {
             cout << "failed!\n";
             return;
          }
       }
       catch (const exception& e)
       {
          cout << "failed!\nError: " << e.what() << endl;
       }
       cout << "Ok\n";
    }
    
    
    void test_eigen ()
    {
       using techsoft::gmslice;
       using techsoft::UTRIANG;
       using techsoft::LTRIANG;
    
       try
       {
          cout << "Testing eigen value.....";
    
          size_t n = 4;
          dMatrix a(n,n), ev(n,n);
          dVector d(n), e(n);
          bool ret,ret2;
    
          a.rand();
          a = a[gmslice(UTRIANG)];
          a[gmslice(LTRIANG)] = ~a;
    
          ret = a.eigen( e);
          ret2 = a.eigen( d, ev);
          if (!ret || !ret2 || !isVecEq( d, e))
          {
             cout << "failed!\n";
             return;
          }
    
          dMatrix c(n,n,0.0);
          c.diag() = d;
          if ((a * ev) != (ev * c))
          {
             cout << "failed!\n";
             return;
          }
    
          a.rand();
          ret = a.eigen( d, e);
          ret2 = a.eigen( d, e, ev);
          if (!ret || !ret2)
          {
             cout << "failed!\n";
             return;
          }
    
          for (size_t i=0; i < n; i++)
          {
             c(i,i) = d[i];
             if (e[i] > double(0) && i < n-1)
                c(i,i+1) = e[i];
             else if (e[i] < double(0) && i > 0)
                c(i,i-1) = e[i];
          }
          if ((a * ev) != (ev * c))
          {
             cout << "failed!\n";
             return;
          }
       }
       catch (const exception& e)
       {
          cout << "failed!\nError: " << e.what() << endl;
       }
       cout << "Ok\n";
    }
    
    
    void test_io ()
    {
       using techsoft::abs;
       #if !defined(_MSC_VER)
       using std::FILE;
       using std::fopen;
       using std::fclose;
       using std::fread;
       using std::fwrite;
       #endif
    
       size_t i,m=4,n=3;
       dMatrix a(m,n), b(m,n);
       dVector v(m), v2(m);
       double tmp;
       cout << "Testing matrix I/O.....";
       a.rand();
    
       std::ofstream outFile( TestFile);
       if (!outFile)
       {
          std::cerr << "\nError opefing test file.";
           return;
       }
       outFile.clear();
    
       std::ifstream inFile( TestFile);
       if (!inFile)
       {
           std::cerr << "\nError opefing test file.";
           return;
       }
       inFile.tie( &outFile);
       outFile << std::scientific << std::setprecision( 18);
    
       outFile << a[0][0] << endl;
       inFile >> b[0][0];
       tmp = a[0][0] - b[0][0];
       if (abs( tmp) > epsilon(tmp))
       {
          cout << "failed!\n" << a[0][0] << endl << b[0][0] ;
          return;
       }
    
       const dMatrix c = a;
       outFile << c(1,1) << endl;
       inFile >> b(1,1);
       tmp = c(1,1) - b(1,1);
       if (abs( tmp) > epsilon( tmp))
       {
          cout << "failed!\n";
          return;
       }
    
       v = a(0);
       outFile << a(0) << endl;
       inFile >> a(0);
       v2 = a(0);
       a(0) = v;
    
       if (!isVecEq( v, v2))
       {
          cout << "failed!\n";
          return;
       }
    
       outFile << a[1] << endl;
       inFile >> b[1];
       for (i=0; i < n; i++)
       {
          tmp = a[1][i] - b[1][i];
          if (abs( tmp) > epsilon(tmp))
          {
             cout << "failed!\n";
             return;
          }
       }
    
       outFile << c[1] << endl;
       inFile >> b[1];
       for (i=0; i < n; i++)
       {
          tmp = c[1][i] - b[1][i];
          if (abs( tmp) > epsilon(tmp))
          {
             cout << "failed!\n";
             return;
          }
       }
    
       outFile << a << endl;
       inFile >> b;
       if (a != b)
       {
          cout << "failed!\n";
          return;
       }
       outFile.close();
       inFile.close();
    
       FILE *fp = fopen( TestFile, "w+b");
       if (fp == NULL)
       {
          std::cerr << "\nError opefing test file.";
           return;
       }
    
       a.rand();
       fwrite( &a[0][0], a.size(), a.typesize(), fp);
       rewind( fp);
       fread( &b[0][0], b.size(), b.typesize(), fp);
       fclose( fp);
       if (a != b)
       {
          cout << "failed!\n";
          return;
       }
        cout << "Ok\n";
    }

--------------------

## Makefile ##

matrix : test.cpp
        g++ -o matrix test.cpp matrix.cpp 
    clean:
        rm -f matrix