c++版矩阵基本操作，行列式，逆（不限矩阵大小）

来源：互联网发布：手机必备软件大全编辑：程序博客网时间：2024/05/21 17:07
原本是为了编程实现线性回归的，想想，里面太多矩阵操作，尤其是求逆。以前学数值分析时，也用到过列主元高斯消去求解线性方程组，LU分解求解线性方程组。这次，同样是用高斯消去法求矩阵行列式的值，用LU分解求解矩阵的逆，效率上程序执行起来还行，比用python跑一边速度快，结果一致，这也潜在说明python库中矩阵求逆的实现应该也是用的LU分解。至于矩阵的其他一些操作，基本上算简单，当然面的稀疏性矩阵的话，采用三元组的形式表示，运算起来会更好，但这里不考虑，可以放到数据结构数组的表示方式那一章中。下面给出c++实现的代码
#include <iostream>#include <stdlib.h>#include <string>#include <math.h>#include "loadData.h"#include <fstream>#include <sstream>#include <stack>using namespace std;#define MAX_SIZE_OF_TRAINING_SET 100#define MAX_NUMIT 100#define ATTR_NUM 3#define MAX 1000000#define MIN -100000#define MAX_MATRIX_COL 1000#define MAX_MATRIX_ROW 100class Matrix{public:    double **mat;    int col,row;public:    int loadMatrix(Matrix *matrix,dataToMatrix dtm)    {        int i,j;        Data *p;        p=dtm.dataSet->next;        matrix->mat=(double **)malloc(sizeof(double*)*dtm.col);        for(i=0; i<dtm.col&&p!=NULL; i++)        {            matrix->mat[i]=(double *)malloc(sizeof(double)*dtm.row);            for(j=0; j<dtm.row; j++)            {                matrix->mat[i][j]=p->attr_double[j];            }            p=p->next;        }        matrix->row=dtm.row;        matrix->col=dtm.col;        return 0;    }    int initMatrix(Matrix *matrix,int col,int row)    {        matrix->col=col;        matrix->row=row;        matrix->mat=(double **)malloc(sizeof(double*)*col);        int i=0,j=0;        for(i=0; i<col; i++)        {            matrix->mat[i]=(double *)malloc(sizeof(double)*row);            for(j=0; j<row; j++)                matrix->mat[i][j]=0;        }        return 0;    }    int initMatrix(Matrix *matrix,int col,int row,double lam)    {        matrix->col=col;        matrix->row=row;        matrix->mat=(double **)malloc(sizeof(double*)*col);        int i=0,j=0;        for(i=0; i<col; i++)        {            matrix->mat[i]=(double *)malloc(sizeof(double)*row);            for(j=0; j<row; j++)            {                matrix->mat[i][j]=0;                if(i==j)                matrix->mat[i][j]=lam;            }        }        return 0;    }    int print(Matrix matrix)    {        int i,j;        for(i=0; i<matrix.col; i++)        {            for(j=0; j<matrix.row; j++)            {                cout<<matrix.mat[i][j]<<"  ";            }            cout<<endl;        }    }    int copy(Matrix matrixA,Matrix *matrixB)    {        int i,j;        //matrixB->mat=(double **)malloc(sizeof(double*)*matrixA.col);        for(i=0; i<matrixA.col; i++)        {            //matrixB->mat[i]=(double *)malloc(sizeof(double)*matrixA.row);            for(j=0; j<matrixA.row; j++)            {                matrixB->mat[i][j]=matrixA.mat[i][j];            }        }        matrixB->col=matrixA.col;        matrixB->row=matrixA.row;        return 0;    }    Matrix getOneRow(Matrix matrix,int iRow)    {        Matrix oneRow;        oneRow.col=matrix.col;        oneRow.row=1;        int i=0;        initMatrix(&oneRow,oneRow.col,oneRow.row);        for(i=0; i<oneRow.col; i++)        {            oneRow.mat[i][0]=matrix.mat[i][iRow-1];        }        return oneRow;    }    Matrix getOneCol(Matrix matrix,int iCol)    {        Matrix oneCol;        oneCol.row=matrix.row;        oneCol.col=1;        int i=0;        initMatrix(&oneCol,oneCol.col,oneCol.row);        for(i=0; i<oneCol.row; i++)        {            oneCol.mat[0][i]=matrix.mat[iCol][i];        }        return oneCol;    }    int deleteOneRow(Matrix *matrix,int iRow)    {        int i,j;        for(i=iRow; i<matrix->col; i++)        {            //for()//由于传递来的一般是最后一列，所以只需要列数--即可，不需移动，不写了        }        matrix->row--;    }    void transposematrix(Matrix matrix,Matrix *matrixT)//矩阵形式的转置    {        int i=0,j=0;        matrixT->col=matrix.row;        matrixT->row=matrix.col;        //matrixT->mat=(double **)malloc(sizeof(double *)*matrixT->col);        for(i=0; i<matrixT->col; i++)        {            //matrixT->mat[i]=(double *)malloc(sizeof(double)*matrixT->row);            for(j=0; j<matrixT->row; j++)            {                matrixT->mat[i][j]=matrix.mat[j][i];                //cout<<matrixT->mat[i][j]<<"  ";            }            //cout<<endl;        }    }    int addmatrix(Matrix *addMatrix,Matrix matrix1,Matrix matrix2)    {        if(matrix1.col!=matrix2.col||matrix1.row!=matrix2.row)            return -1;        int i,j;        addMatrix->col=matrix1.col;        addMatrix->row=matrix1.row;        //addMatrix->mat=(double **)malloc(sizeof(double *)*addMatrix->col);        for(i=0; i<matrix1.col; i++)        {            //addMatrix->mat[i]=(double *)malloc(sizeof(double)*addMatrix->row);            for(j=0; j<matrix1.row; j++)            {                addMatrix->mat[i][j]=matrix1.mat[i][j]+matrix2.mat[i][j];                //cout<<addMatrix->mat[i][j]<<"  ";            }            //cout<<endl;        }        return 0;    }    int submatrix(Matrix *addMatrix,Matrix matrix1,Matrix matrix2)    {        if(matrix1.col!=matrix2.col||matrix1.row!=matrix2.row)            return -1;        int i,j;        addMatrix->col=matrix1.col;        addMatrix->row=matrix1.row;        //addMatrix->mat=(double **)malloc(sizeof(double *)*addMatrix->col);        for(i=0; i<matrix1.col; i++)        {            //addMatrix->mat[i]=(double *)malloc(sizeof(double)*addMatrix->row);            for(j=0; j<matrix1.row; j++)            {                addMatrix->mat[i][j]=matrix1.mat[i][j]-matrix2.mat[i][j];                //cout<<addMatrix->mat[i][j]<<"  ";            }            //cout<<endl;        }        return 0;    }    int multsmatrix(Matrix *multsMatrix,Matrix matrix1,Matrix matrix2)//矩阵形式的相乘    {        if(matrix1.row!=matrix2.col)            return -1;        int i,j,k,l;        multsMatrix->col=matrix1.col;        multsMatrix->row=matrix2.row;        //multsMatrix->mat=(double **)malloc(sizeof(double *)*(multsMatrix->col));        for(i=0; i<matrix1.col; i++)        {           // multsMatrix->mat[i]=(double *)malloc(sizeof(double)*(multsMatrix->row));            for(j=0; j<matrix2.row; j++)            {                multsMatrix->mat[i][j]=0;            }        }        for(i=0; i<matrix1.col; i++)        {            for(j=0; j<matrix2.row; j++)            {                for(k=0; k<matrix1.row; k++)                {                    multsMatrix->mat[i][j]+=matrix1.mat[i][k]*matrix2.mat[k][j];                }            }        }        return 0;    }    //行列式    double detmatrix(Matrix matrix)    {        if(matrix.col!=matrix.row)            return -1;        double det=1;        int i=0,j,k;        double max=MIN;        int swap=-1;        double temp;        double aij[MAX_MATRIX_COL][MAX_MATRIX_ROW];        for(k=0; k<matrix.row-1; k++)//k表示第k次消元，一共需要n-1次        {            for(i=0; i<matrix.col; i++)            {                if(matrix.mat[i][k]>max)//每一次消元都是比较第k列的元素，选出第k列中最大的一行                {                    swap=i;                }            }//找到第k次列主元消去的最大行的下标            if(swap==-1||matrix.mat[swap][k]==0)                return -1;//最大主元为0            for(j=0; j<matrix.row; j++)            {                temp=matrix.mat[k][j];                matrix.mat[k][j]=matrix.mat[swap][j];                matrix.mat[swap][j]=temp;            }//第k次消元，选出最大的一行是swap行，与第k行交换            for(i=k+1; i<matrix.col; i++)            {                aij[i][k]=matrix.mat[i][k]/matrix.mat[k][k];// 第k次消元，主元素为第k行第k列，把第k行以下的行都进行消元                for(j=k; j<matrix.row; j++)//对于k行以下的每一行的每一列元素都减去主行与消元因子的乘积                {                    matrix.mat[i][j]-=aij[i][k]*matrix.mat[k][j];                }            }        }        for(i=0; i<matrix.col; i++)        {            det*=matrix.mat[i][i];            /*            for(j=0; j<matrix.row; j++)            {                cout<<aij[i][j]<<"  ";            }            cout<<endl;            */        }        cout<<"det="<<det<<endl;        return det;    }    //高斯消元矩阵求逆,特别注意，LU分解不能进行行列式变换    int nimatrix(Matrix *niMatrix,Matrix matrix)    {        if(matrix.col!=matrix.row)            return -1;        //if(detmatrix(matrix)==0)//这里调用求行列式进行了列主元消去改变了参数矩阵，如何传递不改变是一个问题            //return -1;        int i=0,j,k;        double temp;        Matrix cpMatrix;        Matrix uMatrix;        Matrix lMatrix;        Matrix uniMatrix;        Matrix lniMatrix;        initMatrix(&uniMatrix,matrix.col,matrix.row);        initMatrix(&lniMatrix,matrix.col,matrix.row);        initMatrix(&cpMatrix,matrix.col,matrix.row);        initMatrix(&uMatrix,matrix.col,matrix.row);        initMatrix(&lMatrix,uMatrix.col,uMatrix.row);        copy(matrix,&cpMatrix);        //cout<<"cpMatrix"<<endl;        //print(cpMatrix);        double aij[MAX_MATRIX_COL][MAX_MATRIX_ROW];        for(k=0; k<matrix.row-1; k++)//k表示第k次消元，一共需要n-1次        {            for(i=k+1; i<matrix.col; i++)            {                aij[i][k]=matrix.mat[i][k]/matrix.mat[k][k];// 第k次消元，主元素为第k行第k列，把第k行以下的行都进行消元                for(j=k; j<matrix.row; j++)//对于k行以下的每一行的每一列元素都减去主行与消元因子的乘积                {                    matrix.mat[i][j]-=aij[i][k]*matrix.mat[k][j];                }            }        }        copy(matrix,&uMatrix);        cout<<"uMatrix"<<endl;        print(uMatrix);        for(j=0; j<matrix.row; j++)        {            for(i=j+1; i<matrix.col; i++)            {                temp=0;                for(k=0; k<j; k++)                {                    temp=lMatrix.mat[i][k]*uMatrix.mat[k][j];                }                lMatrix.mat[i][j]=1/uMatrix.mat[j][j]*(cpMatrix.mat[i][j]-temp);            }        }        for(i=0; i<lMatrix.col; i++)        {            for(j=0; j<lMatrix.row; j++)            {                if(i==j)                    lMatrix.mat[i][j]=1;                if(j>i)                    lMatrix.mat[i][j]=0;            }        }        cout<<"lMatrix"<<endl;        print(lMatrix);        Matrix multsMatrix;        multsMatrix.initMatrix(&multsMatrix,lMatrix.col,uMatrix.row);        matrix.multsmatrix(&multsMatrix,lMatrix,uMatrix);        cout<<"lu"<<endl;        print(multsMatrix);        //计算u逆        for(j=0; j<uMatrix.row; j++)        {            for(i=j; i>=0; i--)            {                if(i==j)                    uniMatrix.mat[i][j]=1/uMatrix.mat[i][j];                else                {                    temp=0;                    for(k=j; k>i; k--)                    {                        temp+=uMatrix.mat[i][k]*uniMatrix.mat[k][j];                    }                    uniMatrix.mat[i][j]=-1/uMatrix.mat[i][i]*temp;                }            }        }        cout<<"uniMatrix"<<endl;        print(uniMatrix);        //Matrix multsMatrix;        //matrix.multsmatrix(&multsMatrix,uMatrix,uniMatrix);        //cout<<"multsMatrix"<<endl;        //print(multsMatrix);        //计算l逆        for(j=0; j<lMatrix.row; j++)        {            for(i=0; i<lMatrix.col; i++)            {                if(j==i)                    lniMatrix.mat[i][j]=1;                else                {                    temp=0;                    for(k=j; k<i; k++)                    {                        temp+=(lMatrix.mat[i][k]*lniMatrix.mat[k][j]);                    }                    lniMatrix.mat[i][j]=-temp;                }            }        }        cout<<"lniMatrix"<<endl;        print(lniMatrix);        multsmatrix(&multsMatrix,uniMatrix,lniMatrix);        cout<<"luni"<<endl;        print(multsMatrix);        //initMatrix(niMatrix,multsMatrix.col,multsMatrix.row);        copy(multsMatrix,niMatrix);        multsmatrix(&multsMatrix,cpMatrix,*niMatrix);        cout<<"luluni"<<endl;        print(multsMatrix);        copy(cpMatrix,&matrix);    }};
0 0