1. 程式人生 > >C語言求矩陣的行列式、伴隨矩陣、逆矩陣

C語言求矩陣的行列式、伴隨矩陣、逆矩陣

CSDN大神編寫的求矩陣的行列式,int getA(int arcs[N][N],int n),通過呼叫遞迴函式,按矩陣的第一行進行分解,雖然行列式的計算都學過,但是自己寫起來還是得費一番功夫的,好在有MATLAB可以驗證結果,結果對拿過來就可以直接用。

voidgetAStart(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檔案放上來。