ZuoZhihua
哈尔滨工程大学船舶工程学院
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 !!// 实矩阵求逆
        module procedure bcinv !!// 复矩阵求逆
    end interface
    private :: brinv, bcinv
contains
    !!// 实矩阵求逆核心代码 trans from 徐士良《Fortran常用算法集》
    !!// 作者：zuozhihua 时间：2020/8/18 地点：江西
    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 !!// 复制值
        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,实矩阵奇异！" 
                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
    !!// 复矩阵求逆核心代码 trans from 徐士良《Fortran常用算法集》
    !!// 作者：zuozhihua 时间：2020/8/10 地点：江西
    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)) !!// 实部矩阵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,复矩阵奇异！" 
                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 !!// 实矩阵相乘
        module procedure cmut !!// 复矩阵相乘
        module procedure rcmut !!// 实复矩阵相乘
        module procedure crmut !!// 实复矩阵相乘 
    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 矩形相乘测试

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(*,*) "原实矩阵："
    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(*,*) "结果实矩阵："
    write(*,"(2(4x,es10.3))") ((rmat3(i,j),j=1,2),i=1,2)

    write(*,*) "原复矩阵："
    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(*,*) "结果复矩阵："
    write(*,"(4(4x,es10.3))") ((cmat3(i,j),j=1,2),i=1,2)

    read(*,*)

end program

2.1 VsCode程序运行图

编译环境：GCC Fortran 10.0; VsCode 2020; Win10 Home Edition; ThinkPadE485.

Done

2.2 问题

本程序模块没有验证矩阵是否符合乘法计算条件，如果3×4的矩阵与5×6的矩阵使用本程序模块相乘，也会出现计算结果，但显然是不对的。
如果你想使用这个模块需要注意这一点。或者你可以对模块代码进行更改。
今天是2020/11/2日，有空会把这个问题解决掉。

本程序模块使用单精度，如需使用双精度，请自行调整。

祝好运！

References

[1] 徐士良，Fortran常用算法集。
[2] 白海波，Fortran程序设计权威指南。
[3] 左志华，Fortran复数矩阵求逆。

Fortran矩阵相乘（复矩阵、实矩阵）