Fortran 高斯消元法

最新推荐文章于 2023-05-19 09:08:14 发布

原创最新推荐文章于 2023-05-19 09:08:14 发布 · 2.9k 阅读

22 ·

CC 4.0 BY-SA版权

文章标签：

#线性代数 #矩阵

本文介绍了一种使用Fortran95进行科学计算的高斯消元法实现，通过模块化编程技巧，详细展示了上三角方程回带与下三角方程回带的方法，适用于解决线性方程组问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

@[TOC] Fortran 高斯消元法

Fortran 高斯消元法

最近是学习使用Fortran 95进行一些科学计算的代码编程，看了宋志叶老师所著《科学计算与工程》，磨练一些Fortran的编程技巧，用模块的方法重新编写了Fortran的Gauss消元法。

// An highlighted block
//TriEq.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

//Gauss.f90
include "Tri_eq.f90"

module Gauss
    use tri_eq

    !---------------------------------------------------subroutine comment
    !   Version
    !   Coded By 
    !   Date:2020/03/16
    !----------------------------------------------------Description
    !   purpose: 高斯消元法
    !           A x = b
    !----------------------------------------------------Input
    !   Input parameters:
    !       A(N,N):系数矩阵
    !       b(N)右向量
    !       Ab(N+1,N)：增广矩阵
    !       N:方程维数
    !   Output parameters:
    !       x 方程的根
    !   Common parameters:
    !
    !-----------------------------------------------------
 
    contains

    subroutine solve(A,b,x,N)

        implicit real*8(a-z)

        integer::i,k,N
        real*8::A(N,N), b(N),x(N),temp
        !Aub：为增广矩阵
        real*8::Aup(N,N),bup(N)
        !增广矩阵
        real*8::Ab(N+1,N)

        Ab(1:N,1:N) = A(1:N,1:N)
        Ab(N+1,1:N) = b(:)

        !Gauss eliminate
        do k = 1,N-1
            do i = k+1,N
                temp = Ab(k,i)/Ab(k,k)
                Ab(:,i) = Ab(:,i) - temp * Ab(:,k)
            end do
        end do
        Aup(:,:) = Ab(1:N,1:N)
        bup(:) = AB(1+N,:)

        call uptri(Aup,bup,x,N)

    end subroutine solve

end module

//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 = ",4(/F12.8))
        end if

    end subroutine

end module

//Main.f90

include "Gauss.f90"
include "ReadFile.f90"

program main
    use Gauss
    use readFile

    implicit real*8(a-z)
    integer,parameter::N = 4
    real*8::A(N,N),b(N),x(N)
    integer::status1,status2

    character(len = 30)::inFile,outFile 
    character(len = 30)::msg

    inFile = "fin.txt"
    outFile = "fout.txt"

    open(11,file = outFile,status = 'NEW',action = 'WRITE',IOSTAT = status1, IOMSG = msg)
    close(11)

    status2 = 0
    do while(status2 == 0)
        call readParameter(inFile,A,b,N,status1,status2)

        if(status2 /= 0) then
            exit
        end if

        !求解方程
        call solve(A,b,x,N)
        call writeFile(outFile,X,N)
    end do
    
end program

运行结果如下：
!fin.txt
-13.9211    21.7540     -14.8743    -7.9025
18.3862     -26.0893    -5.6866     4.4451
-4.1683     3.9325      -33.3169    41.7098
-6.0438     6.7108      -32.9591    -23.3378

136.8721
-126.8849
100.4524
95.7901

1 10 7 2
3 3 7 5
8 8 1 3
9 2 8 4

-5
4 
9
6

fout.txt
    方程组的解
   x = 
 -3.39068043
  2.93402079
 -1.88853397
  0.28436704
    方程组的解
   x = 
  0.78118772
 -0.24432045
 -1.02231965
  1.90912714


参考

 [1]: 科学计算与工程，宋志叶