Matrix.h
#ifndef _MATRIX_H
#define _MATRIX_H
#if defined(DEBUG) || defined(_DEBUG)
#define RANGE_CHECK_
#endif
#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
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;
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 ();
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
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; }
size_t rowsize ()
const throw() {
return mPtr->ncol; }
size_t rowno ()
const throw() {
return mPtr->nrow; }
size_t colno ()
const throw() {
return mPtr->ncol; }
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;
matrix<T>
operator+ ()
const {
return *
this; }
matrix<T>
operator- ()
const;
matrix<T>
operator~ ()
const;
matrix<T>
operator! ()
const;
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);
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;
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;
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>;
};
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);
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);
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;
};
}
#endif
Matrix.cpp
#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;
}
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()];
}
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;
}
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;
}
}
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;
}
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
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;
}
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
{
++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;
}
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();
}
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])
{
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;
}
}
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;
}
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];
}
}
#ifndef MATRIX_ROW_SEP
#define MATRIX_ROW_SEP ('\n')
#endif
#ifndef MATRIX_COL_SEP
#define MATRIX_COL_SEP ('\t')
#endif
template <
class T>
inline istream&
operator>> (istream& is, Mref<T> el)
{
T x; is >> x; el = x;
return is;
}
template <
class T>
inline ostream&
operator<< (ostream& os,
const Mref<T>& el)
{
const T x = el;
os << x;
return os;
}
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;
}
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;
}
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;
}
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;
}
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;
}
template <
class T> istream&
operator>> (istream& is, valarray<T>& v)
{
for (size_t i=
0; i < v.size(); i++)
is >> v[i];
return is;
}
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;
}
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);
}
}
#endif
Test.cpp
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <iomanip>
#include"matrix.cpp"
#if defined(__BORLANDC__)
#pragma option -w-aus
#endif
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 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 ();
void test_submat ();
void test_op ();
void test_util ();
void test_dcomp ();
void test_eigen ();
void test_io ();
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 ()
{
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
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 ()
{
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 ()
{
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
