Doolittle分解的算法实现问题

Yjc

各位老师好，我在写用“Doolittle分解法解线性方程组”这一数值算法的实现程序时遇到了问题，程序给出了错误的结果，想请教各位老师这个程序在算法上有什么错误，谢谢！
（附件中是实现doolittle分解法的子程序，这是我写的交互程序的一部分，A,X,B分别代表系数矩阵，解，常数矩阵；ndim代表方阵A的维度）
Doolittle_Decomposition.f90 (2.01 KB, 下载次数: 2)

2024-3-30 11:40 上传

点击文件名下载附件
源代码

1468592801

你在换列的时候是不是写错写成换行了

1468592801

[Fortran] 纯文本查看 复制代码

      module matixsol
      implicit none
      contains
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c matrix solving program
c 2024.1.18
c  type|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的代码，里面也没有换行的操作

Yjc

感谢老师解答，我的想法是在LU分解中结合列选主元，把每列中最大的元素放到主元所在的位置，以增加数值稳定性

Yjc

我发现问题所在了，在选列主元的时候，应该同时对矩阵A和矩阵L进行行变换，使得主元位置元素的绝对值最大。