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