C言語求行列の行列式、伴随行列、逆行列

3466 ワード

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次になると計算できないような気がします.1次を増やすたびにアルゴリズムの複雑さが倍増します.
残念ながら計算速度が遅すぎて、また大量のデータを迅速に処理する必要があって、アルゴリズムを最適化するしかなくて、本文は絶対に怠け者で、既成の絶対に自分で書かないことを見つけることができます!ネット上で検索して、LU分解はもっと効率的なようで、結果は1篇の頼りになるドキュメントを見つけました:クリックしてリンクを開けます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は行列式の1つの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

上記のコードはマトリックスのLU分解を実現して、この1段を書くのはとても苦労して、幸い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
",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
",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
",i,j,dTemp3[i][j]); } }
はい、ここに記録して、残念ながらCSDNは直接添付ファイルをアップロードすることができなくて、さもなくば直接書いた分類関数のCとHファイルを載せます.