matrix.h

这是一个求解线性规划问题的程序

//---------------------------------------------------------------------------
#ifndef MatrixH
#define MatrixH
//---------------------------------------------------------------------------

#include "TValue.h"

#define TME_OUTOFMEMMORY      0
#define TME_INVALIDOPERATION  1
#define TME_INVALIDPARAM      2
#define TME_OUTOFBOUND        3

//矩阵异常类
class TMException
{
public:
    unsigned int mErrorCode;
public:
    TMException(){
    }
    TMException(unsigned int aE){
        mErrorCode=aE;
    }
    operator unsigned int(void){    //重载(int)强制转换:如 if(e) dosomething;
        return mErrorCode;
    }
};

typedef TValue TMVALUE;


//矩阵类
class TMatrix
{
protected:
   //用一维数组存矩阵,访问元素时可不用A[3][2],从而避免了乘法的地址运算:A+3*M+2
   //加快了矩阵运算的速度
   //前缀"m"--成员变量member
    TMVALUE   *  mO;      //首地址
    unsigned int mRowCount,mColCount;    //行mRowCount,列mColCount

private:
    //逆矩阵
    TMatrix Invert(void)
    {
        if (mRowCount!=mColCount)
           throw TMException(TME_INVALIDOPERATION);

        TMatrix lTemp = *this;
        TMatrix lMatR(mRowCount,mColCount);

        unsigned int i,j,k;
        //初始化一个单位矩阵
        for(i=1;i<=mRowCount;i++) {
           for(j=1;j<=mColCount;j++){
              if (i==j)
                lMatR(i,j) = 1;
              else
                lMatR(i,j) = 0;
           }
        }

        //把lTemp变成一个单位矩阵,同时把相应的行变换作用于lMatR;
        //下三角
        for(j=1,i=j;j<=mColCount;j++,i++){
            //找到该列中第一个不为0的元素
            for (k=i;k<=mRowCount;k++)
                if(lTemp(k,j)!=TValue(0))
                    break;
            if (k>mRowCount){ //一列全为0
                i--;
                continue;
            }
            if (i!=k) {
                lTemp.RExchange(i,k);
                lMatR.RExchange(i,k);
            }

            //进行初等变换,每行加上一行的n倍
            for(k=i+1;k<=mRowCount;k++) {
                lMatR.RAddition(k,i,-(lTemp(k,j)/lTemp(i,j)));
                lTemp.RAddition(k,i,-(lTemp(k,j)/lTemp(i,j)));
            }
        }

        //如果不满秩
        if (lTemp(mRowCount,mColCount)==0)
           throw TMException(TME_INVALIDOPERATION);

        //上三角
        for(i=mRowCount;i>=1;i--){
            //进行初等变换,每行加上一行的n倍
            for(k=i-1;k>=1;k--) {
                lMatR.RAddition(k,i,-(lTemp(k,i)/lTemp(i,i)));
                lTemp.RAddition(k,i,-(lTemp(k,i)/lTemp(i,i)));
            }
        }

        //把对角阵化为单位阵
        for(i=1;i<=mRowCount;i++) {
           lMatR.RScale(i,1/lTemp(i,i));
           lTemp.RScale(i,1/lTemp(i,i));
        }

        return lMatR;

    }

public:
    //构造函数,构造一个aRowCount*aColCount的矩阵
    TMatrix(unsigned int aRowCount,unsigned int aColCount)
    {
        if( aRowCount==0 || aColCount==0)
            throw TMException(TME_INVALIDPARAM);

        mRowCount=aRowCount;
        mColCount=aColCount;
        mO=0;

        mO=new TMVALUE[mRowCount*mColCount];
        if(!mO)
            throw TMException(TME_OUTOFMEMMORY);

        unsigned int lMu=mRowCount*mColCount;
        for(unsigned int i=0;i<lMu;i++)
            mO[i] = 0;
    }
    //拷贝构造函数
    TMatrix(TMatrix& aAnother)
    {
        mRowCount=aAnother.mRowCount;
        mColCount=aAnother.mColCount;
        mO=new TMVALUE[mRowCount*mColCount];
        if(!mO)
            throw TMException(TME_OUTOFMEMMORY);

        unsigned int lMu=mRowCount*mColCount;
        for(unsigned int i=0;i<lMu;i++)
            mO[i]=aAnother.mO[i];
    }
    TMatrix(TMatrix* aAnother)
    {
        mRowCount=aAnother->mRowCount;
        mColCount=aAnother->mColCount;
        mO=new TMVALUE[mRowCount*mColCount];
        if(!mO)
            throw TMException(TME_OUTOFMEMMORY);

        unsigned int lMu=mRowCount*mColCount;
        for(unsigned int i=0;i<lMu;i++)
            mO[i]=aAnother->mO[i];
    }
    
    //析构函数,释放内存空间
    ~TMatrix()
    {
       if (mO)
         delete [] mO;
    }

    //重载=号操作符
    TMatrix operator=(TMatrix aAnother)
    {
        if( (mRowCount!= aAnother.mRowCount) || (mColCount!=aAnother.mColCount) )
            throw TMException(TME_INVALIDOPERATION);

        unsigned int lMu=mRowCount*mColCount;
        for(unsigned int i=0;i<lMu;i++)
            mO[i]=aAnother.mO[i];
        return *this;
    }

    //重置维数
    void Reset(unsigned int aRowCount,unsigned int aColCount)
    {
        if( aRowCount==0 || aColCount==0)
            throw TMException(TME_INVALIDPARAM);

        delete [] mO;
        
        mRowCount=aRowCount;
        mColCount=aColCount;
        mO=0;

        mO=new TMVALUE[mRowCount*mColCount];
        if(!mO)
            throw TMException(TME_OUTOFMEMMORY);

        unsigned int lMu=mRowCount*mColCount;
        for(unsigned int i=0;i<lMu;i++)
            mO[i] = 0;
    }

    //重载()操作符,取出第aI行第aJ列的元素,
    //注意返回的是该元素的引用,所以可以用来赋值:A(1,2)=1;
    TMVALUE& operator()(unsigned int aI,unsigned int aJ)
    {
        if( (aI>mRowCount) || (aJ>mColCount) || (aI==0) || (aJ==0) )
            throw TMException(TME_OUTOFBOUND);
        return *(mO+(aI-1)*mColCount+aJ-1);
    }

    unsigned int GetRowCount(void) { return mRowCount; }
    unsigned int GetColCount(void) { return mColCount; }

    //代数系统的重载操作符.
    //关系
    bool operator ==(TMatrix& aAnother)
    {
        if( (mRowCount!=aAnother.mRowCount) || (mColCount!=aAnother.mColCount) )
            throw TMException(TME_INVALIDOPERATION);

        unsigned int lMu=mRowCount*mColCount;
        for(unsigned int i=0;i<lMu;i++)
            if(mO[i]!=aAnother.mO[i]) return false;
        return true;
    }
    bool operator !=(TMatrix& aAnother)
    {
        return !(*this==aAnother);
    }

    //运算
    //加法
    TMatrix operator+(TMatrix& aAnother)
    {
        if( (mRowCount!=aAnother.mRowCount) || (mColCount!=aAnother.mColCount) )
            throw TMException(TME_INVALIDOPERATION);
        TMatrix lMatResult=*this;

        unsigned int lMu=lMatResult.mRowCount * lMatResult.mColCount;
        for(unsigned int i=0;i<lMu;i++)
            lMatResult.mO[i]=mO[i]+aAnother.mO[i];
        return lMatResult;
    }

    //减法
    TMatrix operator-(TMatrix& aAnother)
    {
        if( (mRowCount!=aAnother.mRowCount) || (mColCount!=aAnother.mColCount) )
            throw TMException(TME_INVALIDOPERATION);
        TMatrix lMatResult=*this;

        unsigned int lMu=lMatResult.mRowCount * lMatResult.mColCount;
        for(unsigned int i=0;i<lMu;i++)
            lMatResult.mO[i]=mO[i]-aAnother.mO[i];
        return lMatResult;
    }

    //乘法
    TMatrix operator*(TMatrix& aAnother)
    {
        if( mColCount!=aAnother.mRowCount )
            throw TMException(TME_INVALIDOPERATION);

        TMatrix lMatR(mRowCount,aAnother.mColCount);

        unsigned int lMu=lMatR.mColCount*lMatR.mRowCount;
        for(TMVALUE *lRow=lMatR.mO,*lThisRow=mO;
            lRow<lMatR.mO+lMu                   ;
            lRow+=lMatR.mColCount,lThisRow+=mColCount         ){


            for(TMVALUE *lWhere=lRow,*lAnoCol=aAnother.mO;
                lWhere<lRow+lMatR.mColCount                      ;
                lWhere++,lAnoCol++                        ){


                *lWhere=0;
                for(TMVALUE *lThisWhere=lThisRow ,*lAnoWhere=lAnoCol;
                    lThisWhere<lThisRow+mColCount                          ;
                    lThisWhere++,lAnoWhere+=aAnother.mColCount             ){
                    *lWhere = (*lThisWhere) * (*lAnoWhere) + *lWhere;
                }
            }
        }
        return lMatR;
    }

    //幂运算
    TMatrix operator^(int aIndex){
        if (mRowCount!=mColCount)
           throw TMException(TME_INVALIDOPERATION);

        TMatrix lMatR(mRowCount,mColCount);

        lMatR = *this;

        for(int i=aIndex;i>1;i--) {
            lMatR = lMatR * *this;
        }
        for(int i=aIndex;i<-1;i++) {
            lMatR = lMatR * *this;
        }

        //逆矩阵
        if (aIndex == -1)
           lMatR = lMatR.Invert();

        //单位阵
        if (aIndex == 0) {
           for (unsigned int i=1;i<=mRowCount;i++) {
               for(unsigned int j=1;j<=mColCount;j++) {
                   if (i==j)
                      lMatR(i,j) = 1;
                   else
                      lMatR(i,j) = 0;
               }
           }
        }
        return lMatR;

    }
    //转置
    TMatrix operator~(void){

        TMatrix lMatR(mColCount,mRowCount);

        unsigned int lMu=mColCount*mRowCount;
        for(TMVALUE* lRow=lMatR.mO,*lThisCol=mO;
            lRow<lMatR.mO+lMu                  ;
            lRow+=lMatR.mColCount,lThisCol++          ){

            for(TMVALUE* lWhere=lRow,*lThisWhere=lThisCol;
                lWhere<lRow+lMatR.mColCount                      ;
                lWhere++,lThisWhere+=mColCount                   ){

                *lWhere=*lThisWhere;
            }
        }
        return lMatR;
    }

    //矩阵的秩(即行向量组的秩)
    int Rank(void)
    {
        TMatrix lTemp=*this;

        //先化成阶梯矩阵 (秩 = 行数 - 全为0的行数) 
        lTemp.Ladder();

        int RZeros=0;
        bool RZero;
        for (unsigned int i=1;i<=mRowCount;i++) {
            RZero=true;
            for (unsigned int j=1;j<=mColCount && RZero;j++) {
                if (lTemp(i,j)!=0)
                   RZero = false;
            }
            if (RZero)
              RZeros+=1;
        }

        return mRowCount-RZeros;
    }

    //阶梯矩阵
    void Ladder(void)
    {
        for(unsigned int j=1,i=j;j<=mColCount;j++,i++){
            //找到该列中第一个不为0的元素
            unsigned int k;
            for (k=i;k<=mRowCount;k++)
                if((*this)(k,j)!=0)
                    break;
            if (k>mRowCount){ //一列全为0
                i--;
                continue;
            }
            if (i!=k)
                RExchange(i,k);

            //进行初等变换,每行减去一行的n倍
            for(k=i+1;k<=mRowCount;k++) {
                RAddition(k,i,-((*this)(k,j)/(*this)(i,j)));
            }
        }
    }

    //取某列的向量
    TMatrix CVector(unsigned int aJ)
    {
        if (aJ>mColCount || aJ<1)
           throw TMException(TME_OUTOFBOUND);

        TMatrix lMatR(mRowCount,1);
        for (unsigned int i=1;i<=mRowCount;i++) {
            lMatR(i,1) = (*this)(i,aJ);
        }

        return lMatR;
    }
    //取某行的向量
    TMatrix RVector(unsigned int aI)
    {
        if (aI>mRowCount || aI<1)
           throw TMException(TME_OUTOFBOUND);

        TMatrix lMatR(1,mColCount);
        for (unsigned int i=1;i<=mColCount;i++) {
            lMatR(1,i) = (*this)(aI,i);
        }

        return lMatR;
    }
    //初等变换
    //行初等变换
    //交换两行
    void RExchange(unsigned int aI1,unsigned int aI2)
    {
        TMVALUE *lRow1=&(*this)(aI1,1);
        TMVALUE *lRow2=&(*this)(aI2,1);
        for(TMVALUE * lWhere1=lRow1,*lWhere2=lRow2;
            lWhere1<lRow1+mColCount                       ;
            lWhere1++,lWhere2++                    ){
            TMVALUE lTemp=*lWhere1;
            *lWhere1=*lWhere2;
            *lWhere2=lTemp;
        }
    }
    //一行乘与aS倍
    void RScale(unsigned int aI,TMVALUE aS)
    {
        TMVALUE *lRow=&(*this)(aI,1);
        for(TMVALUE * lWhere=lRow;
            lWhere<lRow+mColCount;
            lWhere++){
            *lWhere=aS*(*lWhere);
        }
    }
    //一行加上另一行的aS倍
    void RAddition(unsigned int aI1,unsigned int aI2,TMVALUE aS)
    {
        TMVALUE *lRow1=&(*this)(aI1,1);
        TMVALUE *lRow2=&(*this)(aI2,1);
        for(TMVALUE * lWhere1=lRow1,*lWhere2=lRow2;
            lWhere1<lRow1+mColCount                       ;
            lWhere1++,lWhere2++                    ){
            *lWhere1=(*lWhere1)+(*lWhere2)*aS;
        }
    }

    //列初等变换
    void CExchange(unsigned int aJ1,unsigned int aJ2)
    {
        TMVALUE *lCol1=&(*this)(1,aJ1);
        TMVALUE *lCol2=&(*this)(1,aJ2);

        unsigned int lMu=mRowCount*mColCount;
        for(TMVALUE * lWhere1=lCol1,*lWhere2=lCol2;
            lWhere1<lCol1+lMu                      ;
            lWhere1+=mColCount,lWhere2+=mColCount                ){
            TMVALUE lTemp=*lWhere1;
            *lWhere1=*lWhere2;
            *lWhere2=lTemp;
        }
    }
    void CScale(unsigned int aJ,TMVALUE aS)
    {
        TMVALUE *lCol=&(*this)(1,aJ);

        unsigned int lMu=mRowCount*mColCount;
        for(TMVALUE * lWhere=lCol;
            lWhere<lCol+lMu;
            lWhere+=mColCount){
            *lWhere=(*lWhere)*aS;
        }
    }
    void CAddition(unsigned int aJ1,unsigned int aJ2,TMVALUE aS)
    {
        TMVALUE *lCol1=&(*this)(1,aJ1);
        TMVALUE *lCol2=&(*this)(1,aJ2);

        unsigned int lMu=mRowCount*mColCount;
        for(TMVALUE * lWhere1=lCol1,*lWhere2=lCol2;
            lWhere1<lCol1+lMu                      ;
            lWhere1+=mColCount,lWhere2+=mColCount                ){
            *lWhere1=(*lWhere1)+(*lWhere2)*aS;
        }
    }
};

//---------------------------------------------------------------------------

#endif