最小二乘法-Fortran版本

!
subroutine Lsqnonneg(A,b,x,m,n)
!
!代码文件调用方式：
!call Lsqnonneg(A,b,x,m,n)
!求解方程组 Ax=b
!其中系数矩阵 A 为 m×n 的矩阵，
!本算法对于 m> n时（即超定方程组）有效。

!m<n 时计算结果不一定有效
!*****************
! 用最小二乘法计算方程组
! 传进来 Amatrix Bvector
use constants
use vars_fangchengzu_shuju
!use qiujie_jieguo
implicit none
integer, intent(in) :: m,n
logical, allocatable :: p(:),r(:)
real, allocatable :: w(:),s(:)
real,intent(in),dimension(m,n) :: A real,intent(in),dimension(m) :: b real,intent(inout),dimension(n) :: x
real, allocatable :: AT(:,:),bAx(:)
! transpose of A
! b-Ax
real, allocatable :: norm_lie(:)
! 放矩阵A中，n列数据的绝对值的和
real :: tol
real :: arg_min
! arg_mig=||A^T(b-Ax)||
real :: max_wi
integer :: i,jc,j(1),iter
integer :: itmax
iter=0;
!ishape=shape(Amatrix)
!m=ishape(1);n=ishape(2)
allocate(p(n),r(n))
p=.false.
r=.true.
allocate(w(n),s(n))
w=0.;s=0.
allocate(at(n,m),bAx(m))
!A=Amatrix;b=Bvector;
at=transpose(A);
!deallocate(Amatrix,Bvector)
x=0.
!********************
! 获得误差的判定值
allocate(norm_lie(n))
norm_lie=0.
do i=1,n
do jc=1,m
norm_lie(i)=norm_lie(i)+abs(A(jc,i));
enddo
enddo
tol=1.0e-9*maxval(norm_lie)*m
itmax=3*n;
! 最大的迭代次数。

!******************
bAx=b-matmul(A,x)
w=matmul(at,bAx)
!**************************************************** do while (any(R) .and. maxval(w)>tol) !***! do while j=MAXLOC(w)
! j=argmax(wi)
p(j(1))=.true.;
r(j(1))=.false.
!******************
iter=iter+1;
call lsq_sp(A,b,x,p,r,m,n,tol)
! 返回了新的 x p, r
!******************
!******************
bAx=b-matmul(A,x)
w=matmul(at,bAx)
!*******************
write(*,*) "iter=",iter
write(*,*) "max_w=",maxval(w)
if(iter>itmax) then
write(*,*)"迭代次数过多，请增大 tol 的值，即增大误差。

"
write(*,*)" tol= ",tol
return
endif
enddo !***! do while
!****************************************************
return
end subroutine Lsqnonneg
!*****************************************************!
subroutine Lsq_sp(A,b,x,p,r,m,n,tol)
! Lsqnonneg 的子程序
!use qiujie_jieguo !x(n）
use imsl
implicit none
real :: alpha
integer :: m,n,i,num_pt
real,intent(in) :: tol;
integer :: j(1);
real,intent(in), dimension(m,n) :: A
real,intent(in), dimension(m) :: b
real,intent(inout), dimension(n) :: x
logical,intent(inout), dimension(n) :: p ,r
real, allocatable :: at(:,:)
real, allocatable :: s(:)
real, allocatable :: sp(:),ap(:,:),apt(:,:),xs(:)!,inv_aa(:,:) ! 其中 xs的维数和s的第二维，也就是p为 true 的个数一致。

! xs 作为记录 x/(x-s) 的数组，然后取最大值并放入 alpha
real,allocatable :: w(:)
allocate(w(n),at(n,m))
w=0.;
allocate(s(n))
s=0.
!************
num_pt=0;
do i=1,n
if(p(i)==.true.) then
num_pt=num_pt+1;
endif
enddo
allocate(sp(num_pt),ap(m,num_pt),apt(num_pt,m), xs(num_pt))
sp=0.;ap=0.;apt=0.;xs=0.;!inv_aa=0.
!************
! 下面计算 sp
num_pt=0;
do i=1,n
if(p(i)==.true.) then
num_pt=num_pt+1;
ap(:,num_pt)=a(:,i)
endif
enddo
apt=transpose(ap);
!inv_aa=.i. (matmul(apt,ap))
sp=matmul(matmul(.i.(matmul(apt,ap)),apt),b)
deallocate(apt,ap)
!************
num_pt=0;
do i=1,n
if(p(i)==.true.) then
num_pt=num_pt+1;
s(i)=sp(num_pt)
endif
enddo
!************
if (minval(sp)<=0.)then
num_pt=0;
do i=1,n
! write(*,*) "num_pt=",num_pt
if(p(i)==.true.) then
num_pt=num_pt+1
xs(num_pt)=x(i)/(x(i)-s(i))
endif
enddo
alpha=-minval(xs);
x=x+alpha*(s-x);
!******************
!bAx=b-matmul(A,x)
w=0.
w=matmul(at,(b-matmul(A,x)))
if(any(R) .and. maxval(w)>tol) then
j=MAXLOC(w)
p(j(1))=.true.;
r(j(1))=.false.
endif ! 更新 RP
!******************
endif
deallocate(sp,xs)
!************
!************
! 下面再次计算 sp
num_pt=0;
do i=1,n
if(p(i)==.true.) then
num_pt=num_pt+1;
endif
enddo
allocate(sp(num_pt),ap(m,num_pt),apt(num_pt,m), xs(num_pt)) sp=0.;ap=0.;apt=0.;xs=0.;!inv_aa=0.
!************
num_pt=0;
do i=1,n
if(p(i)==.true.) then
num_pt=num_pt+1;
ap(:,num_pt)=a(:,i)
endif
enddo
apt=transpose(ap);
!inv_aa=.i. (matmul(apt,ap))
sp=matmul(matmul(.i.(matmul(apt,ap)),apt),b) !************
num_pt=0;
s=0.
do i=1,n
if(p(i)==.true.) then
num_pt=num_pt+1;
s(i)=sp(num_pt)
endif
enddo
x=s;
deallocate(sp,ap,apt, xs,s)
return
end subroutine Lsq_sp。