ZuoZhihua
哈爾濱工程大學(xué) 船舶工程學(xué)院
2020/8/13 星期四 午 永新 ThinkPad E485 Elementary os 5.0
QQ：1325686572
1 源碼
module operator_i
interface operator(.i.) !!// 自定義重載矩陣求逆操作符
module procedure brinv !!// 實(shí)矩陣求逆
module procedure bcinv !!// 復(fù)矩陣求逆
end interface
contains
!!// 實(shí)矩陣求逆核心代碼 trans from 徐士良《Fortran常用算法集》
!!// 作者：zuozhihua 時(shí)間：2020/8/18 地點(diǎn)：江西
function brinv(re) result(r)
real*8,intent(in) :: re(:,:) !!// 原矩陣
real*8 :: r(size(re,1),size(re,1)) !!// 逆矩陣
integer*4 :: flag,n !!// flag判斷奇異性额划；n是矩陣維度
real*8 :: t,d !!// 中間變量
integer*4 :: is(size(re,1)),js(size(re,1)) !!// 中間變量
n=size(re,1) !!// 獲取矩陣維度
r=re !!// 復(fù)制值
flag=1
do k=1,n
d=0.0
do i=k,n
do j=k,n
if (abs(r(i,j)).gt.d) then
d=abs(r(i,j))
is(k)=i
js(k)=j
end if
end do
end do
if (d+1.0.eq.1.0) then
flag=0
write(*,*) "flag=0,實(shí)矩陣奇異踊兜！"
return !!// 矩陣奇異，退出逆矩陣求解
end if
do j=1,n
t=r(k,j)
r(k,j)=r(is(k),j)
r(is(k),j)=t
end do
do i=1,n
t=r(i,k)
r(i,k)=r(i,js(k))
r(i,js(k))=t
end do
r(k,k)=1/r(k,k)
do j=1,n
if (j.ne.k) r(k,j)=r(k,j)*r(k,k)
end do
do i=1,n
if (i.ne.k) then
do j=1,n
if (j.ne.k) r(i,j)=r(i,j)-r(i,k)*r(k,j)
end do
end if
end do
do i=1,n
if (i.ne.k) r(i,k)=-r(i,k)*r(k,k)
end do
end do
do k=n,1,-1
do j=1,n
t=r(k,j)
r(k,j)=r(js(k),j)
r(js(k),j)=t
end do
do i=1,n
t=r(i,k)
r(i,k)=r(i,is(k))
r(i,is(k))=t
end do
end do
end function
!!// 復(fù)矩陣求逆核心代碼 trans from 徐士良《Fortran常用算法集》
!!// 作者：zuozhihua 時(shí)間：2020/8/10 地點(diǎn)：江西
function bcinv(cpx)
complex*16,intent(in) :: cpx(:,:) !!// 原矩陣
complex*16 :: bcinv(size(cpx,1),size(cpx,2)) !!// 逆矩陣
integer*4 :: flag,n !!// flag判斷奇異性；n是矩陣維度
real*8 :: ar(size(cpx,1),size(cpx,1)),ai(size(cpx,1),size(cpx,1)) !!// 實(shí)部矩陣ar族吻；虛部矩陣ai
real*8 :: d,p,t,q,s,b !!// 中間變量
integer*4 :: is(size(cpx,1)),js(size(cpx,1)) !!// 中間變量
n=size(cpx,1)
forall(i=1:n,j=1:n)
ar(i,j) = real(cpx(i,j));ai(i,j) = imag(cpx(i,j))
end forall
flag=1
do k=1,n
d=0.0
do i=k,n
do j=k,n
p=ar(i,j)*ar(i,j)+ai(i,j)*ai(i,j)
if(p.gt.d) then
d=p
is(k)=i
js(k)=j
end if
end do
end do
if(d+1.0.eq.1.0) then
flag=0
write(*,*) "flag=0,復(fù)矩陣奇異！"
return !!// 矩陣奇異珠增，退出逆矩陣求解
end if
do j=1,n
t=ar(k,j)
ar(k,j)=ar(is(k),j)
ar(is(k),j)=t
t=ai(k,j)
ai(k,j)=ai(is(k),j)
ai(is(k),j)=t
end do
do i=1,n
t=ar(i,k)
ar(i,k)=ar(i,js(k))
ar(i,js(k))=t
t=ai(i,k)
ai(i,k)=ai(i,js(k))
ai(i,js(k))=t
end do
ar(k,k)=ar(k,k)/d
ai(k,k)=-ai(k,k)/d
do j=1,n
if(j.ne.k) then
p=ar(k,j)*ar(k,k)
q=ai(k,j)*ai(k,k)
s=(ar(k,j)+ai(k,j))*(ar(k,k)+ai(k,k))
ar(k,j)=p-q
ai(k,j)=s-p-q
end if
end do
do i=1,n
if(i.ne.k) then
do j=1,n
if (j.ne.k) then
p=ar(k,j)*ar(i,k)
q=ai(k,j)*ai(i,k)
s=(ar(k,j)+ai(k,j))*(ar(i,k)+ai(i,k))
t=p-q
b=s-p-q
ar(i,j)=ar(i,j)-t
ai(i,j)=ai(i,j)-b
end if
end do
end if
end do
do i=1,n
if (i.ne.k) then
p=ar(i,k)*ar(k,k)
q=ai(i,k)*ai(k,k)
s=(ar(i,k)+ai(i,k))*(ar(k,k)+ai(k,k))
ar(i,k)=q-p
ai(i,k)=p+q-s
end if
end do
end do
do k=n,1,-1
do j=1,n
t=ar(k,j)
ar(k,j)=ar(js(k),j)
ar(js(k),j)=t
t=ai(k,j)
ai(k,j)=ai(js(k),j)
ai(js(k),j)=t
end do
do i=1,n
t=ar(i,k)
ar(i,k)=ar(i,is(k))
ar(i,is(k))=t
t=ai(i,k)
ai(i,k)=ai(i,is(k))
ai(i,is(k))=t
end do
end do
forall(i=1:n,j=1:n)
bcinv(i,j) = cmplx(ar(i,j),ai(i,j),8)
end forall
end function
end module
program main
use operator_i
implicit none
real*8 :: r(2,2),ri(2,2)
integer*4 :: i,j
complex*16 :: c(4,4),ci(4,4)
c=reshape((/(0.2368,0.1345),(1.1161,1.2671),(0.1582,-0.2836),(0.1968,0.3576),&
(0.2471,0.1678),(0.1254,0.2017),(1.1675,-1.1967),(0.2071,-1.2345),&
(0.2568,0.1825),(0.1397,0.7024),(0.1768,0.3558),(1.2168,2.1185),&
(1.2671,1.1161),(0.1490,0.2721),(0.1871,-0.2078),(0.2271,0.4773)/),(/4,4/))
ci=.i.c
write(*,*) "復(fù)矩陣："
write(*,"(8(4x,es10.3))") ((c(i,j),j=1,4),i=1,4)
write(*,*) "逆矩陣："
write(*,"(8(4x,es10.3))") ((ci(i,j),j=1,4),i=1,4)
write(*,*) "校核："
write(*,"(8(4x,es10.3))") matmul(c,ci)
r=reshape((/1,3,2,4/),(/2,2/))
ri=.i.r
write(*,*) "實(shí)矩陣："
write(*,"(2(4x,es10.3))") ((r(i,j),j=1,2),i=1,2)
write(*,*) "逆矩陣："
write(*,"(2(4x,es10.3))") ((ri(i,j),j=1,2),i=1,2)
write(*,*) "校核："
write(*,"(2(4x,es10.3))") matmul(r,ri)
end program
2 運(yùn)行
程序運(yùn)行圖.png
參考文獻(xiàn)
[1] 徐士良超歌，F(xiàn)ortran常用算法集。
[2] 白海波蒂教，F(xiàn)ortran程序設(shè)計(jì)權(quán)威指南巍举。
[3] zuozhihua，Fortran復(fù)數(shù)矩陣相乘凝垛。
