Fortran数值计算——矩阵操作函数（可自用）
Fortran数值计算——矩阵操作函数（可自用）
引言
Fortran语言与现代计算机语言显得多少有点笨重，比如：用到的变量和数组需要声明大小、矩阵操作需要编制相应的程序才能实现相乘、求逆、等一系列不怎么人性化的操作，那么问题来了：现在还有必要学习Fortran吗？知乎上关于此类问题的回答不下几万个，有的人说：没必要，可以直接放弃，转向Python、Matlab；有的人说：有必要，对于工科生而言，有限元思想往往会伴随着自己的学术生涯，有限元的二次开发语言仍为Fortran，甚至当你沿用前人遗留下来的代码时，会不得不了解Fortran语言本身。我赞成后者的观点，因为木木在进行有限元架构的分析时，发现软件的操作很容易上手，但想要知道其背后的原理却需要更多的时间去学习理论及源程序的实现。为此，重新捧起了彭老的《Fortran95 程序设计》，再一次的学习与巩固。本篇推文旨在讲解Fortran中的矩阵操作（加减、乘、求逆、行列式），可直接移植于自己的相关程序中，用的时候call一下。
矩阵加减
代码区
subroutine matrix_Plus77(a,b,c,m,n) ! Fortran 77 矩阵相加
implicit none
dimension a(m,n),b(m,n),c(m,n)
double PRECISION a,b,c
integer m,n,row,col
do row=1,m
do col=1,n
c(row,col)=a(row,col)+b(row,col)
end do
end do
return
end
subroutine matrix_Plus90(a,b,c,m,n) ! Fortran 90 矩阵相加
implicit none
dimension a(m,n),b(m,n),c(m,n)
double PRECISION a,b,c
integer m,n
c=a+b
return
end总结
以上两种风格的代码均可在Simply Fortran中以.f 的文件进行编译，从以上代码风格中明显感觉到Fortran语言的“笨重性”，自由格式明显更加简洁，在进行元素加减时，直接c=a+b即可，在下面的代码中将更能体现自由格式的简洁性，但为了保留固定格式的风格（UMAT、UEL等相关子程序），还是会给出相应的源代码。
矩阵乘法
代码区
subroutine matrix_mul77(a,ar,ac,b,br,bc,c) ! Fortran77矩阵相乘
implicit none
dimension a(ar,ac),b(br,bc),c(ar,bc)
double PRECISION a,b,c
integer i,j,k,ar,ac,br,bc
do i=1,ar
do j=1,bc
c(i,j)=0
do k=1,ac
c(i,j)=c(i,j)+a(i,k)*b(k,j)
end do
end do
end do
return
end
subroutine matrix_mul90(a,ar,ac,b,br,bc,c) ! Fortran90矩阵相乘
implicit none
integer ar,ac,br,bc
dimension a(ar,ac),b(br,bc),c(ar,bc)
double PRECISION a,b,c
c=MATMUL(a,b)
return
end小结
• 在Fortran90（自由格式）中有专门的库函数：
MATMUL可以实现矩阵的乘法，固定格式需要编制相应的程序实现矩阵相乘，切记，不是对应元素相乘哦~•
c(i,j)=c(i,j)+a(i,k)*b(k,j)VSc(i,j)=c(i,j)+a(i,j)*b(i,j)• 其实Fortran有专门的链接库：IMSL，Simply Fortran不支持，想体验的话可以采用Visual Fortran的专业版哦
矩阵求逆
代码区
SUBROUTINE Matrix_Inv(a,c,n)
C get inverse matrix
! a(n,n) - array of coefficients for matrix A
! n - dimension
! output ...
! c(n,n) - inverse matrix of A
implicit none
integer,intent(in)::n
double precision,intent(in):: a(n,n)
double precision,intent(out):: c(n,n)
double precision L(n,n), U(n,n), b(n), d(n), x(n)
double precision coeff,tem_a(n,n)
integer i, j, k
tem_a = a
L=0.0D0
U=0.0D0
b=0.0D0
do k=1, n-1
do i=k+1,n
coeff=tem_a(i,k)/tem_a(k,k)
L(i,k) = coeff
do j=k+1,n
tem_a(i,j) = tem_a(i,j)-coeff*tem_a(k,j)
end do
end do
end do
do i=1,n
L(i,i) = 1.0D0
end do
do j=1,n
do i=1,j
U(i,j) = tem_a(i,j)
end do
end do
do k=1,n
b(k)=1.0D0
d(1) = b(1)
do i=2,n
d(i)=b(i)
do j=1,i-1
d(i) = d(i) - L(i,j)*d(j)
end do
end do
x(n)=d(n)/U(n,n)
do i = n-1,1,-1
x(i) = d(i)
do j=n,i+1,-1
x(i)=x(i)-U(i,j)*x(j)
end do
x(i) = x(i)/u(i,i)
end do
do i=1,n
c(i,k) = x(i)
end do
b(k)=0.0D0
end do
return
END 小结
• 从以上代码可知，在Simply Fortran中注释符号既支持
C也支持!• 求逆运算时采用的是数值分析中的LU分解法
矩阵行列式
代码区
! 全选主元高斯消去法(Gauss-Jordan)求矩阵的行列式
SUBROUTINE Matrix_Det(A,N,DET)
DIMENSION A(N,N)
DOUBLE PRECISION A,DET,F,D,Q
F=1.0
DET=1.0
DO 100 K=1,N-1
Q=0.0
DO 10 I=K,N
DO 10 J=K,N
IF (ABS(A(I,J)).GT.Q) THEN
Q=ABS(A(I,J))
IS=I
JS=J
END IF
10 CONTINUE
IF (Q+1.0.EQ.1.0) THEN
DET=0.0
RETURN
END IF
IF (IS.NE.K) THEN
F=-F
DO 20 J=K,N
D=A(K,J)
A(K,J)=A(IS,J)
A(IS,J)=D
20 CONTINUE
END IF
IF (JS.NE.K) THEN
F=-F
DO 30 I=K,N
D=A(I,JS)
A(I,JS)=A(I,K)
A(I,K)=D
30 CONTINUE
END IF
DET=DET*A(K,K)
DO 50 I=K+1,N
D=A(I,K)/A(K,K)
DO 40 J=K+1,N
40 A(I,J)=A(I,J)-D*A(K,J)
50 CONTINUE
100 CONTINUE
DET=F*DET*A(N,N)
RETURN
END小结
• 该程序复刻的徐士良的《Fortran常用算法集 第二版》
主程序Test
代码区
C Matrix Operations
C
C 微 信 公 众 号：易木木响叮当
C 知乎：易木木响叮当
C B站：易木木响叮当
C
C By Yimumu 2022/5/26
program main
implicit none
dimension
a(4,4),b(4,4),c_mul77(4,4),c_mul90(4,4),a_inv(4,4),
&c_Plus77(4,4),c_Plus90(4,4)
double PRECISION
a,b,c_mul77,c_mul90,c_Plus77,c_Plus90,det,a_inv
integer i
data a/3.0,5.0,11.0,5.0,
& -3.0,-5.0,8.0,-1.0,
& -2.0,1.0,5.0,-3.0,
& 4.0,8.0,-7.0,-1.0/
data b/1.0,-2.0,0.0,3.0,
& 3.0,-1.0,8.0,-3.0,
& -2.0,5.0,4.0,2.0,
& 0.0,-7.0,1.0,-4.0/
call matrix_Plus77(a,b,c_Plus77,4,4)
call matrix_Plus90(a,b,c_Plus90,4,4)
call matrix_mul77(a,4,4,b,4,4,c_mul77)
call matrix_mul90(a,4,4,b,4,4,c_mul90)
call Matrix_Det(a,4,det)
call Matrix_Inv(a,a_inv,4)
write(*,*)"Matrix A_Det:"
write(*,10)det
10 format(1x,'det=',D15.6)
c 矩阵形式输出
write(*,*)"Matrix A:"
do i=1,4
write(*,*)a(:,i)
end do
write(*,*)"Matrix B:"
do i=1,4
write(*,*)b(:,i)
end do
write(*,*)"Matrix C_plus77:"
do i=1,4
write(*,*)c_Plus77(:,i)
end do
write(*,*)"Matrix C_plus90:"
do i=1,4
write(*,*)c_Plus90(:,i)
end do
write(*,*)"Matrix C_mul77:"
do i=1,4
write(*,*)c_mul77(:,i)
end do
write(*,*)"Matrix C_mul90:"
do i=1,4
write(*,*)c_mul90(:,i)
end do
write(*,*)"Matrix A_Inv:"
do i=1,4
write(*,*)a_inv(:,i)
end do
end计算结果
由以上结果可知，使用两种风格的代码，计算结果是一致的，且都能以.f文件编译成功，其结果与Matlab结果一致，可自行验证。
最后想给大家说的是，在面对二次开发时不要怕，要对代码有感情，不会的就去学，今天学一点，明天学一点，成长在点滴。以上源代码皆可使用在自己的程序中当作一个SUBROUTINE,随时调用。
