// Matrix.cpp : 实现文件
//

#include "stdafx.h"
#include "math.h"
#include "Matrix.h"
#include ".\matrix.h"


#include <assert.h>


// Matrix

Matrix::Matrix()
{
    this->p = new double[1];
    this->Col = 1;
    this->Row = 1;
}

Matrix::~Matrix()
{
    delete []p;
}


// Matrix 成员函数

long Matrix::GetRows(void)
{
    return this->Row;
}

long Matrix::GetCols(void)
{
    return this->Col;
}

Matrix::Matrix(long n)
{
    this->Row = this->Col = n;
    this->p = new double[n * n];
    for(long i = 0 ; i < n ; i++)
        for(long j = 0 ;j < n; j++)
            if(i == j) *(p + n * i + j) = 1;
            else *(p + n * i + j) = 0;
}
Matrix::Matrix(long row ,long col)
{
    this->Row = row;
    this->Col = col;
    this->p    = new double[row * col];
    for(long i = 0 ; i < row ; i++)
        for(long j = 0 ; j < col ; j ++)
            *(p +col * i + j) = 0;

}
Matrix::Matrix(double *arrAddress ,long arrWidth)
{
    long arrHeight = 1;
    this->p = new double[arrWidth * arrHeight];
    for(long i = 0 ; i < arrHeight ; i++)
        for(long j = 0 ; j < arrWidth; j++)
            *(p + arrWidth * i + j ) = *(arrAddress + arrWidth * i + j);
    this->Col = arrWidth ;
    this->Row = arrHeight;
}
Matrix::Matrix(double *arrAddress,long rows,long cols)
{
    this->p = new double[rows * cols];
    for(long i = 0 ; i < rows ; i++)
        for(long j = 0 ; j < cols ; j ++)
            *(p +cols * i + j) = *(arrAddress + cols * i + j );
    this->Col = cols;
    this->Row = rows;
}
Matrix::Matrix(const Matrix &m)//拷贝构造函数
{
    this->Col = m.Col;
    this->Row = m.Row;
    p = new double[this->Col * this->Row];
    for(long  i = 0 ; i < this->Row ; i++)
        for(long j = 0 ; j < this->Col ; j ++)
            *(p + this->Col * i + j) = *(m.p + this->Col * i + j);
}
bool Matrix::IsVector(void)
{
    return !(this->Col == 1 && this->Row ==1);
}

Matrix Matrix::SubMatrix(long offset)
{
    assert(this->Col ==    this->Row && offset <= this->Col && offset >= 0);
    double *t = new double[offset * offset] ;
    for(long i = 0 ; i < offset ; i++)
        for(long j = 0 ; j < offset ; j++)
            *(t + offset * i + j )  = *(p + this->Col * i + j);
    Matrix m(t ,offset ,offset);
    delete []t;
    return m;
}

double Matrix::Arg(void)
{
    assert(this->Row == this->Col);
    double result = 1;
    double k ;
    Matrix m = *this;
    for(long i = 0 ; i < this->Row - 1 ; i++)
    {
        for(long j = i + 1; j < this->Row ; j++)
        {
            k = m[j][i] / m[i][i];
            m[j][i] = 0 ;
            for(long n = i + 1; n < this->Col ; n ++)
            {
                m[j][n] = m[j][n] - k * m[i][n];
            }
        }
    }
    for(long i = 0 ; i < this->Row ; i++)
    {
        for(long j = 0 ; j < this->Col ; j++)
        {
            if(i == j ) result *= m[i][j];
        }
    }
    return result;
}



double * Matrix::operator[](long heightPos)
{
    assert(heightPos >= 0 && heightPos < this->Row);
    return this->p + this->Col * heightPos;
    return NULL;
}

bool Matrix::IsPtv(void)
{
    assert(this->Col == this->Row);//是方阵才能计算
    bool result = true;
    Matrix m;
    for(long i = 1 ; i<= this->Row; i++)
    {
        m = this->SubMatrix(i);
        if(m.Arg() <= 0)
        {
            result = false;
            break;
        }
    }
    return result;
}

Matrix Matrix::T(void)
{
    double *t = new double[this->Col * this->Row] ;
    for(long i = 0 ; i< this->Row ;i++)
        for(long j = 0 ; j < this->Col ; j++)
            *(t + this->Row * j + i) = *(this->p + this->Col * i + j);
    Matrix m(t ,this->Col ,this->Row);
    delete []t;
    return m;
}
/*
Matrix Matrix::operator *(Matrix & m1)
{
    if(!this->IsVector() && m1.IsVector())//左为数右为矩阵
    {
        Matrix m;
        m = this->p[0] * m1;
        return m;
    }
    else if(this->IsVector() && !m1.IsVector())
    {
        Matrix m;
        m = *this * m1[0][0];
        return m;
    }
    else if(!this->IsVector() && m1.IsVector())
    {
        double *t = new double[1] ;
        t[0] = p[0] * m1[0][0];
        Matrix m(t,1,1);
        delete []t;
        return m;
    }
    else if(this->IsVector() && m1.IsVector() && this->Col == m1.Row)
    {
        double sum;
        double *t = new double[this->Row * m1.Col] ;
        for(long i = 0 ; i < this->Row ; i++)
        {
            for(long j = 0 ; j < m1.Col ; j++)
            {
                sum = 0 ;
                for(long k = 0 ; k < this->Col ; k++)
                {
                    sum += (*(p + this->Col * i + k)) * (m1[k][j]);
                }
                *(t + m1.Col * i + j) = sum ;
            }
        }
        Matrix m(t ,m1.Col ,m1.Row);
        delete []t;
        return m;
    }
    else
    {
        assert(0) ; //未知运算
        return *this;
    }
}
*/
/*
Matrix Matrix::operator *(double alpha ,Matrix &m1)
{
    Matrix m = m1;
    for(long i = 0 ; i < m.Row ; i++)
    {
        for(long j = 0 ;j < m.Col; j++)
        {
            m[i][j] = alpha * m1[i][j];
        }
    }
    return m;
}

*/
Matrix Matrix::operator *(double alpha)
{
    for(int i = 0 ; i < this->Row ; i++)
        for(int j = 0 ; j < this->Col ; j++)
            *(this->p + this->Row * j + i) = *(this->p + this->Row * j + i) * alpha;

    return (*this);
}

Matrix Matrix::operator /(Matrix &m1)
{
    Matrix t = this->T();
    t = t * *(this);
    t = t.Inver();
    t = t * this->T();
    t = t * m1;
    return t;
}
Matrix Matrix::operator *(Matrix &m1)
{
    if(this->Col == m1.Row)
    {
        Matrix ttt(*this);
        int mr = this->Row;
        int mc = m1.Col;
        Matrix tt(mr ,mc);
        for(long i = 0 ; i < mr ; i++)
        {
            for(long j = 0 ; j < mc; j++)
            {
                for(int ii = 0 ; ii < this->Col; ii++)
                {
                    
                    tt[i][j] += ttt[i][ii] * m1[ii][j];
                }
            }
        }
        return tt;
    }
    /*
    if(this->IsVector() && m1.IsVector() && this->Col == m1.Row)
    {
        double sum;
        double *t = new double[this->Row * m1.Col] ;
        for(long i = 0 ; i < this->Row ; i++)
        {
            for(long j = 0 ; j < m1.Col ; j++)
            {
                sum = 0 ;
                for(long k = 0 ; k < this->Col ; k++)
                {
                    sum += (*(p + this->Col * i + k)) * (m1[k][j]);
                }
                *(t + m1.Col * i + j) = sum ;
            }
        }
        Matrix m(t ,m1.Col ,m1.Row);
        delete []t;
        return m;
    }
    else
    {
        assert(0) ; //未知运算
        return *this;
    }*/

}

Matrix Matrix::Inver()
{
    Matrix srcm(*this);
    int rhc = srcm.Row;
    if(srcm.Row == srcm.Col)
    {
        long *iss = new long[rhc];
        long *jss = new long[rhc];
        double fdet = 1;
        double f =1;
        for(long k = 0 ; k < rhc ; k++)
        {
            double fmax = 0;
            for(long i = k ; i < rhc ; i++)
            {
                for(long j = k ; j < rhc; j ++)
                {
                    f = ::abs(srcm[i][j]);
                    if(f > fmax)
                    {
                        fmax = f;
                        *(iss + k) = i;
                        *(jss + k) = j;
                    }
                }
            }
            if(*(iss + k) != k)
            {
                f = -f;
                for(long ii = 0 ; ii < rhc ; ii++)
                {
                    double tt = srcm[k][ii];
                    srcm[k][ii] = srcm[*(iss + k) ][ii];
                    srcm[*(iss + k)][ii] = tt;
                }
            }
            if(*(jss + k) != k)
            {
                f = -f;
                for(long ii = 0 ; ii < rhc ; ii++)
                {
                    double tt = srcm[k][ii];
                    srcm[k][ii] = srcm[*(jss + k) ][ii];
                    srcm[*(jss + k)][ii] = tt;
                }
            }
            fdet *= srcm[k][k];
            srcm[k][k] = 1 / srcm[k][k];
            for(long j = 0 ; j < rhc; j ++)
                if(j != k )
                    srcm[k][j] *= srcm[k][k];
            for(long i = 0; i < rhc ; i++)
                if(i != k)
                    for(long j = 0 ; j < rhc; j++)
                        if(j != k)
                            srcm[i][j] = srcm[i][j] - srcm[i][k] * srcm[k][j];
            for(long i = 0 ; i < rhc ; i++)
                if( i != k)srcm[i][k] *= -srcm[k][k];
        }
        for(long k = rhc -1 ; k>= 0 ; k--)
        {
            if(*(jss + k) != k)
                for(long ii = 0 ; ii < rhc; ii++)
                {
                    double tt = srcm[k][ii] ;
                    srcm[k][ii] = srcm[*(jss + k) ][ii];
                    srcm[*(jss + k) ][ii] = tt;
                }
            if(*(iss + k ) != k)
                for(long ii = 0 ; ii < rhc ; ii++)
                {
                    double tt = srcm[k][ii];
                    srcm[k][ii] = srcm[*(iss + k)][ii];
                    srcm[*(iss + k)][ii] = tt;
                }
        }
        delete []iss;
        delete []jss;
    }
    
    return srcm;

}
Matrix Matrix::operator =(Matrix &m)
{
    if(&m == this)return *this;
    this->Row = m.Row;
    this->Col = m.Col;
    delete []this->p;
    p = new double[this->Row * this->Col];
    for(long i = 0 ; i < this->Row ; i++)
    {
        for(long j = 0 ; j < this->Col ; j++)
        {
            *(this->p + this->Col * i + j) = *(m.p + m.Col * i + j);
        }
    }
    return *this;
}

void Matrix::Print(void)
{

    for(long i = 0 ; i < this->Row ; i ++)
    {
        for(long j = 0 ; j < this->Col ; j++)
        {
            printf("%f  " ,*(this->p + this->Col * i + j));
        }
        printf("\n");
    }

}

请问,这个程序是不是只能求方阵的逆阵啊???

求广义逆么？