[Fortran] 纯文本查看 复制代码
!> @details 下三角矩阵求逆(稠密存储)
!> @author Jackdaw毛毛
!> @date 2019 年 8 月 26 日
subroutine LowerTriMatrixInv(mat,matInv)
real(fp),intent(in) :: mat(:,:) !< 输入矩阵,下三角,稠密存储
real(fp),intent(inout) :: matInv(:,:) !< 输出矩阵,下三角,稠密存储,与输入矩阵关系:互为逆矩阵
integer :: i,j,n
n = size(mat,dim=1)
matInv = 0.0_fp
do i = 1, n
matInv(i,i) = 1.0_fp / mat(i,i)
end do
do i = 2, n
do j = 1, i-1
!do k = j, i-1
! matInv(i,j) = matInv(i,j) + mat(i,k)*matInv(k,j)
!end do
matInv(i,j) = dot_product( mat(i,j:i-1),matInv(j:i-1,j) )
matInv(i,j) = matInv(i,j) * (-matInv(i,i))
end do
end do
end subroutine
!> @details 高斯赛德尔迭代法计算 \f$ Ax=b \f$,收敛速度比雅克比迭代更快。
!> @author Jackdaw毛毛
!> @date 2019 年 8 月 26 日
subroutine useGaussSeidelIterSolveAx_eq_b(matrixA,vectorB,vectorX)
real(fp),intent(in) :: matrixA(:,:) !< 系数矩阵(迭代收敛的充分条件:严格对角占优或对称正定)
real(fp),intent(in) :: vectorB(:) !< 右端项
real(fp),intent(inout) :: vectorX(:) !< 解向量
real(fp),allocatable :: G(:,:),tmpM(:,:),tmpB(:),tmpX(:)
integer ::i,n
n = size(matrixA,dim=1)
allocate(G(n,n),tmpM(n,n),tmpX(n),tmpB(n))
! 迭代格式
! x(k+1) = D^-1(b-Ux(k)-Lx(k+1))
! x(k+1) = (L+D)^-1(-Ux(k) + b )
! 构造迭代矩阵
tmpM = 0.0_fp
do i = 1, n
tmpM(i:n,i) = matrixA(i:n,i) ! L+D
end do
call LowerTriMatrixInv(tmpM,G) ! 求逆,存在G里面
tmpM = 0.0_fp
do i = 1, n-1
tmpM(i,i+1:n) = matrixA(i,i+1:n) ! U
end do
tmpB = matmul(G,vectorB) ! 常向量
G = matmul(G,tmpM) ! 迭代矩阵
deallocate(tmpM)
do i = 1, 1000
tmpX = tmpB - matmul(G,vectorX)
if( maxval(abs(vectorX-tmpX)) .lt. 1.e-10_fp ) then
vectorX = vectorX + tmpX
vectorX = vectorX / 2.0_fp
exit
else
vectorX = tmpX
endif
end do
if(i .ge. 1000 ) write(*,*) "GaussSeidel-iteration failed."
deallocate( G,tmpX,tmpB )
end subroutine
!> @details 高斯赛德尔迭代法计算 \f$ Ax=b \f$,测试算例。
!> @author Jackdaw毛毛
!> @date 2019 年 8 月 26 日
subroutine useGaussSeidelIterSolveAx_eq_b_TEST
real(fp) :: a(2,2),b(2),x(2)
a(1,:) = [ 3.0_fp, 1.0_fp ]
a(2,:) = [ 1.0_fp, 2.0_fp ]
b = [ 5.0_fp, 5.0_fp ]
x = 0.0_fp
call useGaussSeidelIterSolveAx_eq_b(a,b,x)
write(*,*) "解析解 x1 = 1.0, x2 = 2.0"
write(*,*) "数值解",x
end subroutine useGaussSeidelIterSolveAx_eq_b_TEST