【Fortran90】求解定解线性方程组的直接法——矩阵变换

流星雨

感谢[url=【Python】高斯消元法与LU分解 https://www.0xaa55.com/forum.php ... 13&fromuid=6285 (出处: 技术宅的结界)]【Python】高斯消元法与LU分解[/url]及作者watermelon在学习过程中提供的重要帮助！
参考数目：
【1】《数值分析》(李乃成、梅立泉编著，科学出版社，ISBN：978-7-03-032192-3)
【2】《Fortran95/2003科学计算与工程》(宋叶志、茅永兴、赵秀杰编著，清华大学出版社，ISBN：978-7-302-24706-7)

主要内容：
项目的主要功能为求解“定解、显参数的线性方程组”，解法为直接解法(矩阵变换)，共包含5种算法：
1、列主元高斯法：参见【1】P20-22；【2】P12-16
2、LU分解法：参见【1】P25-39；【2】P24-28
3、THOMAS追赶法(只适用于求解三对角阵)：参见【1】P35-37；【2】P28-32
4、Gwens QR分解法：参见【1】P52-56；
5、Householder QR分解法：参见【1】P56-65；【2】P106-112

采用模块化编程，共包含两个程序：Make Cards.exe(生成一张空白的输入卡)；SOLVERS FOR LINEAR EQUATIONS.exe(计算向量x，输出在fout.txt中)
源码力求简明易懂，没有复杂的编程技巧，方便初学者交流学习。作者才疏学浅，有不足之处恳请指教。

附上参考书【1】给出的使用建议：
1、列主元高斯消去法具有较好的数值稳定性，算法组织比较容易，对于阶数不高（n<=30）的线性方程组是一个常用的求解方法；
2、由高斯消去法的矩阵形式导出矩阵的A=LU分解，用LU分解解方程组等价于高斯消去法．为了减少计算量，矩阵A=LU分解特别适宜于求解系数矩阵相同，但右端项不同的多个方程组。(Ax=b→LUx=b→Ly=b，其中L为下三角阵。当b不同时，可先直接解出y，再解Ux=y)；
3、求解三对角方程组的追赶法则是矩阵A=LU分解在解三对角方程组中的应用，它具有计算量小、方法简单、计算过程稳定等优点；
4、求解病态方程组最好事先采取预处理措施，或采用高精度计算，或用矩阵的QR分解求解。由Gwens变换和Householder变换分别导出的矩阵的QR分解常用于求解病态方程组，亦可求解超定方程组、求矩阵特征值等数值代数中的许多其他问题(暂不讨论)．

源码快查：

1、列主元高斯法：
!**************************************************************************************************
! File name: SOVER GAUSSIAN ITER OF COLUMN.f90
! Author: Li Xinyu   
! Version: 1.0        
! Date: 2020/9/29
! Description: This is the sover of GAUSSIAN ITER OF COLUMN
! Others: None
! Function List: None
! Sunroutine List: Gausssolve
! Module list:  UPTRI
! History:     
!    1. Date: 2020/9/29
!        Author: Li Xinyu
!        Modification: Began developing the program
!**************************************************************************************************
MODULE GAUSSIAN_COLUMN
        USE UPTRI
        IMPLICIT NONE
       
        !INTEGER,PARAMETER:: DBL = 8
        CONTAINS
        subroutine Gausssolve(A,b,x,N)
        !*************************************************
        ! Version:1.0     
        ! Date: 2020/9/29
        ! Coded by:        Li Xinyu
        ! Calls: UPTRIsover(Aup,bup,x,N)      
        ! Called By: readtxt(A,b,x,N)
        ! Input: A,b,N      
        ! Output:
        ! Return:x     
        ! Others:      
        !*************************************************
                implicit none
                integer :: i, k !循环变量，i行号，k列号
                integer :: N !矩阵行/列数
                integer :: idmax !主元素编号
                real(kind = 8) :: elmax !主元素绝对值
                real(kind = 8) :: temp !暂存同列临行的系数比
                real(kind = 8) :: A(N,N), b(N), x(N) !定义矩阵A，向量b和待求解向量x
                real(kind = 8) :: Aup(N,N), bup(N)
                real(kind = 8) :: Ab(N,N+1) !增广矩阵
                real(kind = 8) :: vtemp1(N+1), vtemp2(N+1) !用作换行的中间数组
               
                Ab(1:N,1:N) = A !对增广Ab的系数区域赋值为A
                Ab(:,N+1) = b !对增广Ab的尾数区域赋值为b
                !===============================================================
                !列主元消去法核心
                do k=1, N-1 !k作为列编号，外循环扫描行
                        elmax = dabs(Ab(k,k)) !记录数值
                        idmax = k !记录行位置
                        do i=k+1, N !先扫描列，i为行号
                                if (dabs(Ab(i,k)) > elmax) then
                                        elmax = dabs(Ab(i,k)) !获取更大的值
                                        idmax = i
                                end if
                        end do !至此已完成了最大元素的查找，下面交换行顺序
                        vtemp1 = Ab(k,
                        vtemp2 = Ab(idmax ,
                        Ab(k, = vtemp2
                        Ab(idmax , = vtemp1
                        !开始消元
                        do i = k+1, N
                                temp = Ab(i,k)/Ab(k,k)
                                Ab(i, = Ab(i, - temp*Ab(k, !k+1行根据k行消元
                        end do
                end do
                !至此，Ab已转换为上三角阵
                !===============================================================
                Aup(:, = Ab(1:N,1:N) !拆分增广阵为系数阵，储存
                bup( = Ab(:,N+1) !拆分增广阵为尾数阵，储存
                !write (*,*) "新系数矩阵"
                !write (*,*) Aup
                !write (*,*) "新尾数矩阵"
                !write (*,*) bup
                call UPTRIsover(Aup,bup,x,N) !调用回带函数
        end subroutine Gausssolve
END MODULE GAUSSIAN_COLUMN
       
2、LU分解法：
!**************************************************************************************************
! File name: SOVER LU.f90
! Author: Li Xinyu   
! Version: 1.0        
! Date: 2020/9/30
! Description: This is the sover of LU
! Others: None
! Function List: None
! Sunroutine List: LUsolve
! Module list:  UPTRI
! History:     
!    1. Date: 2020/9/30
!        Author: Li Xinyu
!        Modification: Began developing the program
!**************************************************************************************************
MODULE LU
        USE UPTRI
        IMPLICIT NONE
        CONTAINS
        subroutine LUsolve(A,b,x,N)
        !*************************************************
        ! Version:1.0     
        ! Date: 2020/9/30
        ! Coded by:        Li Xinyu
        ! Calls: UPTRIsover(R,y,x,N)
        ! Called By: SOLVERSFORLINEAREQUATIONS
        ! Input: A,b,N      
        ! Output:
        ! Return:x     
        ! Others:      
        !*************************************************
                implicit none
                integer :: i, j, k !循环变量，i行号，k列号
                integer :: N !矩阵行/列数
                real(kind = 8) :: L_U(N,N), A(N,N), b(N), x(N), L1(N,N), U1(N,N),y(N) !定义L_U(N,N)、A、分块下上三角阵L_(N,N)、 _U(N,N)，向量b和待求解向量x,中间向量y Ly=b
                real(kind = 8) :: sumlu = 0 !定义求和变量
                L_U = 1 !赋初值，主要是为了把LII=1
                L_U(1, = A(1,:) !填满U矩阵第一行
                L_U(:,1) = A(:,1)/L_U(1,1) !填满I矩阵第一列
                L_U(1,1) = A(1,1) !这里注意，由于是紧凑型，U1会被L1替代，故需要回调U1
                do i=2, N-1
                        sumlu = 0 !归零
                        do k=1, i-1
                                        sumlu = sumlu + L_U(i,k)*L_U(k,i) !求LiUi之积的和
                        end do
                        L_U(i,i) = A(i,i) - sumlu !求LU对角UII                       
                        do j=i+1, N               
                                sumlu = 0 !归零
                                do k=1, i-1
                                        sumlu = sumlu + L_U(i,k)*L_U(k,j) !求右侧UIJ
                                end do
                                L_U(i,j) = A(i,j) - sumlu                               
                                sumlu = 0 !归零
                                do k=1, i-1
                                        sumlu = sumlu + L_U(j,k)*L_U(k,i) !求下侧LJI
                                end do
                                L_U(j,i) = (A(j,i) - sumlu)/L_U(i,i)                               
                        end do
                end do               
                sumlu = 0 !归零
                do k=1, N-1
                        sumlu = sumlu + L_U(N,k)*L_U(k,N)
                end do
                L_U(N,N) = A(N,N) - sumlu
               
                L1(N,N) = 0 !矩阵分割
                do i=1, N
                        L1(i,i) = 1
                        U1(i,i) = L_U(i,i)
                        do j=i+1, N
                                L1(j,i) = L_U(j,i)
                                U1(i,j) = L_U(i,j)
                        end do
                end do
                !LUx=b
                !Ly=b
                !Ux=y
                y(1) = b(1)/L1(1,1) !先算第一行
                do k=2, N !下三角阵，正循环
                        y(k) = b(k)
                        do i=1, k-1
                                y(k) = y(k) - L1(k,i)*y(i)
                        end do
                        y(k) = y(k)/L1(k,k)
                end do !至此，Ux=y对应Aupx=bup
                call UPTRIsover(U1,y,x,N) !调用回带函数
        end subroutine LUsolve
END MODULE LU

3、THOMAS追赶法(只适用于求解三对角阵)：
!**************************************************************************************************
! File name: SOVER THOMAS.f90
! Author: Li Xinyu   
! Version: 1.0        
! Date: 2020/9/29
! Description: This is the sover of THOMAS
! Others: None
! Function List: None
! Sunroutine List: Thomassolve
! Module list: UPTRI
! History:     
!    1. Date: 2020/9/29
!        Author: Li Xinyu
!        Modification: Began developing the program
!**************************************************************************************************
MODULE THOMAS
        IMPLICIT NONE
        CONTAINS
        subroutine Thomassolve(A,f,x,N)
        !*************************************************
        ! Version:1.0     
        ! Date: 2020/10/1
        ! Coded by:        Li Xinyu
        ! Calls:
        ! Called By: SOLVERSFORLINEAREQUATIONS
        ! Input: A,b,N      
        ! Output:
        ! Return:x     
        ! Others:      
        !*************************************************
                implicit none
                integer :: i, N !循环变量，i行号，k列号
                real(kind = 8) :: A(N,N), f(N), x(N), y(N) !定义矩阵A，向量f和待求解向量x，中间向量y
                real(kind = 8) :: l(2:N), u(N)!向量l(f上分解)，u(f下分解)
                real(kind = 8) :: c(1:N-1), b(N), e(2:N) !定义向量c(记录上对角线元素)，向量b(记录对角线元素)，向量e(记录对角线元素)
                !矩阵对角线赋值给向量e，b，c
                do i=1, N
                        b(i) = A(i,i)
                end do
                do i=1, N-1
                        c(i) = A(i,i+1)
                end do
                do i=2, N
                        e(i) = A(i,i-1)
                end do
                !快速LU分解
                u(1) = b(1)
                do i=2, N
                        l(i) = e(i)/u(i-1)
                        u(i) = b(i)-l(i)*c(i-1)
                end do
                !快速回代求y
                y(1) = f(1)
                do i=2, N
                        y(i) = f(i)-l(i)*y(i-1)
                end do
                !快速回代求x
                x(n) = y(n)/u(n)
                do i=N-1, 1, -1
                        x(i) = (y(i)-c(i)*x(i+1))/u(i)
                end do
        end subroutine Thomassolve
END MODULE THOMAS
       
4、Gwens&Householdersolve QR分解法：
!**************************************************************************************************
! File name: SOVER LU.f90
! Author: Li Xinyu   
! Version: 1.0        
! Date: 2020/10/9
! Description: This is the sover of QR
! Others: None
! Function List: None
! Sunroutine List: Gwenssolve,Householdersolve
! Module list:  
! History:     
!    1. Date: 2020/10/9
!        Author: Li Xinyu
!        Modification: Began developing the program
!**************************************************************************************************
MODULE QR
        USE UPTRI
        IMPLICIT NONE
        CONTAINS
        subroutine Gwenssolve(A,b,x,N)
        !*************************************************
        ! Version:1.0     
        ! Date: 2020/10/10
        ! Coded by:        Li Xinyu
        ! Calls: UPTRIsover(R,y,x,N)
        ! Called By: SOLVERSFORLINEAREQUATIONS
        ! Input: A,b,N      
        ! Output:
        ! Return:x     
        ! Others:      
        !*************************************************
                implicit none
                integer :: i, j, k !循环变量，i行号，j列号
                integer :: N !矩阵列/行数
                real(kind = 8) :: b(N), x(N), y(N) !向量b,待求解向量x,中间向量y
                real(kind = 8) :: A(N,N), PA(N,N,N),PP(N,N,N),R(N,N),Q(N,N) !声明输入矩阵A,PA中间过程矩阵,PP中间过程矩阵,R矩阵,Q矩阵
                real(kind = 8) :: P(N,N,N) !吉文斯矩阵空间P
                !---------------------------------------------------------------------------
                P = 0.0 !P初始化
                PP = 0.0 !Q初始化
                !---------------------------------------------------------------------------
                do j=1, N
                        do i=1, N
                                P(j,i,i) = 1.0 !吉文斯矩阵空间P分配初值，其中第一个i为P的代号，第二三个为横纵坐标
                                PP(j,i,i) = 1.0 !构造单位阵PP
                        end do
                end do
                !---------------------------------------------------------------------------
                PA(1,:,:) = A(:,:) !吉文斯矩阵空间P分配初值，其中第一个i为P的代号，第二三个为横纵坐标\
                do j=1, N-1
                        do i=j+1, N !确定列后每行寻找非0元素并构造吉文斯矩阵P
                                if (abs(PA(j,i,j)) >= 1.0E-5) then !判断非零元素
                                        P(j,:,:) = P(N,:,:) !以最后一个单位阵为基准重新构造单位阵PP
                                        !构造吉文斯矩阵P，这里要注意正负号开根分类讨论。须满足“PijAjj + PiiAij = 0”和“Pij^2 + Pjj^2 = 1”。
                                        if (PA(j,i,j)*PA(j,j,j) > 0) then !Aij和Ajj同号则pij<0
                                                P(j,i,j) = -SQRT(1/(1+(PA(j,j,j)/PA(j,i,j))**2)) !<0
                                                P(j,j,i) = -P(j,i,j) !>0
                                                P(j,j,j) = SQRT(1-P(j,i,j)**2) !>0
                                                P(j,i,i) = P(j,j,j) !>0
                                        else if (PA(j,i,j)*PA(j,j,j) <= 0) then !Aij和Ajj异号则pij>0
                                                P(j,i,j) = SQRT(1/(1+(PA(j,j,j)/PA(j,i,j))**2)) !>0
                                                P(j,j,i) = -P(j,i,j) !<0
                                                P(j,j,j) = SQRT(1-P(j,i,j)**2) !>0
                                                P(j,i,i) = P(j,j,j) !>0
                                        end if
                                        PA(j,:,:) = matmul(P(j,:,:),PA(j,:,:)) !计算矩阵乘法PA
                                        PP(j,:,:) = matmul(P(j,:,:),PP(j,:,:)) !计算累乘P
                                        PA(j+1,:,:) = PA(j,:,:) !每替换一次，就要更新下一层PA
                                        PP(j+1,:,:) = PP(j,:,:)        !每替换一次，就要更新下一层PP
                                end if
                        end do
                end do       
                R = PA(j,:,:) !提取矩阵R
                Q = transpose(PP(j,:,:)) !提取转置Q
                y = matmul(transpose(Q),b) !求解向量y y=Qtb，Rx=y
                call UPTRIsover(R,y,x,N) !调用回带函数
        end subroutine Gwenssolve
       
        subroutine Householdersolve(A,b,x,N)  !镜像变换，H = I - 2 uut/(||u||2)^2
        !*************************************************
        ! Version:1.0     
        ! Date: 2020/9/30
        ! Coded by:        Li Xinyu
        ! Calls: UPTRIsover(R,y,x,N)
        ! Called By: SOLVERSFORLINEAREQUATIONS
        ! Input: A,b,N      
        ! Output:
        ! Return:x     
        ! Others:      
        !*************************************************
                implicit none
                integer :: i, j, k !循环变量，i行号，j列号/批次，k为每批内的列号
                integer :: N !矩阵行/列数
                real(kind = 8) :: b(N), x(N), y(N), u(N) !向量b,待求解向量x,中间向量y,u
                real(kind = 8) :: A(N,N), A1(N,N), A2(N,N), R(N,N), Q(N,N) !声明输入矩阵A,A1中间过程矩阵,A2中间过程矩阵,R矩阵,Q矩阵
                real(kind = 8) :: H0(N,N), H1(N,N), H2(N,N) !声明大矩阵H0，小矩阵H1
                real(kind = 8) :: s, du !声明sigma记录二范数，du取u长度平方
                !---------------------------------------------------------------------------
                A1 = A !中间过程，保存HA
                !---------------------------------------------------------------------------
                H1 = 0.0 !小矩阵单位化处理
                do i=1, N
                        H1(i,i) = 1.0
                end do
                !---------------------------------------------------------------------------       
                do j=1, N
                        !---------------------------------------------------------------------------
                        H0 = 0.0 !大矩阵每次都要重置
                        do i=1, N !大矩阵单位化处理
                                H0(i,i) = 1.0
                        end do
                        !---------------------------------------------------------------------------
                        s = 0.0 !s每次都要重置
                        do i=j, N
                                s = s + A1(i,j)**2
                        end do
                        s = dsqrt(s) !计算向量范数
                        !---------------------------------------------------------------------------
                        u = 0.0 !每一列的u都要初始化
                        if (A1(j,j)>=0) then !取得尽可能大的u的范数
                                u(j) = A1(j,j) + s
                        else
                                u(j) = A1(j,j) - s
                        end if
                        !---------------------------------------------------------------------------
                        do i=j+1, N !取u的值用来取长度平方
                                u(i) = A1(i,j)
                        end do
                        du = 0.0 !每一列的du都要初始化
                        do i=j, N
                                du = du + u(i)**2
                        end do
                        !---------------------------------------------------------------------------
                        do i=j, N !镜像变换，H = I - 2 uut/(||u||2)^2
                                do k=j, N
                                        H0(i,k) = -2.0*u(i)*u(k)/du
                                        if (i==k) then
                                                H0(i,k) = 1.0 + H0(i,k) !直接加常数1
                                        end if
                                end do
                        end do
                        !---------------------------------------------------------------------------
                        A2 = matmul(H0,A1)
                        A1 = A2 !计算HA累乘
                        H1 = matmul(H1,H0) !计算HH累乘       
                end do
                !---------------------------------------------------------------------------
                R = A1 !提取矩阵R
                Q = H1 !提取转置Q
                y = matmul(transpose(Q),b) !求解向量y y=Qtb，Rx=y
                call UPTRIsover(R,y,x,N) !启动上三角回代
        end subroutine Householdersolve
END MODULE QR

5、上三角回代：
!**************************************************************************************************
! File name: UPTRI.f90
! Author: Li Xinyu, Zhou Xingguang        
! Version: 1.0        
! Date: 2020/9/29
! Description: To calculate x in up triangle
! Others: None
! Function List: None
! Sunroutine List: Gausssolve
! Module list:  None
! History:     
!    1. Date: 2020/9/29
!        Author: Li Xinyu
!        Modification: Began developing the program
!**************************************************************************************************
MODULE UPTRI
        IMPLICIT NONE
        CONTAINS
        subroutine UPTRIsover(A,b,x,N)
        !*************************************************
        ! Version:1.0     
        ! Date: 2020/9/29
        ! Coded by:        Li Xinyu
        ! Calls: None      
        ! Called By: Gausssolve(A,b,x,N),LUsolve(A,b,x,N)   
        ! Input: A,b,N      
        ! Output: x
        ! Return: x     
        ! Others:      
        !*************************************************
                implicit none
                integer i,j,N !i行号，j列号，N矩阵行、列数
                real(kind = 8) :: A(N,N), b(N), x(N) !定义矩阵A，向量b和待求解向量x
                !下面开始回代
                x(N) = b(N)/A(N,N) !先算最后一行
                do i=N-1, 1, -1 !i做减法反循环
                        x(i) = b(i)
                        do j=i+1, N !一列列的减完
                                x(i) = x(i) - x(j)*A(i,j)
                        end do
                        x(i) = x(i)/A(i,i) !最后除以当前行x(i)的系数
                end do
        end subroutine UPTRIsover
END MODULE UPTRI

6、输入模块：
!**************************************************************************************************
! File name: INPUT.f90
! Authori Xinyu        
! Version:1.0        
! Date:2020/9/29
! Description: This is the module of reading equations from .txt and transfer them to solvers
! Others:  
! Function List:  
! History:     
!    1. Date: 2020/9/29
!        Author: Li Xinyu
!        Modification: Began developing the program
!    2.        Date: 2020/9/30
!        Author: Li Xinyu
!        Modification: Change the rule of reading to read new cards.
!**************************************************************************************************
MODULE INPUT
        IMPLICIT NONE
        INTEGER,PARAMETER :: fileID = 10
        CONTAINS
       
        subroutine readtxt(A,b,x,N,methods)
        !*************************************************
        ! Version:1.0     
        ! Date: 2020/10/9
        ! Coded by:        Li Xinyu
        ! Calls:
        ! Called By: SOLVERSFORLINEAREQUATIONS(main program)
        ! Input: A,b,x,N,methods     
        ! Output: A,b,x,N,methods
        ! Return:  
        ! Others:      
        !*************************************************
                integer i,j,N,methods !读取用循环变量
                !real(kind = 8) Na(1)
                real(kind = 8),allocatable :: A(:,:), b(:), x(:) !定义矩阵A，向量b和待求解向量x
                character test !测试用字符
                open (unit = fileID,file = 'fin.txt')
                !---------------------------------------------------------------------------
                read (fileID,*) test !读取#用作判定，读取列数
                if (test == '#') then
                        read (fileID,'(1I3)') N
                        write (*,*) '矩阵行列数：'
                        write (*,*) N
                else
                        write (*,*) '输入格式错误：'
                end if
                !---------------------------------------------------------------------------
                allocate(A(N,N)) !分配空间
                allocate(b(N))
                allocate(x(N))
                !---------------------------------------------------------------------------
                write (*,*) !空一行
                read (fileID,*) test!跳至下一行
                if (test == '#') then
                        read (fileID,*) ((A(i,j),j=1,N),i=1,N) !读取A矩阵。这里A是正常写法，故按照先读j再读i的顺序(先控制行不变，读列)，与内存排布方式相反。如果为了加快速度，写矩阵的时候就要转置写法
                        write (*,*) "系数矩阵"
                        write (*,*) A
                else
                        write (*,*) '输入格式错误：'
                end if
                !---------------------------------------------------------------------------
                write (*,*) !空一行
                read (fileID,*) test!跳至下一行
                if (test == '#') then
                        read (fileID,*) b !读取b转置尾数向量
                        write (*,*) "尾数矩阵："
                        do i=1,N
                                write (*,'(D15.3)') b(i)
                        end do
                else
                        write (*,*) '输入格式错误：'
                end if
                !---------------------------------------------------------------------------
                write (*,*) !空一行
                read (fileID,*) test!跳至下一行
                if (test == '#') then
                        read (fileID,*) methods !读取解法
                        write (*,*) "解法："
                        if (methods == 1) then !选择解法，1=列主元高斯
                                write (*,*) "列主元高斯"
                        else if (methods == 2) then !选择解法，2=LU
                                write (*,*) "LU"
                        else if (methods == 3) then !选择解法，3=THOMAS
                                write (*,*) "THOMAS"
                        else if (methods == 4) then !选择解法，4=Gwens QR
                                write (*,*) "Gwens QR"
                        else if (methods == 5) then !选择解法，5=Householder QR
                                write (*,*) "Householder QR"
                        end if
                        write (*,*) !空一行
                else
                        write (*,*) '输入格式错误：'
                end if
                !---------------------------------------------------------------------------
                close (fileID) !及时关闭fin.txt
                !return
                !call Gausssolve(A,b,x,N)  !调用回带函数
        end subroutine readtxt
END MODULE INPUT       

7、求解程序
!**************************************************************************************************
! File name: SOLVERS FOR LINEAR EQUATIONS.f90
! Authori Xinyu        
! Version:1.0        
! Date:2020/9/29
! Description: !!The main program!! To read equations, choose solvers, solve, and output results
! Others:  This program        is also an interface to        solutions for linear equations as much as possible
! Function List: None
! Sunroutine List: CHOOSE SOVERS
! Module list: None
!        INPUT('INPUT.f90'). Reading equations and transfer them to solvers.
!        GAUSSIAN COLUMN('SOVER GAUSSIAN ITER OF COLUMN.f90'). Sover of GAUSSIAN ITER OF COLUMN.
!        LU(SOVER LU.f90). Sover of LU
!        THOMAS(SOVER THOMAS.f90) Sover of THOMAS
!        UPTRI(UPTRI.f90). To calculate x in        up triangle
! History:     
!    1. Date: 2020/9/29
!        Author: Li Xinyu
!        Modification: Began developing the program
!    2. Date: 2020/9/30
!        Author: Li Xinyu
!        Modification: Developed the MCARDS module to make readable cards for users. Developed the GAUSSIAN COLUMN module and LU module
!         3. Date: 2020/10/1
!        Author: Li Xinyu
!        Modification: Developed the THOMAS module
!         4. Date: 2020/10/10
!        Author: Li Xinyu
!        Modification: Developed the QR module
!**************************************************************************************************
program SOLVERSFORLINEAREQUATIONS
        use INPUT !读卡模块
        use GAUSSIAN_COLUMN !列主元高斯迭代法模块
        use LU !LU分解模块
        use THOMAS !THOMAS追赶法模块
        use QR !QR分解模块
        use UPTRI !回代模块
        implicit none
        INTEGER,PARAMETER:: fileIDout = 20
        INTEGER :: N = 1 !列数初值       
        INTEGER :: M = 1 !行数初值       
        INTEGER :: methods
        !---------------------------------------------------------------------------
        real(kind = 8),allocatable :: A(:,:), b(:), x(:) !定义矩阵A，向量b和待求解向量x
        call readtxt(A,b,x,N,methods) !调用读取函数
        !---------------------------------------------------------------------------
        open (unit = fileIDout,file = 'fout.txt') !输出
        write (fileIDout,*) "SOLVERS FOR LINEAR EQUATIONS--Matrix Transformation"  
        !---------------------------------------------------------------------------
        if (methods == 1) then !选择解法，1=列主元高斯
                call Gausssolve(A,b,x,N)  !启用列主元高斯法，自动启动回代计算
                write (fileIDout,*) "列主元高斯法"
        else if (methods == 2) then !选择解法，2=LU
                call LUsolve(A,b,x,N)   !启用LU法，自动启动回代计算
                write (fileIDout,*) "LU法"
        else if (methods == 3) then !选择解法，3=THOMAS
                call Thomassolve(A,b,x,N)   !启用THOMAS法，自动启动回代计算
                write (fileIDout,*) "THOMAS法"
        else if (methods == 4) then !选择解法，4=Gwens QR
                call Gwenssolve(A,b,x,N) !启用Gwens QR法，自动启动回代计算
                write (fileIDout,*) "Gwens QR法"
        else if (methods == 5) then !选择解法，5=Householder QR
                call Gwenssolve(A,b,x,N) !启用Householder QR法，自动启动回代计算
                write (fileIDout,*) "Householder QR法"
        end if
        !---------------------------------------------------------------------------
        !下面开始输出
        write (*,*) "x的解为："
        write (*,'(D15.3)') x
        write (fileIDout,'(D15.3)') x
        close (fileIDout) !及时关闭fin.txt
        !---------------------------------------------------------------------------
        deallocate(A)
        deallocate(b)
        deallocate(x)
        !---------------------------------------------------------------------------
        pause
        !---------------------------------------------------------------------------
end program SOLVERSFORLINEAREQUATIONS
       
8、制卡程序
!**************************************************************************************************
! File name: SOLVERS FOR LINEAR EQUATIONS.f90
! Authori Xinyu        
! Version:1.0        
! Date:2020/9/30
! Description: !!The main program!! To make cards and solve methods
! Function List: None
! Sunroutine List: None
! Module list: None
! History:  
!    1. Date: 2020/10/10
!        Author: Li Xinyu
!        Modification: Developed the program to make readable cards for users
!**************************************************************************************************       
        PROGRAM MAKETXT
                IMPLICIT NONE
                INTEGER,PARAMETER:: fileIDmake = 11
               
                open (unit = fileIDmake,file = 'fin.txt') !打开文件
               
                write (fileIDmake,*) "#请在下一行输入系数矩阵行列数："
                write (fileIDmake,*) "#请在下一行输入系数矩阵，每行元素空格隔开(请勿转置)："
                write (fileIDmake,*) "#请在下一行输入尾数转置向量，每行元素空格隔开："
                write (fileIDmake,*) "#请在下一行选择解法，1=列主元高斯；2=LU；3=THOMAS；4=Gwens QR；5=Householder QR"
               
                close (fileIDmake)
               
        end PROGRAM
       

watermelon

宇哥强啊！如果前面部分（背景介绍，QR分解作用等等）能写的更完善些就更好了！
模块化的编程思维还是很不错的。