plu-分解是對lu分解的一種改進，其增加了選主元的操作增加了計算的穩定性，及在第i次迴圈中將
j =w
here
(max(
∣a[i
:n,i
]∣))
j=where(\max(|a[i:n,i]|))
j=wher
e(max(∣a
[i:n
,i]∣
))行和第i行進行交換來比避免對角元素出現0的情況，計算結果
p a=
lu
pa=lu
pa=l
up為置換矩陣，l為下三角矩陣，u為上三角矩陣。選主元操作在計算過程中以乙個一維陣列儲存代替n×n
n\times n
n×n的矩陣，以下為此演算法的fortran**，a為輸入矩陣，l,u,p分別為下三角，上三角，置換矩陣，n為矩陣大小，當erro為0時計算失敗，成功則為1
subroutine plu_decompose(a,l,u,p,n,erro)
real,intent(inout)::a(n,n)
real,intent(out)::l(n,n),u(n,n),p(n,n)
integer,intent(in)::n
integer::pi(n)
integer,intent(out)::erro
integer::i,j,k,maxp,tempp
real::max
do i=1,n
pi(i)=i
end do
do i=1,n-1
max=a(i,i)
maxp=i
do j=i+1,n
if (abs(a(j,i))>max) then
max=abs(a(j,i))
maxp=j
end if
end do
if (max+1.0==1.0) then
erro=0
stop
end if
if (maxp/=i) then
tempp=pi(i)
pi(i)=pi(maxp);
pi(maxp)=tempp
a(i:maxp:maxp-i,:)=a(maxp:i:i-maxp,:)
end if
do j=i+1,n
a(j,i)=a(j,i)/a(i,i)
end do
do j=i+1,n
do k=i+1,n
a(j,k)=a(j,k)-a(j,i)*a(i,k)
end do
end do
end do
do i=1,n
p(i,i)=0
l(i,i)=1
u(i,i)=a(i,i)
do j=1,i-1
p(i,j)=0
p(j,i)=0
l(i,j)=a(i,j)
l(j,i)=0
u(j,i)=a(j,i)
u(i,j)=0
end do
end do
do i=1,n
p(i,pi(i))=1
end do
end subroutine
aa 的逆，**如下，其中第一次得到的矩陣為乙個下三角陣，這樣可以節省一半的計算量
subroutine reverse(a,b,n,erro)
implicit none
real,intent(inout)::a(n,n)
real,intent(out)::b(n,n)
real::l(n,n),u(n,n),p(n,n),c(n,n),temp
integer::n,erro
integer::i,j,k
call plu_decompose(a,l,u,p,n,erro)
if (erro==0) then
return
end if
do i=1,n
do j=1,i-1
c(j,i)=0.0
end do
do j=i,n
temp=0
do k=1,j-1
temp=temp+l(j,k)*c(k,i)
end do
if (i==j)then
c(j,i)=(1.0-temp)/l(j,j)
else
c(j,i)=(0.0-temp)/l(j,j)
end if
end do
end do
do i=1,n
do j=n,1,-1
temp=0
do k=n,j+1,-1
temp=temp+b(k,i)*u(j,k)
end do
b(j,i)=(c(j,i)-temp)/u(j,j)
end do
end do
b=matmul(b,p)
end subroutine
program main
implicit none
real::a(4,4),a_(4,4)
real::l(4,4)
real::w(4,4)
integer :: i,j
integer:: u
data((a(i,j),i=1,4),j=1,4)/ 2,0,2,0.6,3,3,4,-2,5,5,4,2,-1,-2,3.4,-1 /
data((a_(i,j),i=1,4),j=1,4)/ 2,0,2,0.6,3,3,4,-2,5,5,4,2,-1,-2,3.4,-1 /
a=transpose(a)
a_=transpose(a_)
do i=1,4
print *,a(i,:)
end do
call reverse(a,l,4,u)
w=matmul(a_,l)
print *,''
do i=1,4
print *,w(i,:)
end do
pause
end program main
c 矩陣求逆的lu分解實現
c 矩陣求逆的lu分解實現 初學c 嘗試利用基礎知識編寫矩陣求逆。但發現求解伴隨陣的過程非常複雜且難以實現，而我正好看到一篇求三角陣伴隨矩陣的文章，故嘗試程式設計實現。在這種方法下，計算量明顯減小，實現方法，思路適合初學者。參考文獻 三角形矩陣求伴隨矩陣的一種方法 曾月新 求逆矩陣思路 1.求矩陣的...
伴隨矩陣求逆矩陣
在之前的文章 線性代數之矩陣 中已經介紹了一些關於矩陣的基本概念，本篇文章主要就求解逆矩陣進行進一步總結。我們先看例子來直觀的理解什麼是余子式 minor，後邊將都用英文minor，中文的翻譯較亂 這個例子 我們假設矩陣為a 中我們看到a 1,1 的minor就是將a 1,1 所在的行和列刪除後剩下...
矩陣的求逆
最近做乙個加密演算法遇到需要計算矩陣的逆，閒著無聊，記錄一下，以後免得再麻煩。include include include define max 20 define e 0.000000001 計算矩陣src的模 double calculate a double src max int n fo...