Doolittle分解的算法实现问题
各位老师好，我在写用“Doolittle分解法解线性方程组”这一数值算法的实现程序时遇到了问题，程序给出了错误的结果，想请教各位老师这个程序在算法上有什么错误，谢谢！（附件中是实现doolittle分解法的子程序，这是我写的交互程序的一部分，A,X,B分别代表系数矩阵，解，常数矩阵；ndim代表方阵A的维度）
你在换列的时候是不是写错写成换行了 module matixsol
implicit none
contains
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c matrix solving program
c 2024.1.18
ctype|variable | description
c i|mata| linear coefficient matrix
c i|matb| constant array
c o|matx| element array
c the matrixs must be full rank
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine matsolve(mata,matb,matx)
implicit none
c input/output
real*8 :: mata(:,:)
real*8 :: matb(:)
real*8 :: matx(:)
c temporary variates
real*8,allocatable :: matl(:,:)
real*8,allocatable :: matu(:,:)
real*8,allocatable :: maty(:)
integer*4 :: numb
integer*4 :: cycl
ccccccccccccccccccccccccccccccccccc
numb = size( matb )
allocate( matl(numb,numb) )
allocate( matu(numb,numb) )
allocate( maty(numb) )
c initialize
matl = 0.0
matu = 0.0
maty = 0.0
c LU decomposition of matrix
matu(1 , 1:numb) = mata(1 , 1:numb)
matl(1:numb , 1) = mata(1:numb , 1) / matu(1 , 1)
do cycl = 2 , numb
matu(cycl , cycl:numb) = mata(cycl , cycl:numb)
& - matmul(matl(cycl , 1:cycl-1)
& , matu(1:cycl-1 , cycl:numb))
if (cycl .le. numb) then
matl(cycl:numb , cycl) = (mata(cycl : numb , cycl)
& - matmul(matl(cycl:numb
& , 1:cycl-1) , matu( 1:cycl-1
& , cycl))) / matu(cycl , cycl)
end if
end do
c solve Uy=b
maty(1) = matb(1)
do cycl = 2 , numb
maty(cycl) = matb(cycl) - dot_product(matl(cycl , 1:cycl-1)
& , maty(1:cycl-1))
end do
c solve Ux=y
matx(numb) = maty(numb) / matu(numb , numb)
do cycl = numb-1 , 1 , -1
matx(cycl) = (maty(cycl) - dot_product( matu(cycl , cycl + 1
& : numb),matx(cycl + 1 : numb))) / matu(cycl,cycl)
end do
return
end subroutine matsolve
end module matixsol
这是我之前写的LU的代码，里面也没有换行的操作 感谢老师解答，我的想法是在LU分解中结合列选主元，把每列中最大的元素放到主元所在的位置，以增加数值稳定性
我发现问题所在了，在选列主元的时候，应该同时对矩阵A和矩阵L进行行变换，使得主元位置元素的绝对值最大。
页:
[1]