[插值拟合] 【多元线性回归方程】算法求解~救命呀呀呀呀哎呀~~

program multiple_linear_regression
    implicit none
    integer :: n,num, i,j
    real, allocatable :: x(:,:),y(:),xtrans(:,:),xx(:,:),aa(:,:), b(:)
    write(*,*) "请输入自变量的个数："
    read(*,*) n
    write(*,*) "请输入次数："
    read(*,*) num
    allocate( x(num,n+1) )
    allocate( xtrans(n+1,num) )
    allocate( xx(num,num) )
    allocate( aa(num,num) )!aa是用来存放求逆后的矩阵
    allocate( y(num) )
    allocate( b(n+1) )
    x(:,1)=1
    write(*,*) "请输入自变量X矩阵："
    do i=1, num
        read(*,*) (x(i,j),j=2,n+1)
    end do
    write(*,*) "请输入因变量y的值："
    read(*,*) y(:)
    !转置矩阵
    xtrans=transpose(x)
    !矩阵相乘
    xx=matmul(x,xtrans)
    !（X'X）矩阵求逆
    call nijuzhen(xx,num,aa)
    !求B=（X'X）的逆矩阵 * X'* Y
    b=matmul( matmul(aa,xtrans),y )
    write(*,*) b
    
    end
    
                  
    subroutine nijuzhen(a,n,aa)!高斯消去法求逆矩阵
        implicit none
        integer :: i,k,j,n
        real :: ae(n,n+n),a(n,n),aa(n,n)!a矩阵右侧添加一个单位矩阵
        ae(:,:)=0!!!!
        do i=1,n
            ae(i,1:n)=a(i,:)
            ae(i,n+i)=1
        end do
        !ae化为上三角矩阵
        do k=2,n
            do i=k,n
            ae(i,:)=ae(i,:)-ae(k-1,:)*ae(i,k-1)/ae(k-1,k-1)
            end do
        end do
        !将第i行第i列元素化成1
        do i=1,n
            a(i,:)=a(i,:)/a(i,i)
        end do
        !矩阵变换将左侧化为单位矩阵
        do j=2,n
            do i=1,j-1
                if( a(i,j)/=0 ) then
                    a(i,:)=a(i,:)-a(i,j)*a(j,:)
                end if
            end do
        end do
        !矩阵右侧即得到逆矩阵
        do i=1,n
            aa(i,:)=ae(i,(n+1):(n+n))
        end do
    end