本来是写一个模版矩阵类,后来有添加了几个线性代数相关的函数,高斯消去法和LU分解;刚刚又加了一个QR分解函数,还想继续添加……
/************************************************************************/
/* 模版矩阵类,支持基本的数据类型,和stl的复数类,如果自定义的类型需 */
/* 要支持运算符:+,-,*,/,+=,-=,*=,/=,=,常用的数学函数 */
/************************************************************************/
/* By JJJ */
/* QQ:771588193 */
/* 2013/11/23 */
/************************************************************************/
#pragma once
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <assert.h>
using namespace std;
namespace myblas
{
/************************************************************************/
/* 简单的动态数组 */
/************************************************************************/
template<typename Ty>
class basic_vector
{
public:
basic_vector():v_data(NULL),v_len(0) {}
explicit
basic_vector(int n):v_data(NULL),v_len(0) {
alloc(n);
}
basic_vector(const basic_vector &vec):v_data(NULL),v_len(0) {
assign(vec);
}
basic_vector(basic_vector &&_tmp):v_data(_tmp.v_data),v_len(_tmp.v_len) {
_tmp.v_data = NULL; _tmp.v_len = 0;
}
~basic_vector() {
release();
}
//
inline void alloc(int n)
{
v_data = new Ty[n]();
if (NULL!=v_data)
v_len = n;
}
inline void release()
{
if (NULL!=v_data)
{
delete[] v_data;
v_data = NULL;
v_len = 0;
}
}
bool assign(const basic_vector &vec)
{
if (NULL!=vec.v_data)
{
alloc(vec.v_len);
if (NULL!=v_data)
{
int i;
for (i=0; i<v_len; i++)
v_data[i] = vec.v_data[i];
return true;
}
}
return false;
}
inline int get_len()
{
return v_len;
}
//
basic_vector& operator = (const basic_vector &vec)
{
assign(vec);
}
inline Ty &operator () (int i)
{
#ifdef _CHECK_BOARDER_
assert(0<=i);
assert(i<v_len);
#endif
return v_data[i];
}
//
protected:
Ty *v_data;
int v_len;
};
typedef basic_vector<int> array_int;
ostream& operator << (ostream &out,array_int &arr)
{
int i,j;
int n = arr.get_len();
out<<endl;
for (i=0; i<n; i++)
{
for (j=0; j<n; j++)
if (arr(i)==j)
out<<1<<'\t';
else
out<<0<<'\t';
out<<endl;
}
return out;
}
//
/************************************************************************/
/* 矩阵模版类 */
/************************************************************************/
typedef enum matrix_type {
NONE =0, IDENTIFY, ONEVAL, UPPER_TRIANGULAR, LOWER_TRIANGULAR
}Mat_Ty;
template<typename T>
class basic_matrix
{
public:
//
//construct and destruct
basic_matrix():m_data(NULL),m_row(0),m_col(0) {}
explicit
basic_matrix(int n,Mat_Ty mty =NONE,T val =T(1)):m_data(NULL),m_row(0),m_col(0) {
if (allocate(n,n)) set_type(mty,val);
}
explicit
basic_matrix(int r,int c,Mat_Ty mty =NONE,T val =T(1)):m_data(NULL),m_row(0),m_col(0) {
if (allocate(r,c)) set_type(mty,val);
}
basic_matrix(const basic_matrix &src):m_data(NULL),m_row(0),m_col(0) {
assign(src);
}
explicit
basic_matrix(const basic_matrix &src, int r1,int c1,int r2,int c2):m_data(NULL),m_row(0),m_col(0) {
assign(src,r1,c1,r2,c2);
}
basic_matrix( basic_matrix &&_right):m_data(_right.m_data),m_row(_right.m_row),m_col(_right.m_col) {
_right.m_data = NULL; _right.m_row = 0; _right.m_col = 0;
}
~basic_matrix() {
release();
}
//
//manage memory and data
inline bool allocate(int r, int c)
{
if (r>0 && c>0)
{
T *tmp = new T[r * c]();
if (NULL!=tmp)
{
m_data = tmp;
m_row = r;
m_col = c;
return true;
}
}
return false;
}
inline void release()
{
if (NULL!=m_data)
{
delete[] m_data;
m_data = NULL;
m_row = m_col = 0;
}
}
inline void resize(int r,int c,bool Keep_oldvalue =false)
{
if (r<=0 || c<=0)
return;
if (r*c>m_row*m_col)
{
T *tmp = new T[r*c]();
if (NULL!=tmp)
{
if (Keep_oldvalue)
for (int i=0; i<m_row*m_col; i++)
tmp[i] = m_data[i];
delete[] m_data;
m_data = tmp;
m_row = r;
m_col = c;
}
}
else
{
m_row = r;
m_col = c;
}
}
bool assign(const basic_matrix &other)
{
if (NULL!=other.m_data)
{
if (NULL==m_data && (m_row*m_col)<(other.m_row*other.m_col))
{
release();
allocate(other.m_row,other.m_col);
}
else
{
m_row = other.m_row;
m_col = other.m_col;
}
if (NULL!=m_data)
{
int i;
for (i=0; i<m_row*m_col; ++i)
m_data[i] = other.m_data[i];
return true;
}
}
return false;
}
bool assign(const basic_matrix &src,int r1,int c1,int r2,int c2)
{
bool test = ( NULL!=src.m_data && r1<r2 && c1<c2 &&
r1>=0 && c1>=0 && r2<src.m_row && c2<src.m_col );
if (test)
{
int nr = r2 - r1 +1;
int nc = c2 - c1 +1;
if (NULL==m_data && (m_row*m_col)<(nr*nc))
{
release();
allocate(nr,nc);
}
else
{
m_row = nr;
m_col = nc;
}
if (NULL!=m_data)
{
int s_col = src.m_col;
const T *s_data = src.m_data;
int i(0),j(0),k(0);
for (i=r1-1; i<r2; ++i)
for (j=c1-1; j<c2; ++j)
m_data[k++] = s_data[i*s_col+j];
return true;
}
}
return false;
}
const T& max_abs()
{
int i,k;
for (k=0,i=1; i<m_row*m_col; i++)
if ( abs(m_data[k]) < abs(m_data[i]) )
k = i;
return m_data[k];
}
const T& min_abs()
{
int i,k;
for (k=0,i=1; i<m_row*m_col; i++)
if ( abs(m_data[k]) > abs(m_data[i]) )
k = i;
return m_data[k];
}
inline static void swap(basic_matrix &a, basic_matrix &b)
{
int t;
t = a.m_row; a.m_row = b.m_row; b.m_row = t;
t = a.m_col; a.m_col = b.m_col; b.m_col = t;
T *td;
td = a.m_data; a.m_data = b.m_data; b.m_data = td;
}
inline void swap_row(int r1,int r2)
{
#ifdef _CHECK_BOARDER_
assert(0<=r1 && r1<m_row);
assert(0<=r2 && r2<m_row);
#endif
if ( r1!=r2 )
{
int i;
T tmp;
for (i=0; i<m_col; i++)
{
tmp = m_data[r1*m_col+i];
m_data[r1*m_col+i] = m_data[r2*m_col+i];
m_data[r2*m_col+i] = tmp;
}
}
}
inline void swap_col(int c1,int c2)
{
#ifdef _CHECK_BOARDER_
assert(0<=c1 && c1<m_col);
assert(0<=c2 && c2<m_col);
#endif
if ( c1!=c2 )
{
int i;
T tmp;
for (i=0; i<m_row; i++)
{
tmp = m_data[i*m_col+c1];
m_data[i*m_col+c1] = m_data[i*m_col+c2];
m_data[i*m_col+c2] = tmp;
}
}
}
inline void set_one_val(T val)
{
int i;
for (i=0; i<(m_row * m_col); ++i)
{
m_data[i] = val;
}
}
inline void set_ident(T val)
{
int n = min(m_row,m_col);
set_one_val(T(0));
int i;
for (i=0; i<n; i++)
{
m_data[i*m_col+i] = val;
}
}
inline void set_type(Mat_Ty mty,T val =T(1))
{
switch (mty)
{
case NONE: break;
case ONEVAL: set_one_val(val); break;
case IDENTIFY: set_ident(val); break;
default: break;
}
}
inline int get_row()
{
return m_row;
}
inline int get_col()
{
return m_col;
}
inline int get_size()
{
return (m_row * m_col);
}
inline bool is_empty()
{
return (NULL==m_data);
}
//
// operators
inline T& operator () (int i, int j =0)
{
#ifdef _CHECK_BOARDER_
// just for debug
assert(NULL!=m_data);
assert(0<=i);
assert(i<m_row);
assert(0<=j);
assert(j<m_col);
#endif // _CHECK_BORDER_
return m_data[i*m_col+j];
}
basic_matrix &operator = (const basic_matrix &other)
{
assign(other);
return *this;
}
basic_matrix &operator = (basic_matrix &&_right)
{
swap(*this,_right);
return *this;
}
basic_matrix operator + (const basic_matrix &other)
{
return add_another(*this,other,true);
}
basic_matrix operator - (const basic_matrix &other)
{
return add_another(*this,other,false);
}
basic_matrix &operator += (const basic_matrix &other)
{
add_self(*this,other,true);
return (*this);
}
basic_matrix &operator -= (const basic_matrix &other)
{
add_self(*this,other,false);
return (*this);
}
basic_matrix operator * (const basic_matrix &other)
{
return mat_multi(*this, other);
}
basic_matrix operator * (const T &&alpha)
{
if (NULL!=m_data)
{
basic_matrix dest(m_row,m_col);
if (NULL!=dest.m_data)
{
for (int i=0; i<m_row*m_col; ++i)
{
dest.m_data[i] = m_data[i] * alpha;
}
}
return move(dest);
}
else
return basic_matrix();
}
friend basic_matrix operator * (const T &&alpha, const basic_matrix &src)
{
if (NULL!=src.m_data)
{
int r = src.m_row;
int c = src.m_col;
basic_matrix dest(r,c);
if (NULL!=dest.m_data)
{
for (int i=0; i<r*c; ++i)
{
dest.m_data[i] = alpha * src.m_data[i];
}
}
return move(dest);
}
else
return basic_matrix();
}
basic_matrix operator ^ (int n)
{
if (m_row != m_col)
{
cout<<"Inputs must be a scalar and a square matrix."<<endl;
return basic_matrix();
}
if (0==n)
{
return basic_matrix(m_row,m_col,IDENTIFY);
}
else if (1==n)
{
return (*this);
}
else if (-1==n)
{
return inverse(*this);
}
else if (1<n)
{
basic_matrix prod = mat_multi(*this,*this);
for (n-=2; n; n--)
prod = mat_multi(prod,*this);
return move(prod);
}
else // if (-1>n)
{
basic_matrix tmp = inverse(*this);
basic_matrix prod = mat_multi(tmp,tmp);
for (n+=2; n; n++)
prod = mat_multi(prod,tmp);
return move(prod);
}
}
basic_matrix operator ~()
{
return transpose(*this);
}
//
//math compute function
static basic_matrix transpose(const basic_matrix &a)
{
if (NULL!=a.m_data)
{
if (1==a.m_row && 1==a.m_col)
{
return basic_matrix(a);
}
else
{
basic_matrix b(a.m_col,a.m_row);
T *a_data = a.m_data;
T *b_data = b.m_data;
int I_loop = a.m_row;
int J_loop = a.m_col;
int i,j;
for (i=0; i<I_loop; ++i)
for (j=0; j<J_loop; ++j)
b_data[j*I_loop+i] = a_data[i*J_loop+j];
return move(b);
}
}
else
return basic_matrix();
}
// return a +/- b
static basic_matrix add_another(const basic_matrix &a,const basic_matrix &b,bool OP_ADD)
{
bool test = (
NULL!= a.m_data &&
NULL!=b.m_data &&
a.m_row==b.m_row &&
a.m_col==b.m_col);
if (test)
{
basic_matrix c(a.m_row,a.m_col);
if (NULL!=c.m_data)
{
const T *a_data = a.m_data;
const T *b_data = b.m_data;
T *c_data = c.m_data;
int i;
if (OP_ADD)
{
for (i=0; i<(a.m_row * a.m_col); ++i)
{
c_data[i] = a_data[i] + b_data[i];
}
}
else
{
for (i=0; i<(a.m_row * a.m_col); ++i)
{
c_data[i] = a_data[i] - b_data[i];
}
}
}
return move(c);
}
else
{
cout<<"The matrices dimensions must agree."<<endl;
return basic_matrix();
}
}
// a = a +/- b
static bool add_self(basic_matrix &a,const basic_matrix &b,bool OP_ADD)
{
bool test = (
NULL!= a.m_data &&
NULL!=b.m_data &&
a.m_row==b.m_row &&
a.m_col==b.m_col);
if (test)
{
T *a_data = a.m_data;
const T *b_data = b.m_data;
int i;
if (OP_ADD)
{
for (i=0; i<(a.m_row * a.m_col); ++i)
a_data[i] += b_data[i];
}
else
{
for (i=0; i<(a.m_row * a.m_col); ++i)
a_data[i] -= b_data[i];
}
}
else
cout<<"The matrices dimensions must agree."<<endl;
return test;
}
// compute c = a * b
static basic_matrix mat_multi(const basic_matrix &a,const basic_matrix &b)
{
bool test = (NULL!= a.m_data && NULL!=b.m_data && a.m_col==b.m_row );
if (test)
{
basic_matrix c(a.m_row,b.m_col);
c.set_one_val(T(0));
if (NULL!=c.m_data)
{
const T *A = a.m_data;
const T *B = b.m_data;
T *C = c.m_data;
int I_loop = a.m_row;
int J_loop = b.m_col;
int K_loop = a.m_col;
int i,j,k;
for (i=0; i<I_loop; ++i)
for (k=0; k<K_loop; ++k)
for (j=0; j<J_loop; ++j)
C[i*J_loop+j] += A[i*K_loop+k] * B[k*J_loop+j];
}
return move(c);
}
else
{
cout<<"Inner matrix dimensions must agree."<<endl;
return basic_matrix();
}
}
//
protected:
T *m_data; //data pointer
int m_row,m_col; //matrix size
};
//
// Extraction operator for matrix elements.
template<typename T>
istream& operator >>
(istream& _is, basic_matrix<T>& _mat)
{
int r = _mat.get_row();
int c = _mat.get_col();
int i,j;
for (i=0; i<r; i++)
{
cout<<"input rows "<<i+1<<": "<<endl;
for (j=0; j<r; j++)
{
_is >> _mat(i,j);
}
}
return _is;
}
// Insertion operator for matrix elements.
template<typename T,class _Elem,class _Tr>
basic_ostream<_Elem, _Tr>& operator <<
(basic_ostream<_Elem, _Tr>& _os, basic_matrix<T>& _mat)
{
if (_mat.is_empty())
return (_os << endl << "No data in this matrix" << endl );
// display all elements of basic_matrix<_Ty>
basic_ostringstream<_Elem, _Tr, allocator<_Elem> > _Sstr;
_Sstr.flags(_os.flags());
_Sstr.imbue(_os.getloc());
_Sstr.precision(_os.precision());
int r = _mat.get_row();
int c = _mat.get_col();
int i,j;
_Sstr << '[' << endl;
for (i=0; i<r; i++)
{
for (j=0; j<c; j++)
{
_Sstr << fixed << _mat(i,j) << '\t';
}_Sstr << endl;
}_Sstr << ']' <<endl;
basic_string<_Elem, _Tr, allocator<_Elem> > _Str = _Sstr.str();
return (_os << _Str.c_str());
}
template<typename T,class _Elem,class _Tr> inline
basic_ostream<_Elem, _Tr>& operator <<
(basic_ostream<_Elem, _Tr>& _os, basic_matrix<T>&& _mat)
{
return (_os << forward< basic_matrix<T>& >(_mat) );
}
//
// out/in binary files
template<typename T>
ifstream& operator >> (ifstream &_ifs, basic_matrix<T> &_mat)
{
// need the data type (T) with simple struct.
int size = _mat.get_row() * _mat.get_col() * sizeof(T);
_ifs.read((char *)&_mat(0),size);
return _ifs;
}
template<typename T>
ofstream& operator << (ofstream &_ofs, basic_matrix<T> &_mat)
{
// need the data type (T) with simple struct.
int size = _mat.get_row() * _mat.get_col() * sizeof(T);
_ofs.write((char *)&_mat(0),size);
return _ofs;
}
/************************************************************************/
/* functions */
/************************************************************************/
// convert value type
template<typename T1,typename T2>
bool convert(basic_matrix<T2> &dst,basic_matrix<T1> &src)
{
bool test = ( src.is_empty() || dst.get_row()==src.get_row() || dst.get_col()==src.get_col() );
if (test)
{
int i,j;
for (i=0; i<src.get_row(); i++)
for (j=0; j<src.get_col(); j++)
dst(i,j) = T2(src(i,j));cout<<"convertint..."<<endl;
}
return test;
}
//
template<typename T>
basic_matrix<T> inverse(basic_matrix<T> &a)
{
if (!a.is_empty() && a.get_row()==a.get_col())
{
basic_matrix<T> E(a.get_row());
E.set_ident(T(1));
return gauss_jordan(a,E);
}
else
{
cout<<"matrix must be square"<<endl;
return basic_matrix<T>();
}
}
template<typename T>
basic_matrix<T> inverse(basic_matrix<T> &&a)
{
return inverse( forward< basic_matrix<T>& > (a));
}
// Gauss eliminate
template<typename T>
basic_matrix<T> gauss_elimi(basic_matrix<T> A,basic_matrix<T> b)
{
bool test = ( !A.is_empty() && A.get_row()==A.get_col() && A.get_row()==b.get_row() );
if (test)
{
int n = A.get_row();
int col_b = b.get_col();
int i,j,k;
T C,tmp;
// P is permutation matrix
array_int P(n);
for (i=0; i<n; i++)
P(i) = i;
for (i=0; i<n; ++i)
{
// choose the pivot in columns
k = i;
for (j=i+1; j<n; j++)
if ( abs(A(P(j),i)) > abs(A(P(k),i)) )
k = j;
if ( abs(A(P(k),i))==abs(T(0)) )
return basic_matrix<T>();
if (k!=i)
swap(P(i),P(k));
for (j=i+1; j<n; j++)
{
C = A(P(j),i) / A(P(i),i);
for (k=i; k<n; k++)
A(P(j),k) -= C * A(P(i),k);
for (k=0; k<col_b; k++)
b(P(j),k) -= C * b(P(i),k);
}
}
basic_matrix<T> X(b.get_row(),b.get_col(),ONEVAL,T(0));
for (i=n-1; i>=0; --i)
{
C = T(1) / A(P(i),i);
for (j=0; j<col_b; j++)
{
tmp = T(0);
for (k=i+1; k<n; k++)
tmp += A(P(i),k) * X(k,j);
X(i,j) = C * ( b(P(i),j) - tmp );
}
}
return move(X);
}
else
{
cout<<"matrix dimensions must agree."<<endl;
return basic_matrix<T>();
}
}
// Gauss-Jordan elimination
template<typename T>
basic_matrix<T> gauss_jordan(basic_matrix<T> A,basic_matrix<T> X)
{
bool test = ( !A.is_empty() && A.get_row()==A.get_col() && A.get_row()==X.get_row() );
if (test)
{
int r = A.get_row();
int col_a = A.get_col();
int col_x = X.get_col();
array_int P(r);
int i,j,k,index;
T C;
// 下三角变换
for (i=0; i<r-1; ++i)
{
// choose the pivot in rows
index = 0;
for (k=1; k<col_a; k++)
{
if ( abs(A(i,k)) > abs(A(i,index)) )
index = k;
}
P(i) = index;
// then elimination
for (j=i+1; j<r; ++j)
{
C = A(j,index) / A(i,index);
for (k=0; k<col_a; ++k)
A(j,k) -= C * A(i,k);
for (k=0; k<col_x; ++k)
X(j,k) -= C * X(i,k);
}
}
index = 0;
for (k=1; k<col_a; k++)
{
if ( abs(A(i,k)) > abs(A(i,index)) )
index = k;
}
P(i) = index;
// up elimination
for (i=r-1; i>0; --i)
{
index = P(i);
for (j=i-1; j>=0; --j)
{
C = A(j,index) / A(i,index);
for (k=0; k<col_a; ++k)
A(j,k) -= C * A(i,k);
for (k=0; k<col_x; ++k)
X(j,k) -= C * X(i,k);
}
}
for (i=0; i<r; ++i)
{
C = T(1) / A(i,P(i));
for (j=0; j<col_x; ++j)
X(i,j) *= C;
}
// reorder
basic_matrix<T> Y(r,col_x);
for (i=0; i<r; i++)
{
for (j=0; j<col_x; j++)
{
Y(P(i),j) = X(i,j);
}
}
return move(Y);
}
else
{
cout<<"In function gauss_jordan, matrix dimensions must agree."<<endl;
return basic_matrix<T>();
}
}
// LU: Doolittle's method
template<typename T>
array_int doolittle(basic_matrix<T> &a, basic_matrix<T> &L, basic_matrix<T> &U)
{
assert(!a.is_empty() && a.get_row()==a.get_col());
int n = a.get_row();
if (U.get_row()!=n || U.get_col()!=n)
U.resize(n,n);
if (L.get_row()!=n || L.get_col()!=n)
L.resize(n,n);
U.set_one_val(T(0));
L.set_one_val(T(0));
int r,i,j,k;
T tmp;
array_int P(n); // P(n): permutation matrix
for (i=0; i<n; i++)
P(i) = i;
basic_vector<T> S(n); // S(n): maximum element in each row
for (i=0; i<n; i++)
{
tmp = a(i,0);
for (j=1; j<n; j++)
if ( abs(tmp) < abs(a(i,j)) )
tmp = a(i,j);
S(i) = tmp;
}
for (r=0; r<n; r++)
{
// choose pivot
k = r;
for (j=r+1; j<n; j++)
if ( (abs(a(P(k),r))/abs(S(k))) < (abs(a(P(j),r))/abs(S(j))) )
k = j;
if (k!=r)
swap(P(k),P(r));
// U is an upper triangular matrix
for (j=r; j<n; j++)
{
tmp = T(0);
for (k=0; k<r; k++)
tmp += L(P(r),k) * U(k,j);
U(r,j) = a(P(r),j) - tmp;
}
// L is a lower triangular matrix
L(P(r),r) = T(1);
for (i=r+1; i<n; i++)
{
tmp = T(0);
for (k=0; k<r; k++)
tmp += L(P(i),k) * U(k,r);
L(P(i),r) = ( a(P(i),r) - tmp ) / U(r,r);
}
}
/*
int r,k,i,j;
T s;
for (r=0; r<n; r++)
{
//
for (j=r; j<n; j++)
{
s = T(0);
for (k=0; k<r; k++)
s += L(r,k) * U(k,j);
U(r,j) = a(r,j) - s;
}
//
for (i=r+1; i<n; i++)
{
s = T(0);
for (k=0; k<r; k++)
s += L(i,k) * U(k,r);
L(i,r) = (a(i,r) - s) / U(r,r);
}
}
*/
return move(P);
}
//
// LU: method from Doolittle,
template<typename T>
array_int doolittle(basic_matrix<T> &a)
{
assert(!a.is_empty());
assert(a.get_row()==a.get_col());
int n = a.get_row();
array_int P(n);
int r,k,i,j;
T tmp;
for (i=0; i<n; i++)
P(i) = i;
for (r=0; r<n; r++)
{
k = r;
for (j=r+1; j<n; j++)
{
if ( abs(a(j,r)) < abs(a(k,r)))
k = j;
}
if (k!=r)
swap(P(k),P(r));
a.swap_row(k,r);
//
for (j=r; j<n; j++)
{
tmp = T(0);
for (k=0; k<r; k++)
tmp += a(r,k) * a(k,j);
a(r,j) = a(r,j) - tmp;
}
//
for (i=r+1; i<n; i++)
{
tmp = T(0);
for (k=0; k<r; k++)
tmp += a(i,k) * a(k,r);
a(i,r) = (a(i,r) - tmp) / a(r,r);
}
}
return move(P);
}
// deter of matrix
template<typename T>
T deter(basic_matrix<T> a)
{
if (!a.is_empty() && a.get_row() == a.get_col())
{
int n = a.get_row();
array_int P(n);
bool sig = false;
int i,j,k,index;
T C;
for (i=0; i<n; i++)
P(i) = i;
//
for (i=0; i<n-1; i++)
{
index = i;
for (k=i+1; k<n; k++)
if ( abs(a(i,index)) < abs(a(i,k)) )
index = k;
if (index!=i)
{
swap(P(i),P(index));
sig = !sig;
}
for (j=i+1; j<n; j++)
{
C = a(P(j),i) / a(P(i),i);
for (k=0; k<n; k++)
a(P(j),k) -= C * a(P(i),k);
}
}
C = sig?T(-1):T(1);
for (int i=0; i<n; ++i)
{
C *= a(P(i),i);
}
return C;
}
else
{
cout<<"matrix must be square!"<<endl;
return T(0);
}
}
// solve equation
template<typename T>
basic_matrix<T> solve(basic_matrix<T> &A,basic_matrix<T> &b)
{
assert(!A.is_empty());
assert(A.get_row()==A.get_col());
assert(A.get_row()==b.get_row());
int n = A.get_row();
int col_b = b.get_col();
int i,j,k;
T C,tmp;
basic_matrix<T> L(n),U(n);
basic_matrix<T> Y(n,b.get_col());
basic_matrix<T> X(n,b.get_col());
array_int P = doolittle(A,L,U);
for (i=0; i<n; i++)
{
for (j=0; j<col_b; j++)
{
tmp = T(0);
for (k=0; k<i; k++)
tmp += L(P(i),k) * Y(k,j);
Y(i,j) = b(P(i),j) - tmp ;
}
}
for (i=n-1; i>=0; --i)
{
C = U(i,i);
for (j=0; j<col_b; j++)
{
tmp = T(0);
for (k=i+1; k<n; k++)
tmp += U(i,k) * X(k,j);
X(i,j) = (Y(i,j) - tmp) / C;
}
}
return move(X);
}
// solve equation
template<typename T>
basic_matrix<T> solve2(basic_matrix<T> a,basic_matrix<T> b)
{
int n = a.get_row();
int col_b = b.get_col();
int i,j,k;
T C,tmp;
basic_matrix<T> Y(n,col_b);
array_int P = doolittle(a);
for (i=0; i<n; i++)
for (j=0; j<col_b; j++)
Y(i,j) = b(P(i),j);
for (i=0; i<n; i++)
{
for (j=0; j<col_b; j++)
{
tmp = T(0);
for (k=0; k<i; k++)
tmp += a(i,k) * b(k,j);
b(i,j) = Y(i,j) - tmp ;
}
}
for (i=n-1; i>=0; --i)
{
C = a(i,i);
for (j=0; j<col_b; j++)
{
tmp = T(0);
for (k=i+1; k<n; k++)
tmp += a(i,k) * Y(k,j);
Y(i,j) = (b(i,j) - tmp) / C;
}
}
return move(Y);
}
//
// QR: the data type used operator <,>,<=,>= and function sqrt(..)
template<typename T>
bool QR(basic_matrix<T> A,basic_matrix<T> &Q,basic_matrix<T> &R)
{
if (A.is_empty())
return false;
if (Q.is_empty())
Q.resize(A.get_row(),A.get_row());
Q.set_ident(T(1));
int a_row = A.get_row();
int a_col = A.get_col();
int N = min(a_row, a_col);
int i,j,k,t;
T sum,dr,cr,hr;
basic_vector<T> ur(a_row), pr(a_col), wr(a_row);
for (t=0; t<N; t++)
{
sum = T(0);
for (k=t; k<a_row; k++)
sum += A(k,t) * A(k,t);
dr = sqrt(sum);
cr = A(t,t)>T(0)?(-dr):dr;
hr = cr *( cr - A(t,t) );
// ur
for (i=0; i<t; i++)
ur(i) = T(0);
ur(t) = A(t,t) - cr;
for (i=t+1; i<a_row; i++)
ur(i) = A(i,t);
// wr
for (i=0; i<a_row; i++)
{
sum = T(0);
for (j=0; j<a_row; j++)
sum += Q(i,j) * ur(j);
wr(i) = sum;
}
// Q
for(i=0; i<a_row; i++)
for(j=0; j<a_row; j++)
Q(i,j) -= wr(i) * ur(j) / hr;
// pr
for(i=0; i<a_col; i++)
{
sum=0;
for(j=0; j<a_row; j++)
sum += A(j,i) * ur(j);
pr(i) = sum / hr;
}
// R
for(i=0;i<a_row;i++)
for(j=0;j<a_col;j++)
A(i,j) -= ur(i) * pr(j);
}
basic_matrix<T>::swap(A,R);
return true;
}
}
//
我测试用的cpp文件,比较乱……
#include <stdlib.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include <complex>
#include <cmath>
using namespace std;
#ifndef _CHECK_BOARDER_
#define _CHECK_BOARDER_
#endif // !_CHECK_BORDER_
#include "basic_matrix.h"
using namespace myblas;
typedef double realnum;
typedef realnum data_t;
typedef basic_matrix<data_t> matrix;
//
void rand(matrix &a)
{
srand(clock());
int i,j;
data_t C = data_t(1<<9);
for (i=0; i<a.get_row(); i++)
for (j=0; j<a.get_col(); j++)
a(i,j) = data_t(rand()&((1<<10)-1))/C;
}
realnum test(matrix &a, matrix &b)
{
if (a.get_row()!=b.get_row() && a.get_col()!=b.get_col())
return -1;
realnum esp = 0.0;
for (int i = 0; i < a.get_row(); i++)
{
for (int j = 0; j < a.get_col(); j++)
{
esp = max(esp,abs(a(i,j)-b(i,j)));
}
}
return esp;
}
realnum test(matrix &&a, matrix &&b){ return test(forward<matrix&>(a),forward<matrix&>(b)); }
realnum test(matrix &a, matrix &&b){ return test(forward<matrix&>(a),forward<matrix&>(b)); }
realnum test(matrix &&a, matrix &b){ return test(forward<matrix&>(a),forward<matrix&>(b)); }
matrix rand(int n,int m)
{
matrix a(n,m);
rand(a);
return move(a);
}
matrix rand(int n)
{
matrix a(n);
rand(a);
return move(a);
}
int main(int argc,char *argv[])
{
int m,n;
if (2==argc)
{
argv++;
m=atoi(argv[1]);
n = m;
}
else if (2<argc)
{
m=atoi(argv[1]);
n=atoi(argv[2]);
}
else
cin>>m>>n;
{
matrix a=rand(m,n),q,r;
/*
ifstream fin;
fin.open("sample",ios::binary | ios::in);
fin>>a;
fin>>b;
fin>>x0;
fin.close();
*/
/*
x = solve(a,b);
cout<<"solve 1: "<<abs((a*x-b).max_abs())<<endl;
x = solve2(a,b);
cout<<"solve 2: "<<abs((a*x-b).max_abs())<<endl;
x = gauss_elimi(a,b);
cout<<"solve 3: "<<abs((a*x-b).max_abs())<<endl;
x = gauss_jordan(a,b);
cout<<"solve 4: "<<abs((a*x-b).max_abs())<<endl;
*/
cout.precision(10);
/*
cout<<"a:"<<a<<endl;
cout<<inverse(a*a*a*a)<<endl;
cout<<(a^-4)<<endl;
*/
if (QR(a,q,r))
{
int t = 16;
cout<<"q:"<<q<<endl;
cout<<"r:"<<r<<endl;
cout<<"q*r:"<<q*r<<endl;
cout<<q*r<<endl;
for (int i=0; i<t; i++)
QR(r*q,q,r);
cout<<"r*q: "<<r*q<<endl;
}
/*
ofstream fout;
fout.open("data",ios::binary | ios::out);
fout<<a;
fout<<b;
fout<<x;
fout.close();
*/
}
// system("pause");
return 0;
}