ZuoZhihua
哈爾濱工程大學(xué) 船舶工程學(xué)院
2020/8/27 星期四 午 永新 ThinkPad E485 Windows Home
微信：zozha96
1 Fortran矩陣求逆槽袄、相乘模塊
module operation_i
!!// operator ".i." is used for solve the inverse MATRIX of a given MATRIX.
public
interface operator(.i.) !!// 自定義重載矩陣求逆操作符
module procedure brinv !!// 實(shí)矩陣求逆
module procedure bcinv !!// 復(fù)矩陣求逆
end interface
private :: brinv, bcinv
contains
!!// 實(shí)矩陣求逆核心代碼 trans from 徐士良《Fortran常用算法集》
!!// 作者：zuozhihua 時(shí)間：2020/8/18 地點(diǎn)：江西
function brinv(re) result(r)
real,intent(in) :: re(:,:) !!// 原矩陣
real :: r(size(re,1),size(re,1)) !!// 逆矩陣
integer :: flag,n !!// flag判斷奇異性；n是矩陣維度
real :: t,d !!// 中間變量
integer :: 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,intent(in) :: cpx(:,:) !!// 原矩陣
complex :: bcinv(size(cpx,1),size(cpx,2)) !!// 逆矩陣
integer :: flag,n !!// flag判斷奇異性喷舀；n是矩陣維度
real :: ar(size(cpx,1),size(cpx,1)),ai(size(cpx,1),size(cpx,1)) !!// 實(shí)部矩陣ar砍濒；虛部矩陣ai
real :: d,p,t,q,s,b !!// 中間變量
integer :: 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))
end forall
end function
end module
module operation_x
!!// operator ".x." is used for the multiplicative operation between MATRIXS.
public
interface operator(.x.) !!// 自定義重載矩陣相乘操作符
module procedure rmut !!// 實(shí)矩陣相乘
module procedure cmut !!// 復(fù)矩陣相乘
module procedure rcmut !!// 實(shí)復(fù)矩陣相乘
module procedure crmut !!// 實(shí)復(fù)矩陣相乘
end interface
private :: rmut, cmut, rcmut, crmut
contains
!!// real m1,m2
function rmut(m1,m2) result(ret)
real,intent(in) :: m1(:,:),m2(:,:)
real :: ret(size(m1,1),size(m2,2))
ret=matmul(m1,m2)
end function
!!// cmplx*16 m1,m2
function cmut(m1,m2) result(ret)
complex,intent(in) :: m1(:,:),m2(:,:)
complex :: ret(size(m1,1),size(m2,2))
ret=matmul(m1,m2)
end function
!!// cmplx*16 & real
function rcmut(m1,m2) result(ret)
real,intent(in) :: m1(:,:)
complex,intent(in) :: m2(:,:)
complex :: ret(size(m1,1),size(m2,2))
ret=matmul(m1,m2)
end function
!!// real & cmplx*16
function crmut(m1,m2) result(ret)
real,intent(in) :: m2(:,:)
complex,intent(in) :: m1(:,:)
complex :: ret(size(m1,1),size(m2,2))
ret=matmul(m1,m2)
end function
end module
2 矩形相乘測(cè)試
program main
use operation_x
real :: rmat1(2,2),rmat2(2,2),rmat3(2,2)
complex :: cmat1(2,2),cmat2(2,2),cmat3(2,2)
rmat1=1; rmat2=2
cmat1=(3,3); cmat2=(4,4); cmat2(2,2)=(3,7)
rmat3=rmat1.x.rmat2; cmat3=cmat1.x.cmat2
write(*,*) "原實(shí)矩陣："
write(*,"(2(4x,es10.3))") ((rmat1(i,j),j=1,2),i=1,2)
write(*,"(2(4x,es10.3))") ((rmat2(i,j),j=1,2),i=1,2)
write(*,*) "結(jié)果實(shí)矩陣："
write(*,"(2(4x,es10.3))") ((rmat3(i,j),j=1,2),i=1,2)
write(*,*) "原復(fù)矩陣："
write(*,"(4(4x,es10.3))") ((cmat1(i,j),j=1,2),i=1,2)
write(*,"(4(4x,es10.3))") ((cmat2(i,j),j=1,2),i=1,2)
write(*,*) "結(jié)果復(fù)矩陣："
write(*,"(4(4x,es10.3))") ((cmat3(i,j),j=1,2),i=1,2)
read(*,*)
end program
2.1 VsCode程序運(yùn)行圖
編譯環(huán)境：GCC Fortran 10.0; VsCode 2020; Win10 Home Edition; ThinkPadE485.
Done
2.2 問(wèn)題
本程序模塊沒(méi)有驗(yàn)證矩陣是否符合乘法計(jì)算條件爸邢，如果3×4的矩陣與5×6的矩陣使用本程序模塊相乘，也會(huì)出現(xiàn)計(jì)算結(jié)果拿愧，但顯然是不對(duì)的杠河。
如果你想使用這個(gè)模塊需要注意這一點(diǎn)〗焦迹或者你可以對(duì)模塊代碼進(jìn)行更改券敌。
今天是2020/11/2日，有空會(huì)把這個(gè)問(wèn)題解決掉柳洋。
本程序模塊使用單精度待诅，如需使用雙精度，請(qǐng)自行調(diào)整熊镣。
祝好運(yùn)卑雁！
References
[1] 徐士良，F(xiàn)ortran常用算法集绪囱。
[2] 白海波测蹲，F(xiàn)ortran程序設(shè)計(jì)權(quán)威指南。
[3] 左志華毕箍，Fortran復(fù)數(shù)矩陣求逆弛房。
