mymat.h

电力系统潮流计算 采取PQ分解法 用VC++完成编程。

#ifndef MYMAT_H
#define MYMAT_H

#include<iostream>
#include<iomanip>
#include<new>
using namespace std;


template <class T>
class MyMat{
	friend ostream & operator<<(ostream &,MyMat<T> &);
	
private:
	int nn;
	int mm;
	T**v;

public:
	MyMat();
	MyMat(int n,int m);
	MyMat(const T& a,int n,int m);
	MyMat(const T* a,int n,int m);
	MyMat(const MyMat& rhs);
    MyMat& operator=(const MyMat &rhs);
    MyMat& operator=(const T & a);
	inline T*operator[](const int i);
	inline const T* operator[](const int i )const;
	inline int Get_Row() const;
	inline int Get_Col() const;
	MyMat conj();
	MyMat inv();

	bool operator==(const MyMat&rhs);
	bool operator!=(const MyMat&rhs);
	MyMat operator+(const MyMat&rhs);
	MyMat operator-(const MyMat&rhs);
	MyMat operator*(const MyMat&rhs);
    MyMat operator+=(const MyMat&rhs);
	MyMat operator-=(const MyMat&rhs);


	~MyMat();

};

template <class T>
MyMat<T>::MyMat():nn(0),mm(0),v(0){}

template <class T>
MyMat<T>::MyMat(int n,int m):nn(n),mm(m),v(new T*[n]){
	v[0]=new T[m*n];
	for(int i=1;i<n;i++)
		v[i]=v[i-1]+m;
}

template <class T>
MyMat<T>::MyMat(const T& a,int n,int m):nn(n),mm(m),v(new T*[n]){
	int i,j;
	v[0]=new T[m*n];
	for(i=1;i<n;i++)
		v[i]=v[i-1]+m;
	for(i=0;i<n;i++)
		for(j=0;j<m;j++)
			v[i][j]=a;
}

template <class T>
MyMat<T>::MyMat(const T* a,int n,int m):nn(n),mm(m),v(new T*[n]){
	int i,j;
	v[0]=new T[m*n];
	for(i=1;i<n;i++)
		v[i]=v[i-1]+m;
	for(i=0;i<n;i++)
		for(j=0;j<m;j++)
			v[i][j]=*a++;
}

template <class T>
MyMat<T>::MyMat(const MyMat& rhs):nn(rhs.nn),mm(rhs.mm),v(new T*[nn]){
	int i,j;
	v[0]=new T[mm*nn];
	for(i=1;i<nn;i++)
		v[i]=v[i-1]+mm;
	for(i=0;i<nn;i++)
		for(j=0;j<mm;j++)
			v[i][j]=rhs[i][j];
}

template <class T>
MyMat<T>& MyMat<T>::operator=(const MyMat<T> &rhs){
	if(this!=&rhs){
		int i,j;
		if(nn!=rhs.nn||mm!=rhs.mm){
			if(v!=0){
				delete[](v[0]);
				delete[](v);
			}
			nn=rhs.nn;
			mm=rhs.mm;
			v=new T*[nn];
			v[0]=new T[mm*nn];
		}
		for(i=1;i<nn;i++)
			v[i]=v[i-1]+mm;
		for(i=0;i<nn;i++)
			for(j=0;j<mm;j++)
				v[i][j]=rhs[i][j];
	}
	return *this;
}

template <class T>
MyMat<T>& MyMat<T>::operator =(const T& a){
	for(int i=0;i<nn;i++)
		for(int j=0;j<mm;j++)
			v[i][j]=a;
		return *this;
}

template <class T>
inline T* MyMat<T>::operator[](const int i){
	return v[i];
}

template <class T>
inline const T* MyMat<T>::operator[](const int i)const{
	return v[i];
}

template <class T>
inline int MyMat<T>::Get_Row()const{
	return nn;
}

template <class T>
inline int MyMat<T>::Get_Col()const{
	return mm;
}

template <class T>
MyMat<T> MyMat<T>::conj(){
	MyMat<T> output(this->Get_Col(),this->Get_Row());
	int i,j;
	for(i=0;i<this->Get_Col();i++)
		for(j=0;j<this->Get_Row();j++)
			output[i][j]=(*this)[j][i];
		return output;
}

template <class T>
MyMat<T> MyMat<T>::inv(){
	int i,j;
	MyMat<T> output(this->Get_Row(),this->Get_Row());
	if(this->Get_Row()!=this->Get_Col())
	{
		cout<<"不是方阵，不能求逆！";
		exit(1);
	}
	else{
			T* A=new T[(this->Get_Row())*(this->Get_Row())];
			for(int ii=0;ii<this->Get_Row();ii++)
				for(int jj=0;jj<this->Get_Row();jj++)
					A[ii*(this->Get_Row())+jj]=(*this)[ii][jj];
		for(j=0;j<this->Get_Row();j++){
			T* B=new T[this->Get_Row()];
			for(int p=0;p<this->Get_Row();p++)
				B[p]=0;
			B[j]=1;
			for(i=0;i<this->Get_Row();i++){
				output[i][j]=G_S(A,B,this->Get_Row(),i);

			}
		}
	}
	return output;
}


template <class T>
bool MyMat<T>::operator==(const MyMat&rhs){
	int i,j;
	if(this->nn!=rhs.nn||this->mm!=rhs.mm)
		return false;
	else{
		for( i=0;i<nn;i++)
		
			for( j=0;j<mm;j++){
			if(this->v[i][j]!=rhs.v[i][j])
				break;
		}
		if(i>nn-1&&j>mm-1)
			return true;
		else
			return false;

	}
}


template <class T>
bool MyMat<T>::operator!=(const MyMat&rhs){
	int i,j;
	if(this->nn!=rhs.nn||this->mm!=rhs.mm)
		return true;
	else{
		for( i=0;i<nn;i++)
			for( j=0;j<mm;j++){
			if(this->v[i][j]!=rhs.v[i][j])
				break;
		}
		if(i>nn-1&&j>mm-1)
			return false;
		else
			return true;

	}
}

template <class T>
ostream&  operator<<(ostream &output, MyMat<T>&Mat){
	
    for(int i=0;i<Mat.Get_Row();i++){
		for(int j=0;j<Mat.Get_Col();j++)
			output<<left<<setw(20)<<Mat[i][j];
		output<<endl;
	}
	return output;
}


template <class T>
MyMat<T> MyMat<T>::operator+(const MyMat<T>&rhs){
	MyMat<T> output(rhs.Get_Row(),rhs.Get_Col());
	if(this->Get_Row()!=rhs.Get_Row()||this->Get_Col()!=rhs.Get_Col()){
		cout<<"矩阵维数不等，无法计算！";
		exit(1);
	}
	else
		if(this->Get_Row()==0||this->Get_Col()==0)
			cout<<"空矩阵！";
		else{
			for(int i=0;i<this-> Get_Row();i++)
				for(int j=0;j<this->Get_Col();j++)
				output[i][j]=(*this)[i][j]+rhs[i][j];
			
		}
		return output;
}


template <class T>
MyMat<T> MyMat<T>::operator-(const MyMat<T>&rhs){
	MyMat<T> output(rhs.Get_Row(),rhs.Get_Col());
	if(this->Get_Row()!=rhs.Get_Row()||this->Get_Col()!=rhs.Get_Col()){
		cout<<"矩阵维数不等，无法计算！";
		exit(1);
	}
	else
		if(this->Get_Row()==0||this->Get_Col()==0)
			cout<<"空矩阵！";
		else{
			for(int i=0;i<this-> Get_Row();i++)
				for(int j=0;j<this->Get_Col();j++)
				output[i][j]=(*this)[i][j]-rhs[i][j];
			
		}
		return output;
}

template <class T>
MyMat<T> MyMat<T>::operator*(const MyMat<T>&rhs){
	int i,j,k;
	MyMat<T> output(this->Get_Row(),rhs.Get_Col());
	
	if(this->Get_Col()!=rhs.Get_Row()){
		cout<<"无法相乘！";
		exit(1);
	}
	else{
		
		for(i=0;i<Get_Row();i++)
			for(j=0;j<rhs.Get_Col();j++)
				output[i][j]=0;
        for(i=0;i<Get_Row();i++)
			for(j=0;j<rhs.Get_Col();j++){
		for(int k=0;k<rhs.Get_Row();k++)
			output[i][j]+=(*this)[i][k]*rhs[k][j];
			}
	}
	return output;

}

template <class T>
MyMat<T> MyMat<T>::operator+=(const MyMat<T>&rhs){
	
	if(this->Get_Row()!=rhs.Get_Row()||this->Get_Col()!=rhs.Get_Col()){
		cout<<"矩阵维数不等，无法计算！";
		exit(1);
	}
	else
		if(this->Get_Row()==0||this->Get_Col()==0)
			cout<<"空向量！";
		else{
			for(int i=0;i<this->Get_Row();i++)
				for(int j=0;j<this->Get_Col();j++)
				(*this)[i][j]+=rhs[i][j];
			
		}
		return *this;
}

template <class T>
MyMat<T> MyMat<T>::operator-=(const MyMat<T>&rhs){
	
	if(this->Get_Row()!=rhs.Get_Row()||this->Get_Col()!=rhs.Get_Col()){
		cout<<"矩阵维数不等，无法计算！";
		exit(1);
	}
	else
		if(this->Get_Row()==0||this->Get_Col()==0)
			cout<<"空向量！";
		else{
			for(int i=0;i<this->Get_Row();i++)
				for(int j=0;j<this->Get_Col();j++)
				(*this)[i][j]-=rhs[i][j];
			
		}
		return *this;
}

template <class T>
T G_S(T *ptra,T *ptrb,int di,int n){
	int i,j,k,ii;
	T max;
	T*PA=new T[di*di];
	T *PB=new T[di];
	for(i=0;i<di;i++)
		for(j=0;j<di;j++)
			PA[i*di+j]=ptra[i*di+j];
	for(i=0;i<di;i++)
		PB[i]=ptrb[i];
	T*tempA=new T[di];
	T tempB;
	for(i=0;i<di;i++){
		max=fabs(PA[i*di+i]);
		for( ii=i;ii<di;ii++){
			if(max<fabs(PA[ii*di+i]))
				max=fabs(PA[ii*di+i]);
		}
		for(ii=i;max!=fabs(PA[ii*di+i]);ii++)
			;
		for(j=0;j<di;j++){
			tempA[j]=PA[ii*di+j];
			PA[ii*di+j]=PA[i*di+j];
			PA[i*di+j]=tempA[j];
		}
		tempB=PB[ii];
		PB[ii]=PB[i];
		PB[i]=tempB;
		T temp,scal;
		for(k=i+1;k<di;k++){
		temp=PA[k*di+i];
		scal=temp/PA[i*di+i];
		for(j=0;j<di;j++)
			PA[k*di+j]-=scal*PA[i*di+j];
		PB[k]-=scal*PB[i];
		}
	}
	T *X=new T[di];
	X[di-1]=PB[di-1]/PA[(di-1)*di+di-1];
	for(i=di-2;i>=0;i--){
		X[i]=PB[i];
		for(j=di-1;j>i;j--)
			X[i]-=PA[i*di+j]*X[j];
		X[i]/=PA[i*di+i];
	}
	return X[n];
}


template <class T>
MyMat<T>::~MyMat(){
	if(v!=0){
		delete[](v[0]);
		delete[] v;
	}
}
#endif