DFP算法
<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,p ,gradt1,x1
real,dimension(:,:),allocatable::hessin ,H ,G ,U
real::x0,tol
integer::n ,iter,i,j
print*,'请输入变量的维数'
read*,n
allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),p(n),gradt1(n),x1(n))
allocate(hessin(n,n),H(n,n),G(n,n),U(n,n))
print*,'请输入初始向量x'
read*,x
print*,'请输入hessin矩阵'
read*,hessin
print*,'请输入矩阵b'
read*,b
iter=0
tol=0.000001</P>
<P> do i=1,n
do j=1,n
if (i==j)then
H(i,j)=1
else
H(i,j)=0
endif
enddo
enddo
100 gradt=matmul(hessin,x)+b
if(sqrt(dot_product(gradt,gradt))<tol)then
!print*,'极小值点为：',x
!print*,'迭代次数:',iter
goto 101
endif
dir=matmul(H,gradt)
x0=golden(x,dir,hessin,b)
x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
s=x1-x
y=gradt1-gradt
call vectorm(s,G)
U=G
call vectorm(matmul(H,y),G)
H=H+1/dot_product(s,y)*U-1/dot_product(matmul(H,y),y)*G
x=x1
iter=iter+1
if(iter>=10*n)then
print*,"out"
goto 101
endif
print*,"第",iter,"次运行结果为", "方向为",dir,"步长",x0
print*,x,"f(x)=",f(x,hessin,b)
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,j)=p(i)*p(j)
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
end function golden
101 end</P>
<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,p ,gradt1,x1
real,dimension(:,:),allocatable::hessin ,H ,G ,U
real::x0,tol
integer::n ,iter,i,j
print*,'请输入变量的维数'
read*,n
allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),p(n),gradt1(n),x1(n))
allocate(hessin(n,n),H(n,n),G(n,n),U(n,n))
print*,'请输入初始向量x'
read*,x
print*,'请输入hessin矩阵'
read*,hessin
print*,'请输入矩阵b'
read*,b
iter=0
tol=0.000001</P>
<P> do i=1,n
do j=1,n
if (i==j)then
H(i,j)=1
else
H(i,j)=0
endif
enddo
enddo
100 gradt=matmul(hessin,x)+b
if(sqrt(dot_product(gradt,gradt))<tol)then
!print*,'极小值点为：',x
!print*,'迭代次数:',iter
goto 101
endif
dir=matmul(H,gradt)
x0=golden(x,dir,hessin,b)
x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
s=x1-x
y=gradt1-gradt
call vectorm(s,G)
U=G
call vectorm(matmul(H,y),G)
H=H+1/dot_product(s,y)*U-1/dot_product(matmul(H,y),y)*G
x=x1
iter=iter+1
if(iter>=10*n)then
print*,"out"
goto 101
endif
print*,"第",iter,"次运行结果为", "方向为",dir,"步长",x0
print*,x,"f(x)=",f(x,hessin,b)
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,j)=p(i)*p(j)
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
end function golden
101 end
</P>
<P>本程序由Fortran 90编写，在Visual Fortran 5上编译通过，本程序由沙沙提供！
</P> 这是什么啊
页:
[1]