linqus 发表于 2005-11-5 21:33

[求助]求fortran版拟newton BFGS程序

想找fortran版的拟newton BFGS程序研究一下,<BR>不知坛里哪位可以提供一下,<BR>非常感谢!<BR><BR>或者,详细的程序算法及流程。<BR><BR>

lmj_1111 发表于 2006-7-15 19:22

bfgs

BFGS算法<P>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点;
    !!!输入函数信息,输出函数的稳定点及迭代次数;
    !!!iter整型变量,存放迭代次数;
    !!!x为n维变量,初始值由用户输入;gradt实型变量,存放函数梯度;
    !!!dir实型变量,存放搜索方向;
    program main
    real,dimension(:),allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,:),allocatable::hessin ,B1 ,G,G1
    real::x0,tol
    integer::n ,iter,i,j
    print*,'请输入变量的维数'
    read*,n
    allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),gradt1(n),x1(n))
    allocate(hessin(n,n),B1(n,n),G(n,n),G1(n,n))
    print*,'请输入初始向量x'
    read*,x
    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
    read*,b
    iter=0
tol=0.00001</P>
<P> do i=1,n
    do j=1,n
       if (i==j)then
       B1(i,j)=1
    else
       B1(i,j)=0
    endif
    enddo
enddo   
    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))<tol)then
      !print*,'极小值点为:',x
   !print*,'迭代次数:',iter
   goto 101
    endif
call gaussj(B1,n,(-1)*gradt)
dir=gradt
    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
s=x1-x
y=gradt1-gradt
call vectorm(gradt,G)
G1=G
call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
gradt=gradt1
    iter=iter+1
if(iter>10*n)then
    print*,"out"
    goto 101
endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir
    goto 100
    contains</P>
<P>    !!!子程序,返回函数值   
    function f(x,A,b) result(f_result)
    real,dimension(:),intent(in)::x,b
    real,dimension(:,:),intent(in)::A
    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
    end function f
!!!子程序,矩阵与向量相乘
subroutine vectorm(p,G)
real,dimension(:),intent(in)::p
real,dimension(:,:),intent(out)::G
n=size(p)
do i=1,n
    !do j=1,n
       G(i,:)=p(i)*p
    !enddo
enddo
end subroutine

    !!!精确线搜索0.618法子程序 ,返回步长;
    function golden(x,d,A,b) result(golden_n)
    real::golden_n
    real::x0
    real,dimension(:),intent(in)::x,d
    real,dimension(:),intent(in)::b
    real,dimension(:,:),intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
    parameter(r=0.618)
    tol=0.0001
    dx=0.1
    x0=1
    x1=x0+dx
    f0=f(x+x0*d,A,b)
    f1=f(x+x1*d,A,b)
    if(f0<f1)then
4       dx=dx+dx
      x2=x0-dx
      f2=f(x+x2*d,A,b)
      if(f2<f0)then
         x1=x0
      x0=x2
      f1=f0
      f0=f2
      goto 4
      else
         a1=x2
      b1=x1
      endif
    else
2       dx=dx+dx
      x2=x1+dx
      f2=f(x+x2*d,A,b)
      if(f2>=f1)then
         b1=x2
      a1=x0
      else
         x0=x1
      x1=x2
      f0=f1
      f1=f2
      goto 2
      endif
    endif
    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
    f1=f(x+x1*d,A,b)
    f2=f(x+x2*d,A,b)
3   if(abs(b1-a1)<=tol)then
      x0=(a1+b1)/2
    else
      if(f1>f2)then
      a1=x1
      x1=x2
      f1=f2
      x2=a1+r*(b1-a1)
      f2=f(x+x2*d,A,b)
      goto 3
   else
      b1=x2
      x2=x1
      f2=f1
      x1=a1+(1-r)*(b1-a1)
      f1=f(x+x1*d,A,b)
      goto 3
   endif
    endif
    golden_n=x0
    endfunction golden</P>
<P>
    !!!A为二维数,返回值为的A逆;b为一维数组,返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
    integer n,nmax
    real a(n,n),b(n)
    parameter(nmax=50)
    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
    real big,dum,pivinv
    do j=1,n
       ipiv(j)=0
    enddo
    do i=1,n
       big=0.
       do j=1,n
       if(ipiv(j)/=1)then
          do k=1,n
          if(ipiv(k)==0)then
          if(abs(a(j,k))>=big)then
         big=abs(a(j,k))
         irow=j
         icol=k
       endif
       else if(ipiv(k)>1)then
          pause'singular matrix in gaussj'
       endif
       enddo
    endif
    enddo
    ipiv(icol)=ipiv(icol)+1
    if(irow/=icol)then
       do l=1,n
          dum=a(irow,l)
       a(irow,l)=a(icol,l)
       a(icol,l)=dum
       enddo
       dum=b(irow)
       b(irow)=b(icol)
       b(icol)=dum
    endif
    indxr(i)=irow
    indxc(i)=icol
    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
    do l=1,n
      a(icol,l)=a(icol,l)*pivinv
    enddo
    b(icol)=b(icol)*pivinv
    do ll=1,n
       if(ll/=icol)then
          dum=a(ll,icol)
       a(ll,icol)=0
       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
       b(ll)=b(ll)-b(icol)*dum
       endif
    enddo
    enddo
    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
       do k=1,n
          dum=a(k,indxr(l))
       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
       enddo
    endif
    enddo
    end subroutine gaussj
101 end
</P>
<P>本程序用Fortran 90编写,在Virual Fortran 5上编译通过,本程序由沙沙提供!</P>
页: [1]
查看完整版本: [求助]求fortran版拟newton BFGS程序