[求助]求fortran版拟newton BFGS程序
想找fortran版的拟newton BFGS程序研究一下，<BR>不知坛里哪位可以提供一下，<BR>非常感谢！<BR><BR>或者，详细的程序算法及流程。<BR><BR>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]