C语言求矩阵的行列式、伴随矩阵、逆矩阵

CSDN大神编写的求矩阵的行列式，int getA(int arcs[N][N],int n)，通过调用递归函数，按矩阵的第一行进行分解，虽然行列式的计算都学过，但是自己写起来还是得费一番功夫的，好在有MATLAB可以验证结果，结果对拿过来就可以直接用。

void getAStart(int arcs[N][N],int n,int ans[N][N])，求矩阵的伴随矩阵，也需要调用求矩阵的行列式。原文地址：

点击打开链接http://blog.csdn.net/abcjennifer/article/details/6693612

验证算法计算速度，当我用C调用去求一个10阶矩阵的行列式和伴随矩阵时，跑10条数据大概用7秒的时间（自己其他函数占用很小时间），但是当矩阵变成15阶时，感觉就计算不出来了，每增加一阶，算法的复杂度都成倍的增长。

无奈计算速度太慢，又需要快速处理大量数据，只好优化算法了，本文绝对算的上是懒人，能找到现成的绝对不自己写！网上搜了下，LU分解貌似效率更高一点，结果找到一篇靠谱的文档：点击打开链接http://www.docin.com/p-690103638.html

文献对矩阵分解的原理进行了详细介绍，当初矩阵分析学的不好，也没仔细看，只可惜找到它太晚，以至于找到他时，我已经照着下面的MATLAB代码把LU分解的C代码写完了，

MATLAB代码如下：

function [L,U]=myLU(A)  
 
%实现对矩阵A的LU分解，L为下三角矩阵  
 

 
[n,n]=size(A);  
 
L=zeros(n,n);  
 
U=zeros(n,n);  
 
for i=1:n  
 
L(i,i)=1;  
 
end  
 
    for k=1:n  

        for j=k:n  

        U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)');  


        %U(k,j)
        end  

        for i=k+1:n  

        L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k);  
        %L(i,k)
        end  



    end  

短短几行代码让人用C写的那叫一个痛苦。

参考那个文档把LU 分解后L矩阵（下三角）和U矩阵（上三角）求逆的直接展过来用，对着MATLAB结果，终于对了，不容易啊。用C写程序觉对是需要耐心，MATLAB求行列式一个det就可以了，求逆inv,求伴随是啥命令给忘了，不过直接掉det(A)*inv(A),  得出来的就是伴随矩阵。

验证结果很简单，A=L*U,那么inv(A)=inv(U)*inv(L);

搞定

void myLU(double A[COLUMN][COLUMN],double Low[COLUMN][COLUMN],double Up[COLUMN][COLUMN],int n)
{
	int i,j,k,t,q;
	//double Low[COLUMN][COLUMN]={0};
	double temp1[COLUMN]={0};
	double temp2[COLUMN]={0};
	double temp3=0;
	for (i=0;i<n;i++){
		for(j=0;j<n;j++){
			if(i==j)
		Low[i][i]=1;
		}
	}
	/*Low[0][0]=1;Low[1][1]=1;Low[2][2]=1;Low[3][3]=1;Low[4][4]=1;
	Low[5][5]=1;Low[6][6]=1;Low[7][7]=1;Low[8][8]=1;*/
	
	for(k=0;k<n;k++){
		for(j=k;j<n;j++){

		//U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)');  
		/*fetch_matrix_L(Low,temp1,k);
		fetch_matrix_U(Up,temp2,k,j)*/;
		if(k==0){
			temp3=0;
		}
		else{
				for(i=0;i<=k-1;i++){
				temp1[i]=Low[k][i];//L(k,1:k-1)
				}		
				for(i=0;i<=k-1;i++){
				temp2[i]=Up[i][j];//U
				}		
				temp3=0;
				for(i=0;i<=k-1;i++){
				temp3+=temp1[i]*temp2[i];//点乘，求和
				}
		}
		
			
		Up[k][j]=A[k][j]-temp3;
		//printf("%d----%d-----%f\n",k,j,Up[k][j]);
		}
		//下三角
		memset(temp1,0,sizeof(double));
		memset(temp2,0,sizeof(double));
		temp3=0;
		for(t=k+1;t<n;t++){
		/*fetch_matrix_L(Low,temp1,k);
		fetch_matrix_U(Up,temp2,k,k);*/
			if(k==0){
				temp3=0;
			}
			else
			{
				for(i=0;i<=k-1;i++){
				temp1[i]=Low[t][i];//L(k,1:k-1)
				}		
				for(i=0;i<=k-1;i++){
				temp2[i]=Up[i][k];//U
				}		
				temp3=0;
				for(i=0;i<=k-1;i++){
				temp3+=temp1[i]*temp2[i];//点乘，求和
				}

			}
			Low[t][k]=(A[t][k]-temp3)/Up[k][k];
			//printf("%d----%d-----%f\n",k,j,Low[t][k]);
		}
	
	}


}

上述代码实现矩阵的LU分解，写这一段可费了不少劲，好在有MATLAB，结果不对，就一步步的跟，总能找到错误在哪里？大不了MATLAB里 打个断点,C里面我设个断点。每走一步观察一下结果。

矩阵分解完后后，便是求各自的逆运算了。也懒得单独写一个函数了，copy代码。

//--------------test-------------------------
	//dTemp=dTemp_test;复制测试数组,zuduan
	/*for(i=0;i<9;i++){
	memcpy(dTemp[i],dTemp_test[i], sizeof(double)*9);
	}*/
//-------------------------------------------
	//det_A=getA(dTemp,piRowLenth[3]);//|A|,被det_AAA取代
	//printf("%1.5e\n",det_A);//real data
//------------test LU 分解-------------------
	
	myLU(dTemp,Temp11,Temp22,9);//注意temp1与temp2下面是否复用
	det_AA=Temp11[0][0];
	for(i=1;i<9;i++){
	det_AA*=Temp11[i][i];
	}
	det_AAA=Temp22[0][0];
	for(i=1;i<9;i++){
	det_AAA*=Temp22[i][i];
	}
	det_AAA=det_AAA*det_AA;
//	printf("%1.5e\n",det_AAA);
	det_A=det_AAA;


//-----求L和U(Temp22)的逆------------------------------------------
	for(i=0;i<9;i++){
	u[i][i]=1.0/Temp22[i][i];
		for(k=i-1;k>=0;k--){
		dDataSun=0;
		for(j=k+1;j<=i;j++)
		dDataSun=dDataSun+Temp22[k][j]*u[j][i];
		u[k][i]=-dDataSun/Temp22[k][k];
		
		}
	}
//--------求L（temp11）------------------------------------
	for (i=0;i<9;i++){
		l[i][i]=1;
		for(k=i+1;k<9;k++){
			 for(j=i;j<=k-1;j++)
				l[k][i]=l[k][i]-Temp11[k][j]*l[j][i];
				}
	
	}

//--------------求L与U乘积
	for(i=0;i<9;i++){
		for(j=0;j<9;j++){
			for(k=0;k<9;k++)
		dTemp3[i][j]+=u[i][k]*l[k][j];
		//	printf("%d---%d----%f\n",i,j,dTemp3[i][j]);
		}
	}

好了，记录到这，可惜CSDN不能直接上传附件，不然直接把写的分类函数的C和H文件放上来。