矩阵乘法及两种求逆方法

在做TPS时用到了两种矩阵求逆的方法，都是适合大型矩阵求逆的方法，一个是高斯约旦法，还有是依照大一线性代数中的增广矩阵法求解，复杂度都是n的三次方。

#include “matrix_diy.h“

using namespace std;

double ** muti（double **Mint M_rowint M_coldouble **Nint N_rowint N_col）//矩阵乘法
{
	double **ans;
	ans=new double *[M_col];  //分配空间
	for（int i=0;i		ans[i]=new double[N_row];
	if（M_row!=N_col）    //出错
		return NULL;
	for（int i=0;i	{  
        for（int j=0;j		{  
            ans[i][j]=0;//初始化 
            for（int k=0;k                ans[i][j]+=M[i][k]*N[k][j];  
            }  
        }  
    } 
	return（ans）;
}

double **T（double **Mint M_rowint M_col）//矩阵转置
{
	double **ans;
	ans=new double *[M_row];  //分配空间
	for（int i=0;i		ans[i]=new double[M_col];
	for（int i=0;i		for（int j=0;j			ans[i][j]=M[j][i];
	return （ans）;
}

double **inverse1（double **M1int M1_col）//矩阵的逆用增广矩阵求解
{
	int n=M1_col;
	double**ans;     //为结果分配空间，作为返回值
	ans=new double *[n];
	for（int i=0;i		ans[i]=new double[n];
	double**M;     //为结果分配空间，作为中间值
	M=new double *[n];
	for（int i=0;i		M[i]=new double[2*n];   //增广矩阵
	for（int i=0;i	{
		for（int j=0;j		{
			M[i][j]=M1[i][j]; 
		}
		for（int j=n;j<2*n;j++）
		{
			if（i==j-n）M[i][j]=1;
			else M[i][j]=0;
		}
	}
	int flag1=0;
	int flag2=0;
	for（int k=0;k	{
		if（M[k][k]==0.0）  //零
		{
			for（int s=k+1;s			{
			   if（M[s][k]==0.0）continue;
			   else  
				{// 找到不为零的数
					flag2=s;
					double tmp=0;
					//和该行做行交换
		            for（int j=k;j<2*n;j++）
	                 {
		                 tmp=M[k][j];
		                 M[k][j]=M[flag2][j];
		                 M[flag2][j]=tmp;
		              }
					break;
			    }
			}
			if（M[k][k]==0.0）continue;
		}
		if（fabs（M[k][k]）==0.0）  //近似于零
			cout<<“ERROR“<<“奇异“<		for（int i=k+1;i		{
			double q=M[i][k]/M[k][k];
			for（int j=k;j<2*n;j++）
				M[i][j]=M[i][j]-M[k][j]*q;    // 先变成上三角矩阵
			for（int j=0;j			{
			    if（fabs（M[i][j]）<1e-10）      //过小则近似于0
					M[i][j]=0;
			}
		}
	}

	for（int k=n-1;k>0;k--）
	{
		if（M[k][k]==0.0）  //主元为零
		{
			for（int s=k-1;s>0;s--）
			{
			   if（M[s][k]!=0.0）
				{// 找到不为零的数
					flag1=s;
					double tmp=0;
					//和该行做行交换
		            for（int j=0;j<2*n;j++）
	                 {
		                 tmp=M[k][j];
		                 M[k][j]=M[flag1][j];
		                 M[flag1][j]=tmp;
		              }
					break;
			    }
			}
			if（fabs（M[k][k]）==0.0）continue;
		}
		for（int i=k-1;i>=0;i--）
		{
			double q=M[i][k]/M[k][k];
			for（int j=k;j<2*n;j++）
				M[i][j]=M[i][j]-M[k][j]*q;    // 再变成全是主元的矩阵
			for（int j=0;j			{
			    if（fabs（M[i][j]）<1e-10）      //过小则近似于0
					M[i][j]=0;
			}
		}
	}
	for（int i=0;i	{
		double ik=M[i][i];
		for（int j=0;j<2*n;j++）
			M[i][j]=M[i][j]/ik;                       //左边的矩阵变成I矩阵了，则右边的矩阵即为逆
	}
	for（int i=0;i		for（int j=0;j			ans[i][j]=M[i][j+n];

  for （int i=0;i

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件         377  2018-11-12 20:25  matrix_diy.h
     文件        5822  2018-11-20 21:06  matrix_diy.cpp