最小二乘法-Fortran版本.docVIP

! subroutine Lsqnonneg(A,b,x,m,n) !http://QFT !代码文件调用方式： !call Lsqnonneg(A,b,x,m,n) !求解方程组 Ax=b !其中系数矩阵 A 为 m×n 的矩阵， !本算法对于 m n时（即超定方程组）有效。 !mn 时计算结果不一定有效 !***************** ! 用最小二乘法计算方程组 ! 传进来 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(iteritmax) then write(*,*)迭代次数过多，请增大 tol 的值，即增大误差。 write(*,*) tol= ,tol return endif enddo !***! do while !**************************************************** return end subroutine Lsqnonneg !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !*****************************************************! !