高斯消元法 (2)---主元消除法

最新推荐文章于 2022-12-22 18:09:28 发布

ROCK_LOTUS

最新推荐文章于 2022-12-22 18:09:28 发布

阅读量508

点赞数

文章标签：线性代数矩阵

本文链接：https://blog.csdn.net/ROCK_LOTUS/article/details/105158347

版权

采用主元消除的方法，求解矩阵方程

GaussianElimation.f90

! GaussianElimation.f90
!-------------------------------------------------------GaussianElimation.f90
!----------------------------------------------------------------------------
!------------------该模块包含了高斯消元的主要模块---------------------------------
!--------

!引用迭代求解模块
include "Tri_eq.f90"

!m_gauss模块，需要调用tri_eq模块进行求解
module m_gauss
    use tri_eq

    contains

!solve:求解器
!-----------------------------------------------------------------------
!   输入：
!        A:系数矩阵
!        b:右向量
!        x:待求解的结果
!        N:维数

    subroutine solve(A,b,x,N)
        
        implicit real*8(a-z)
        integer::N 
        real*8::A(N,N),b(N),x(N)

        !增广矩阵Ab
        real*8::Ab(N+1,N)

        !增广矩阵赋值
        Ab(1:N,:) = A
        Ab(N+1,:) = b
        
        !求解行最简式
        call RowReduce(Ab,N)
        
        !将增广矩阵中对应的系数矩阵与右向量重新赋值给A,b
        A = Ab(1:N,:)
        b = Ab(N+1,:)

        !调用上三角求解
        call uptri(A,b,x,N)

    end subroutine

!RowReduce:化简为行最简形式
!   输入
!       Ab:增广矩阵
!       N：维数
    subroutine RowReduce(Ab,N)
        integer::i,j,N,md       !md用于标记某一列中最大值对应的行数，以便于进行交换
        real*8::Ab(N+1,N),temp

        do i = 1,N-1
            md = Max_id(Ab(i,i:N),N-i+1)+(i-1)
            if(md == 0) then    !如果某一列对应的值全为零，则跳入下一步
                cycle
            end if
            !交换最大值对应行和当前行
            call swap(Ab,N,i,md)
            !call DisplayArray(Ab,N)

            !用当前行，对其余所有行进行消元处理
            do j = i+1,N 
                temp = -Ab(i,j)/Ab(i,i)
                call Cijk(Ab,N,i,j,temp)
            end do
        end do

        !此处便于获得对角矩阵，便于求解
        !通过反带，消除上三角系数
        do j = N,2,-1
            do i = j-1,1,-1
                temp = -Ab(j,i)/Ab(j,j)
                call Cijk(Ab,N,j,i,temp)
            end do
        end do

    end subroutine

    
    !交换两行
    subroutine swap(A,N,i,j)
        implicit real*8 (a-z)
        integer::N,i,j
        real*8::A(N+1,N),temp(N+1)

        temp = A(:,i)
        A(:,i) = A(:,j)
        A(:,j) = temp
        
        return
    end subroutine

!将A的第i行乘上系数k加到第j行
    subroutine Cijk(A,N,i,j,k)
        implicit real*8 (a-z)
        integer::N,i,j
        real*8::k 
        real*8::A(N+1,N),temp(N+1)

        temp = A(:,i)*k
        A(:,j) = A(:,j) + temp
    end subroutine

!用于展示数组
    subroutine DisplayArray(A,N)
        implicit real*8(a-z)
        integer::N,i,j
        real*8 A(N+1,N)

        do i = 1,N 
            do j = 1,N+1
                write(*,'(T1,F12.8,$)') A(j,i)
            end do
            write(*,*)
        end do
    end subroutine

!一个函数方法，用于求解某一行中最大的数对应的坐标
    integer function Max_id(b,N)
        integer::i,N
        real*8::Max_Value,b(N)
        Max_id = 1
        Max_Value = abs(b(Max_id))
        do i = 2,N 
            if(Max_Value < abs(b(i))) then
                Max_id = i
            end if
        end do
        return
    end function
    
end module

Tri_eq.f90

module tri_eq
!---------------------------------------------------module comment
!   Version
!   Coded By 
!   Date:2020/03/16
!---------------------------------------------------Description
!   用于解上三角方程回带方法的模块
!---------------------------------------------------Contains
!   contains
!   1. uptri  方法函数
!   2. downtri 方法函数
!   
contains

subroutine uptri(A,b,x,N)
!---------------------------------------------------subroutine comment
!   Version
!   Coded By 
!   Date:2020/03/16
!----------------------------------------------------Description
!   purpose: 上三角回带方法
!           A x = b
!----------------------------------------------------Input
!   Input parameters:
!       A(N,N):系数矩阵
!       b(N)右向量
!       N:方程维数
!   Output parameters:
!       x 方程的根
!   Common parameters:
!
!-----------------------------------------------------
    implicit real(kind = 8) (a-z)
    integer::i,j,N 
    real*8 A(N,N),b(N),x(N)
    x(N) = b(N)/A(N,N)
    !回带
    do i = n-1,1,-1
        x(i) = b(i)
        do j  = i+1,N 
            x(i) = x(i) - A(j,i)*x(j)
        end do
        x(i) = x(i)/A(i,i)
    end do

    end subroutine uptri

subroutine downtri(A,b,x,N)
!---------------------------------------------------subroutine comment
!   Version
!   Coded By 
!   Date:2020/03/16
!----------------------------------------------------Description
!   purpose: 下三角回带方法
!           A x = b
!----------------------------------------------------Input
!   Input parameters:
!       A(N,N):系数矩阵
!       b(N)右向量
!       N:方程维数
!   Output parameters:
!       x 方程的根
!   Common parameters:
!
!-----------------------------------------------------
    implicit real(kind = 8) (a-z)
    integer::i,j,N 
    real*8 A(N,N),b(N),x(N)
    x(1) = b(1)/A(1,1)
    !回带
    do i = 2,N
        x(i) = b(i)
        do j  = 1,i-1 
            x(i) = x(i) - A(j,i)*x(j)
        end do
        x(i) = x(i)/A(i,i)
    end do

    end subroutine downtri

end module tri_eq

ReadFile.f90

module readFile
!---------------------------------------------------module comment
!   Version
!   Coded By 
!   Date:2020/03/16
!---------------------------------------------------Description
!   用于解上三角方程文件读写的模块
!---------------------------------------------------Contains
!   contains
!   1. readParameter 方法函数
!   2. writeFile 方法函数
!   
    contains

    subroutine readParameter(FileName,A,b,N,status1,status2)
        !Description:用于读取输入
        !variable declaration
        !FileName:输入文件的名字
        !msg:错误信息
        !N:维数
        !status1:打开文件状态，status2:读文件状态
        !A(N,N),b(N):系数矩阵
        character(len = 30)::FileName
        character(len = 30)::msg
        integer::N
        integer::status1,status2
        real(kind = 8) A(N,N),b(N)

        !openFile
        open(unit = 10, file = FileName, status = 'OLD', action = 'READ', IOSTAT = status1, IOMSG = msg)

        if(status1 == 0) then 
            !如果打开文件成功

            !read the parameter
            read(10,*,IOSTAT = status2) ((A(i,j),i = 1,N),j = 1,N)
            read(10,*,IOSTAT = status2) b

            !status2 /= 0，读取文件失败或达到文件末尾
            if(status2 /= 0) then
                write(*,*) "File reading Error or End"
                close(10)
                return
            end if

        end if

    end subroutine

    subroutine writeFile(FileName,X,N)   
        !Description:用于写入结果
        !variable declaration
        !FileName:写入文件的名字
        !msg:错误信息
        !N:维数
        !X：结果
 
        character(len = 30) FileName
        character(len = 30) msg 
        integer::status1, status2
        integer::N
        real(kind = 8)::X(N)
       
        !打开文件，此处使用status = 'OLD'，表示该文件已经存在
        open(11,file = FileName,status = 'OLD',action = 'WRITE',IOSTAT = status1, IOMSG = msg)

        if(status1 == 0) then
            write(11,101,IOSTAT = status2) X
            101 format(T5,"方程组的解",/,T4,"x = ",6(/F12.8))
        end if

    end subroutine

end module

main.f90

include"GaussianElimation.f90"
include"ReadFile.f90"

program main
    use m_gauss
    use readFile

    !declaration
    integer,parameter::N = 6
    integer::status1,status2        !分别用于表示文件打开和读取的状态

    real*8::A(N+1,N),Aup(N,N),bup(N),x(N)

    character(len = 30)::inputFile,outputFile,msg

    inputFile = "fin.txt"
    outputFile = "fout.txt"
    open(11,file = outputFile,status = 'NEW',action = 'WRITE',IOSTAT = status1, IOMSG = msg)

    call readParameter(inputFile,Aup,bup,N,status1,status2)
    
    A(1:N,:) = Aup
    A(1+N,:) = bup

    write(*,10)
    10 format("Original Matrix")
    call DisplayArray(A,N)

    write(*,20)
    20 format("Row Reduce")
    call RowReduce(A,N)
    call DisplayArray(A,N)

    Aup = A(1:N,:)
    bup = A(1+N,:)

    write(*,"(/)")

    write(*,30)
    30 format("Solve Equation Ax = b")

    call solve(Aup,bup,x,N)

    write(*,40) x
    40 format(T5,'x = ',/,T4,6(/F12.3))
    call writeFile(outputFile,x,N)
    
end program main

输入文件fin.txt

-3.3435 -0.4946 0.3834  -4.9537 -2.4013 -3.5446
1.0198  -4.1618 4.9613  2.7491  3.0007  -3.6393
-2.3703 -2.7102 -4.2182 3.1730  -0.6859 3.6929
1.5408  4.1334  -0.5732 3.6869  4.1065  0.7970
1.8921  -3.4762 -3.9335 -4.1556 -3.1815 0.4986
2.4815  3.2582  4.6190  -1.0022 -2.3620 -3.5505

9.0537
0.0228 
-8.4177
-4.6380
10.5575
9.8252

输出文件 fout.txt

    方程组的解
   x = 
  1.26515941
  0.11026391
 -1.24520947
 -0.43379925
 -0.99093338
 -2.62011810