矩阵运算.txt

矩阵运算

MATRIX.H 
#include "iostream.h" 
class MATRIX 
{ 
    double **data;//矩阵数据 
    int n,m;//n行m列 
public: 
    MATRIX(); 
    ~MATRIX(); 
    void input(); 
    void output(); 
    MATRIX(MATRIX&); 
    MATRIX(int,int); 
    MATRIX operator+(const MATRIX&); 
    MATRIX operator-(const MATRIX&); 
    MATRIX operator*(const MATRIX&); 
    MATRIX operator*(double); 
    MATRIX operator/(double); 
    friend MATRIX operator*(double,MATRIX); 
    friend double det(MATRIX);//行列式值 
    friend MATRIX yzs(MATRIX,int,int);//代数余子式 
    friend MATRIX inv(MATRIX);//矩阵逆 
    double get(int,int);//取第i行第j列的值 
    friend double Vmax(MATRIX M1);//求向量最大分量 
    friend double tzz(MATRIX M1);//幂法求最大特征值 
    void put(double,int,int);//将一个值输入到第i行j列 
    void LU(MATRIX,MATRIX);//LU分解放入传入的两个矩阵中 
    MATRIX(double**,int,int);//将一个2维数组转化为矩阵 
}; 
double abs(double);//绝对值 
MATRIX.CPP 
#include "MATRIX.h" 
MATRIX::MATRIX() 
{ 
    n=0;m=0; 
} 
MATRIX::MATRIX(int a,int b) 
{ 
    n=a; 
    m=b; 
    data=new (double*[n+1]); 
    int i,j; 
    for(i=1;i<=n;i++) 
        data[i]=new (double[m+1]); 
    for(i=1;i<=n;i++) 
        for(j=1;j<=m;j++) 
            data[i][j]=0; 
} 
MATRIX::MATRIX(MATRIX& M) 
{ 
    this->m=M.m; 
    this->n=M.n; 
    data=new (double*[n+1]); 
    int i,j; 
    for(i=1;i<=n;i++) 
        data[i]=new (double[m+1]); 
    for(i=1;i<=this->n;i++) 
        for(j=1;j<=this->m;j++) 
            this->data[i][j]=M.data[i][j]; 
} 
MATRIX::~MATRIX() 
{ 
    int i; 
    if(n!=0&&m!=0) 
    { 
        for(i=1;i<=n;i++) 
            delete data[i]; 
        delete data; 
    } 
} 
MATRIX MATRIX::operator +(const MATRIX& M1) 
{ 
    if(M1.n==this->n&&M1.m==this->m) 
    { 
        MATRIX ret(this->n,this->m); 
        for(int i=1;i<=M1.n;i++) 
            for(int j=1;j<=M1.m;j++) 
                ret.data[i][j]=M1.data[i][j]+this->data[i][j]; 
        return ret; 
    } 
    else 
    { 
        cout<<"error"<<endl; 
        return MATRIX(); 
    } 
} 
MATRIX MATRIX::operator -(const MATRIX& M1) 
{ 
    if(this->n==M1.n&&this->m==M1.m) 
    { 
        MATRIX ret(this->n,this->m); 
        for(int i=1;i<=M1.n;i++) 
            for(int j=1;j<=M1.m;j++) 
                ret.data[i][j]=this->data[i][j]-M1.data[i][j]; 
        return ret; 
    } 
    else 
    { 
        cout<<"error"<<endl; 
        return MATRIX(); 
    } 
} 
MATRIX MATRIX::operator *(const MATRIX& M1) 
{ 
    if(this->n==M1.m&&this->m==M1.n) 
    { 
        MATRIX ret(this->n,M1.m); 
        double sum; 
        for(int i=1;i<=this->n;i++) 
            for(int j=1;j<=M1.m;j++) 
            { 
                sum=0; 
                for(int k=1;k<=this->m;k++) 
                    sum+=this->data[i][k]*M1.data[k][j]; 
                ret.data[i][j]=sum; 
            } 
        return ret; 
    } 
    else 
    { 
        cout<<"error"<<endl; 
        return MATRIX(); 
    } 
} 
MATRIX MATRIX::operator /(double d) 
{ 
    if(d!=0) 
    { 
        MATRIX ret(this->n,this->m); 
        for(int i=1;i<=this->n;i++) 
            for(int j=1;j<=this->m;j++) 
                ret.data[i][j]=this->data[i][j]/d; 
        return ret; 
    } 
    else 
    { 
        cout<<"error"<<endl; 
        return MATRIX(); 
    } 
} 
MATRIX MATRIX::operator *(double d) 
{ 
    if(d!=0) 
    { 
        MATRIX ret(this->n,this->m); 
        for(int i=1;i<=this->n;i++) 
            for(int j=1;j<=this->m;j++) 
                ret.data[i][j]=this->data[i][j]*d; 
        return ret; 
    } 
    else 
    { 
        cout<<"error"<<endl; 
        return MATRIX(); 
    } 
} 
double det(MATRIX M1) 
{ 
    if(M1.n==1&&M1.m==1) 
        return M1.data[1][1]; 
    if(M1.n==M1.m) 
    { 
        int sign=1; 
        double sum=0; 
        for(int i=1;i<=M1.n;i++) 
        { 
            MATRIX M2(M1.n-1,M1.m-1); 
            for(int j=1;j<=M2.n;j++) 
                for(int k=1;k<=M2.m;k++) 
                { 
                    if(k<i) 
                        M2.data[j][k]=M1.data[j+1][k]; 
                    else 
                        M2.data[j][k]=M1.data[j+1][k+1]; 
                } 
            sum+=sign*M1.data[1][i]*det(M2); 
            sign=sign*(-1); 
        } 
        return sum; 
    } 
    else 
    { 
        cout<<"error"<<endl; 
        return 0; 
    } 
} 
void MATRIX::input() 
{ 
    for(int i=1;i<=this->n;i++) 
        for(int j=1;j<=this->m;j++) 
        { 
            cout<<i<<"行"<<j<<"列"<<endl; 
            cin>>this->data[i][j]; 
        } 
} 
void MATRIX::output() 
{ 
    for(int i=1;i<=this->n;i++) 
    { 
        for(int j=1;j<=this->m;j++) 
        { 
            cout<<this->data[i][j]<<"  "; 
        } 
        cout<<endl; 
    } 
} 
MATRIX yzs(MATRIX M1,int a,int b) 
{ 
    if(a<=M1.n&&b<=M1.m) 
    { 
        MATRIX M2(M1.n-1,M1.m-1); 
        for(int i=1;i<=M2.n;i++) 
            for(int j=1;j<=M2.m;j++) 
            { 
                if(i>=a) 
                { 
                    if(j>=b) 
                    { 
                        M2.data[i][j]=M1.data[i+1][j+1]; 
                    } 
                    else 
                    { 
                        M2.data[i][j]=M1.data[i+1][j]; 
                    } 
                } 
                else 
                { 
                    if(j>=b) 
                    { 
                        M2.data[i][j]=M1.data[i][j+1]; 
                    } 
                    else 
                    { 
                        M2.data[i][j]=M1.data[i][j]; 
                    } 
                } 
            } 
        return M2; 
    } 
    else 
    { 
        return MATRIX(); 
    } 
} 
MATRIX inv(MATRIX M1) 
{ 
    if(M1.n!=M1.m) 
    { 
        return MATRIX(); 
    } 
    MATRIX M2(M1.n,M1.m); 
    for(int i=1;i<=M2.n;i++) 
        for(int j=1;j<=M2.m;j++) 
            M2.data[i][j]=det(yzs(M1,i,j)); 
    return (M2/det(M1)); 
} 
double MATRIX::get(int i,int j) 
{ 
    return this->data[i][j]; 
} 
double Vmax(MATRIX M1) 
{ 
    if(M1.m!=1) 
        return 0; 
    double max=M1.data[1][1]; 
    for(int i=1;i<=M1.n;i++) 
        if(M1.data[i][1]>max) 
            max=M1.data[i][1]; 
    return max; 
} 
double tzz(MATRIX M1) 
{ 
    if(M1.n!=M1.m) 
        return 0; 
    double t1,t2; 
    MATRIX m_nVector1(M1.n,1),m_nVector2(M1.n,1); 
    for(int i=1;i<=M1.n;i++) 
        m_nVector1.data[i][1]=M1.data[i][1]; 
    m_nVector2=M1*m_nVector1; 
    t1=Vmax(m_nVector2); 
    m_nVector1=m_nVector2/t1; 
    for(i=0;i<100000;i++) 
    { 
        m_nVector2=M1*m_nVector1; 
        t2=Vmax(m_nVector2); 
        if(abs(t1-t2)<0.00005) 
            return t1; 
        m_nVector1=m_nVector2/t2; 
        t1=t2; 
    } 
    return 0; 
} 
double abs(double t1) 
{ 
    if(t1<0) 
        return (0-t1); 
    else 
        return t1; 
} 
MATRIX operator*(double a,MATRIX M1) 
{ 
    return M1*a; 
} 
void MATRIX::LU(MATRIX L,MATRIX U) 
{ 
    if(n!=m) 
    { 
        cout<<"error"<<endl; 
        return; 
    } 
    double **l,**u; 
    l=new (double*[n+1]); 
    u=new (double*[n+1]); 
    int i,j,k; 
    for(i=1;i<=n;i++) 
    { 
        l[i]=new (double[n+1]); 
        u[i]=new (double[n+1]); 
    } 
    for(i=1;i<=n;i++) 
        for(j=1;j<=n;j++) 
        { 
            l[i][j]=0; 
            u[i][j]=0; 
        } 
    for(i=1;i<=n;i++) 
        l[i][i]=1; 
    for(i=1;i<=n;i++) 
        u[1][i]=data[1][i]; 
    for(i=2;i<=n;i++) 
        l[i][1]=data[i][1]/u[1][1]; 
    double sum=0; 
    for(int r=2;r<=n;r++) 
    { 
        for(i=r;i<=n;i++) 
        { 
            sum=0; 
            for(k=1;k<=r-1;k++) 
                sum+=l[r][k]*u[k][i]; 
            u[r][i]=data[r][i]-sum; 
            sum=0; 
            for(k=1;k<=r-1;k++) 
                sum+=l[i][k]*u[k][r]; 
            l[i][r]=(data[i][r]-sum)/u[r][r]; 
        } 
    } 
    MATRIX MT1(n,n),MT2(n,n); 
    for(i=1;i<=n;i++) 
        for(j=1;j<=n;j++) 
        { 
            MT1.put(l[i][j],i,j); 
            MT2.put(u[i][j],i,j); 
        } 
    L=MT1; 
    U=MT2; 
    return; 
} 
void MATRIX::put(double c,int a,int b) 
{ 
    this->data[a][b]=c; 
    return; 
} 
MATRIX::MATRIX(double** c,int a,int b) 
{ 
    this->n=a; 
    this->m=b; 
    data=new (double*[n+1]); 
    int i,j; 
    for(i=1;i<=n;i++) 
        data[i]=new (double[m+1]); 
    for(i=1;i<=n;i++) 
        for(j=1;j<=m;j++) 
            data[i][j]=c[i][j];