gFortran:LU分解求矩阵逆&编译工程

2023-05-15#电脑 #Fortran #学习

计算流体力学课的编程作业。做起来发现：真的难啊！自己的基础是真的差啊！

gFortran

fortran算是上古语言了，网上的资料和代码不多，但许多科学计算软件都用这个，包括前段时间报道的将要失去维护的粒子物理计算软件FORM。

这玩意的编译器有很多，我用的是gFortran，比起Intel Fortran相对轻量。（参考：选择编译器）用这个编译工程时，发现学习成本较高，所幸找到篇类似的博客，解决了。（不过要更方便的话，还是建议写Makefile，可惜我还没学会）方法如下：

建议用英文路径，别放在windows下的C盘用户文件夹下，避免奇怪的错误。将所有代码内的反斜杠(\)批量替换成斜杠(/)。如果用正则表达式，记得转义(“\“–>”/“)。
gfortran -c main.f -o main.o
默认生成main.o。-c命令表示编译但不链接文件。可以加上-o选项，指定输出文件。eg. gfortran -c main.f -o main2.o
gfortran -o ns2d_main main2.o flux.f force.f funtion.f grid.f high_order.f init.f input.f ismoo.f multigrid.f output.f solve_time.f step.f turblence.f viscous.f boundary.f
将main.o链接到ns2d_main。这个名字可以自行指定。若编译报错，看看缺的是什么函数，把对应的源码文件写上就行。参考
./ns2d_main
运行当前路径下的可执行文件。linux下任何文件都可运行。不写./的话，bash会找对应的命令，只有写在alias文件里才能找到。
这样可以给不同版本指定不同文件名，能编译好几版共存，比较方便。

LU分解求矩阵逆

在写这个的时候遇到不少问题。虽然纸上写起来不难，但写成程序，逻辑问题比较大。Fortran的列优先操作让人很不习惯，默认下标从1开始也使得移植C/C++程序比较麻烦。所幸有算法网的教程和C语言的程序可以参考。这里许多教程都是有问题的，包括网上的计算器。我调了很久才发现，当时调不出还挺郁闷的。

关于数学操作的过程，建议用英文搜，找国外大学的课件，写得清晰、大方。

分配

从文件读入数据，然后分配可变大小的数组。注意，这只是程序片段。

integer:: m,n,l,i,j,k
real*8::s
real*8,allocatable ::LM(:,:),LL(:,:),C(:,:),U(:,:),UM(:,:),UI(:,:),A(:,:)
INTEGER::ERR_MESSAGE
open(1,file="matri2.txt")
read(1,*)!第1行标示参数名，读入后丢弃
read(1,*)m,n
ALLOCATE(U(n,m),LL(n,m),C(n,m),Um(n,m),LM(n,m),UI(n,m),A(n,m),STAT=ERR_MESSAGE)!row first
IF(ERR_MESSAGE.NE.0) then 
    print *,"ALLOCATION ERROR"
end if

要先声明再分配。然后读入数组即可。

LU分解

在原数组上做LU分解。让$A=LU$。L数组的名字是LL，U数组的是U。fortran里没有大小写之分。

do i=1,m
    LL(i,i)=1
end do
do i=1,m-1
    do j=i+1,m
        LL(i,j)=U(i,j)/U(i,i)
        do k=1,m
            U(k,j)=U(k,j)-LL(i,j)*U(k,i)
        end do
    end do
end do

把分解完的矩阵再乘起来，检查是否能得原矩阵。然后分别对这俩矩阵求逆，得LM、UM。L是下三角，把对上三角矩阵求逆的代码里的行列调换即可，也就是求$(L^T)^{-1}$。这里调了很久，下标差1，移植的难度就加大了。一开始不得法，搞了通宵也不成，后来每次输出当前的下标，和C++的程序对照着调，直到每次访问的位置一致，下标差1。

需要注意的是，即使L的对角线上元素是1，直接对非对角元素取负也不行。4*4的矩阵就会有问题。

! L-1 for diagonal=1
! Note: simply change items negative is wrong.
! do i=1,m
!     lm(i,i)=1.0
!     do j=1,i-1
!         lm(j,i)=-LL(j,i)
!     end do
! end do
do i=1,m
    lm(i,i)=1.0/ll(i,i)
    do j=i-1,1,-1
        s=0.0!注意看警告，输出全是0时小心数据类型问题
        do k=j+1,i!调死我了! 和C++程序对照ijk调
            s=s+ll(j,k)*lm(k,i)!upper triangle
            ! write(*,*)i,j,k
        end do
        lm(j,i)=-s/ll(j,j)
    end do
end do
print *,"L^-1 Matrix"
do i=1,m
    write(*,*) Lm(:,i)!第一个指标优先输出，或者a(i,:)，第二个指标优先输出
end do
! U-1 first Transverse then treat as low triangle
! do i=m,1,-1
!     um(i,i)=1.0/u(i,i)
!     do j=i-1,1,-1
!         s=0
!         do k=m,i,-1
!             um(k,j)=um(k,j)-um(k,i)*u(k,j)
!         end do
!     end do

!     do j=i+1,m
!         um(j,i)=um(j,i)/u(i,i)
!     end do
! end do
do i=1,m
    um(i,i)=1.0/u(i,i)
    do j=i-1,1,-1
        s=0.0!注意看警告，输出全是0时小心数据类型问题
        do k=j+1,i!调死我了! 和C++程序对照ijk调
            s=s+u(k,j)*um(i,k)
            ! write(*,*)i,j,k
        end do
        um(i,j)=-s/u(j,j)
    end do
end do
print *,"U^-1 Matrix"
do l=1,m
    write(*,*) um(:,l)
end do

注释掉的是我自己写的，但跑不通。

调试的过程中，难免没有头绪，此时要多输出调试信息，比如数组下标和当前状态的矩阵。不然就会浪费大量时间在推测变量取值上。也要注意分解问题，比如最后求逆矩阵时求不对，就先检查L、U的逆是否有问题，算一下$LL^{-1}$，发现是L的逆的求法问题。

求逆矩阵：$A^{-1}=U^{-1}L^{-1}$。再乘回来看是否是单位阵即可。可能有浮点误差，不过无伤大雅。这里$A^{-1}$的变量名是UI。

print *,"U^-1*L^-1=A^-1"
do i=1,m
    do j=1,m
        do k=1,m
            UI(j,i)=ui(j,i)+um(k,i)*lm(j,k)
        end do
    end do
end do
do i=1,m
    write(*,*) ui(:,i)!第一个指标优先输出，或者a(i,:)，第二个指标优先输出
end do

在调试时，一直盯着屏幕发呆，甚至通宵并没有用。警惕这种虚假自律，多走动，多想想，再回来写。千万别生气，也别摧残自己的身体。

弦月之舞