matrix.cc

计算f t)=exp -|t|)的FFT

#include <iostream>
#include <math.h>
#include "matrix.h"
using namespace std;
void printArray(int i,Complex *array){
	for(int j=0;j<i;j++){
		array[j].print();
		cout << "\t";
	}
	cout << endl;
}

Matrix::Matrix(int d):cross(0){
	this->dimension=d;
	this->elements=new Complex*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new Complex[d];
		for(j=0;j<d;j++)
			this->elements[i][j]=0;
	}	
}

Matrix::Matrix(int d,Complex **es):cross(0){
	this->dimension=d;
	this->elements=new Complex*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new Complex[d];
		for(j=0;j<d;j++)
			this->elements[i][j]=es[i][j];
	}
}

Matrix::Matrix(int d,Complex e,int c):cross(0){ //0 for every element, 1 just for cross
	this->dimension=d;
	this->elements=new Complex*[d];
	int i,j;
	if(!c)
		for(i=0;i<d;i++){
			this->elements[i]=new Complex[d];
			for(j=0;j<d;j++)
				this->elements[i][j]=e;	
		}
	else{
		this->cross=1;
		for(i=0;i<d;i++){
			this->elements[i]=new Complex[d];
			for(j=0;j<d;j++)
				if(i==j)this->elements[i][j]=e;	
		}
	}
}

Matrix::Matrix(int d,Complex *&es,Fill_Way fw):cross(0){
	if(fw!=crs&&fw!=row&&fw!=col)d=int(sqrt(d));	
	this->dimension=d;
	this->elements=new Complex*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new Complex[d];	
		switch(fw){
			case crs:
				this->cross=1;
				for(j=0;j<d;j++){
					if(i==j)this->elements[i][j]=es[i];
					else this->elements[i][j]=0;	
				}
				break;
			case row:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[j];
				}
				break;
			case col:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i];
				}
				break;
			case all:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i*d+j];
				}
				break;
			default:
				cerr << "Warning: unknown fill way, use all as default!" << endl;
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i*d+j];
				}		
		}				
	}
}

Matrix::Matrix(int d,double **es):cross(0){
	this->dimension=d;
	this->elements=new Complex*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new Complex[d];
		for(j=0;j<d;j++)
			this->elements[i][j]=es[i][j];
	}
}

Matrix::Matrix(int d,double e,int c):cross(0){ //0 for every element, 1 just for cross
	this->dimension=d;
	this->elements=new Complex*[d];
	int i,j;
	if(!c)
		for(i=0;i<d;i++){
			this->elements[i]=new Complex[d];
			for(j=0;j<d;j++)
				this->elements[i][j]=e;	
		}
	else{
		this->cross=1;
		for(i=0;i<d;i++){
			this->elements[i]=new Complex[d];
			for(j=0;j<d;j++)
				if(i==j)this->elements[i][j]=e;	
		}
	}
}

Matrix::Matrix(int d,double *es,Fill_Way fw):cross(0){
	if(fw!=crs&&fw!=row&&fw!=col)d=int(sqrt(d));	
	this->dimension=d;
	this->elements=new Complex*[d];
	int i,j;
	for(i=0;i<d;i++){
		this->elements[i]=new Complex[d];	
		switch(fw){
			case crs:
				this->cross=1;
				for(j=0;j<d;j++){
					if(i==j)this->elements[i][j]=es[i];
					else this->elements[i][j]=0;	
				}
				break;
			case row:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[j];
				}
				break;
			case col:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i];
				}
				break;
			case all:
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i*d+j];
				}
				break;
			default:
				cerr << "Warning: unknown fill way, use all as default!" << endl;
				for(j=0;j<d;j++){
					this->elements[i][j]=es[i*d+j];
				}		
		}				
	}
}

Matrix::Matrix(const Matrix &m1,const Matrix &m2,const Matrix &m3,const Matrix &m4):cross(0){	
	if(m1.dimension-m2.dimension||m1.dimension-m3.dimension||m2.dimension-m3.dimension||
	   m1.dimension-m4.dimension||m2.dimension-m4.dimension||m3.dimension-m4.dimension){
		cerr <<	"Warning: dimension of 4 matrix are not equal, the this matrix!" << endl;
		exit(1);
	}
	int d=m1.dimension;
	this->dimension=2*d;
	this->elements=new Complex*[2*d];
	int i,j;
	for(i=0;i<2*d;i++){
		this->elements[i]=new Complex[2*d];		
		for(j=0;j<2*d;j++){
			if(i<d){
				if(j<d)this->elements[i][j]=m1.elements[i][j];
				else this->elements[i][j]=m2.elements[i][j-d];
			}
			else{
				if(j<d)this->elements[i][j]=m3.elements[i-d][j];
				else{ this->elements[i][j]=m4.elements[i-d][j-d];		
				}
			}
		}
	}//this->elements[0][0].print();
	//this->setCross();	
}

int Matrix::getDimension(void){
	return this->dimension;
}

Complex Matrix::getElement(int i,int j){
	if(i>=this->dimension||j>=this->dimension){
		cerr << "Warning: subscript out of boundary, return the first element!" << endl;
		return this->elements[0][0];
	}
	return this->elements[i][j];
}

void Matrix::setElement(Complex e,int i,int j){
	if(i>=this->dimension||j>=this->dimension){
		cerr << "Warning: subscript out of boundary, no elements set!" << endl;
		return;
	}	
	this->elements[i][j]=e;
	this->setCross();
}

void Matrix::setElement(double e,int i,int j){
	if(i>=this->dimension||j>=this->dimension){
		cerr << "Warning: subscript out of boundary, no elements set!" << endl;
		return;
	}	
	this->elements[i][j]=e;
	this->setCross();
}

int Matrix::operator==(Matrix &m){
	if(this->dimension-m.dimension)return 0;
	int i,j;
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			if(this->elements[i][j]-m.getElement(i,j)==Complex(0,0))return 0;
		}
	}
	return 1;
}

void Matrix::setCross(void){
	int i,j;
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			if(i-j){
				if(!(this->elements[i][j]==Complex(0,0))){
					this->cross=0;return;
				};
			}
		}
	}
	this->cross=1;
}

int Matrix::getCross(void){
	return this->cross;
}

Complex *Matrix::getRow(int i){
	if(i>=this->dimension){
		cerr << "Warning: subscript out of boundary, return the first row!" << endl;
		return this->elements[0];
	}
	return this->elements[i];
}

Complex *Matrix::getCol(int j){
	Complex *tempCol;
	int i;
	tempCol=new Complex[this->dimension];
	if(j>=this->dimension){
		cerr << "Warning: subscript out of boundary, return the first row!" << endl;
		return this->getCol(0);
	}
	for(i=0;i<this->dimension;i++){
		tempCol[i]=this->elements[i][j];
	}
	return tempCol;
}

void Matrix::setRow(Complex *es,int i){
	int j;
	if(i>=this->dimension){
		cerr << "Warning: subscript out of boundary, no rows set!" << endl;
		return;
	}
	for(j=0;j<this->dimension;j++)
		this->setElement(es[j],i,j);
	this->setCross();
}

void Matrix::setCol(Complex *es,int j){
	int i;
	if(j>=this->dimension){
		cerr << "Warning: subscript out of boundary, no colomns set!" << endl;
		return;
	}	
	for(i=0;i<this->dimension;i++)
		this->setElement(es[i],i,j);
	this->setCross();
}

Complex *Matrix::getCrossE(void){
	Complex *tempCross;
	int i,j;
	tempCross=new Complex[this->dimension];
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			if(i==j){
				tempCross[i]=this->elements[i][j];
				break;
			}
		}
	}
	return tempCross;
}

Matrix Matrix::trans(void){
	int i,j;
	Matrix *r=new Matrix(this->dimension);
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			r->setElement(this->elements[i][j],j,i);
		}
	}
	return *r;
}

Complex *Matrix::operator*(double *s){
	int i,j;
	Complex *temp=new Complex[this->dimension];
	for(i=0;i<this->dimension;i++){
		temp[i]=0;
		for(j=0;j<this->dimension;j++){
			temp[i]=temp[i]+(this->elements[i][j]*s[j]);
		}
	}
	return temp;
}

Complex *Matrix::operator*(Complex *s){
	int i,j;
	Complex *temp=new Complex[this->dimension];
	for(i=0;i<this->dimension;i++){
		temp[i]=0;
		for(j=0;j<this->dimension;j++){
			temp[i]=temp[i]+(this->elements[i][j]*s[j]);
		}
	}
	return temp;
}

Matrix Matrix::operator*(const Matrix &m){
	if(this->dimension-m.dimension){
		cerr << "Warning: dimensions of two matrix not equal, return the left matrix!" << endl;
		return *this;
	}
	int i,j;
	Matrix *r=new Matrix(this->dimension); 
	Complex *es,*cs;
	es=new Complex[this->dimension];
	cs=new Complex[this->dimension];
	if(this->cross){					//左乘乘行 
		cs=this->getCrossE();
		for(i=0;i<this->dimension;i++){
			//es=m.getRow(i);
			for(j=0;j<this->dimension;j++){
				es[j]=cs[i]*m.elements[i][j];
			}
			r->setRow(es,i);
		}
	}
	else if(m.cross){					//右乘乘列 
		//cs=m.getCrossE();
		for(i=0;i<this->dimension;i++){
			es=this->getCol(i);
			for(j=0;j<this->dimension;j++){
				if(i==j)es[j]=m.elements[i][j]*es[j];
			}
			r->setCol(es,i);
		}
	}
	else{								//直接乘法 
		int s;
		for(s=0;s<this->dimension;s++){
			for(i=0;i<this->dimension;i++){
				for(j=0;j<this->dimension;j++){
					r->setElement(r->getElement(i,j)+this->elements[i][s]*m.elements[s][j],i,j);
				}		
			}			
		}
	}
	r->setCross();
	return *r;
}

Matrix Matrix::operator-(Matrix &m){
	Matrix *r=new Matrix(this->dimension);
	int i,j;
	for(i=0;i<this->dimension;i++){
		for(j=0;j<this->dimension;j++){
			r->setElement(this->elements[i][j]-m.getElement(i,j),i,j);
		}
	}
	r->setCross();
	return *r;
}

void Matrix::setCrossE(Complex *es){
	int i,j;	
	for(i=0;i<this->dimension;i++)
		for(j=0;j<this->dimension;j++)
			if(i==j)this->setElement(es[i],i,j);
	this->setCross();
}

void Matrix::print(int d){
	int i,j;
	if(d>this->dimension)this->print();
	else
		for(i=0;i<d;i++){
			for(j=0;j<d;j++){
				(this->elements[i][j]).print();
				cout << "\t";
			}
			cout << endl;
		}
}

void Matrix::print(void){
	this->print(this->dimension);
}

Matrix::~Matrix(void){
	for(int i=0;i<this->dimension;i++)
		delete [] this->elements[i];
	delete [] this->elements;		
}