matrix.h

自己写的矩阵向量库

#include <iostream>
#include <stdlib.h>
#include <algorithm>
#include <math.h>
#include "Vtr.h" 

class Mtx
{
private:
	int nrows;
	int ncols;
	double** ets;
public:
	Mtx(int n,int m, double**);
	Mtx(int n, int m, double d=0);
	Mtx(const Mtx&);
	~Mtx();

	Mtx& operator=(const Mtx&);
	Mtx& operator+=(const Mtx&);
	Mtx& operator-=(const Mtx&);
	Vtr operator*(const Vtr&) const;
	double* operator[](int i) const {return ets[i];}
	double& operator()(int i, int j) {return ets[i][j];}

	friend Mtx operator+(const Mtx&) ;
	Mtx operator+(const Mtx&) const;
	Mtx operator* (const Mtx&) const;
	friend Mtx operator-(const Mtx&);
	friend Mtx operator-(const Mtx&, const Mtx&);

	int CG(Vtr& x, const Vtr& b, double& eps, int& iter);
};


Mtx::Mtx(int n, int m, double** dbp)
{
	nrows = n;
	ncols = m;
	ets = new double* [nrows];
	for(int i=0; i<nrows; i++)
	{
		ets[i]=new double [ncols];
		for(int j=0; j<ncols; j++)  ets[i][j]=dbp[i][j];
	}
}

Mtx::Mtx(int n,int m,double a)
{
	ets=new double* [nrows=n];
	ncols=m;
	for(int i=0; i<nrows; i++)
	{
		ets[i]=new double [ncols];
		for(int j=0; j<ncols; j++)  ets[i][j]=a;
	}
}

Mtx::Mtx(const Mtx& mat)
{
	ets=new double* [nrows=mat.nrows];
	ncols=mat.ncols;
	for(int i=0; i<nrows; i++)
	{
		ets[i]=new double [ncols];
		for(int j=0; j<ncols; j++) ets[i][j]=mat[i][j];
	}
}

inline Mtx::~Mtx()
{
	for(int i=0; i<nrows; i++) delete[] ets[i];
	delete[] ets;
}

Mtx& Mtx::operator =(const Mtx& mat)
{
	if(this != &mat)
	{
		if(nrows != mat.nrows || ncols != mat.ncols)
			error("bad matrix sizes");
		for(int i=0; i<nrows;i++)
			for(int j=0; j<ncols; j++)
				ets[i][j]=mat[i][j];
	}
	return *this;
}

Mtx& Mtx::operator +=(const Mtx& mat)
{
    if(nrows != mat.nrows || ncols != mat.ncols)
		error("bad matrix sizes");
    for(int i=0; i<nrows;i++)
		for(int j=0; j<ncols; j++)
			ets[i][j]+=mat[i][j];
	return *this;
}

Mtx& Mtx::operator -=(const Mtx& mat)
{
    if(nrows != mat.nrows || ncols != mat.ncols)
		error("bad matrix sizes");
    for(int i=0; i<nrows;i++)
		for(int j=0; j<ncols; j++)
			ets[i][j] -= mat[i][j];
	return *this;
}

inline Mtx operator +(const Mtx& mat) {return mat;}

inline Mtx operator-(const Mtx& mat)
{
	return Mtx(mat.nrows, mat.ncols)-mat;
}

Mtx Mtx::operator +(const Mtx& mat) const
{
    if(nrows != mat.nrows || ncols != mat.ncols)
		error("bad matrix sizes");
	Mtx sum = *this;
	sum += mat;
    return sum;
}

Mtx operator -(const Mtx& m1, const Mtx& m2) 
{
    if(m1.nrows != m2.nrows || m1.ncols != m2.ncols)
		error("bad matrix sizes");
	Mtx sum = m1;
	sum -= m2;
    return sum;
}

Vtr Mtx::operator *(const Vtr& v) const
{
	if(ncols != v.size())
		error("matrix and vector sizes do not match");
	Vtr tm(nrows);
	for(int i=0; i<nrows; i++)
		for(int j=0; j<ncols; j++)
			tm[i] += ets[i][j]*v[j];
	return tm;
}

Mtx Mtx::operator*(const Mtx& b) const
{
	if(ncols != b.nrows)
		error("matrix sizes do not match");
	Mtx m(nrows, b.ncols);
	for(int i=0; i<nrows; i++)
		for(int j=0; j<b.ncols; j++)
			for(int k=0; k<ncols; k++)
				m[i][j]=m[i][j]+ets[i][k]*b[k][j];
			return m;
}

int Mtx::CG(Vtr& x, const Vtr& b, double&eps, int&iter)
{
	if(nrows != b.size())
		cout << "matrix and vecter sizes do not match\n";
	const int maxiter=iter;
	Vtr r=b-(*this)*x;
	Vtr p=r;
	double zr=dot(r,r);
	const double stp=eps*b.twonorm();

	if (r.twonorm()==0) {
		eps=0;
		iter=0;
		return 0;
	}

	for(iter=0; iter<maxiter; iter++)
	{
		Vtr mp=(*this)*p;
		double alpha=zr/dot(p,p);
		x+=alpha*p;
		r-=alpha*mp;
		if (r.twonorm()<=stp) break;
		double zrold=zr;
		zr=dot(r,r);
		p=r+(zr/zrold)*p;

	}

	eps=r.twonorm();
	if(iter==maxiter) return 1;
	else return 0;

}
/*
int main()
{
	int n=300;
	Mtx a(n,n);
	for(int i=0; i<n; i++)
		for(int j=0; j<n; j++) a[i][j]=1/(i+j+1.0);
    
    Vtr t(n);
	Vtr x(n);
	for(i=0; i<n; i++)
		t[i]=1/(i+3.14);

	int iter=300;
	double eps=1.0e-9;

	int ret=a.CG(x,a*t, eps, iter);
    Vtr y(n);
	y = x-t;
	if(ret==0) 
		cout <<"CG returned successfully\n";
	cout <<iter << "iterations are used in CG methed. \n";
	cout << "Residual="<<eps <<".\n";
	cout << "Tow-norm of exact vecter = "
		 << y.twonorm() << '\n';
	
    return 0;
}
*/