DFP.rar_ fortran DFP_dfp

!!!本程序适用于求解形如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))&lt;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&gt;=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&lt;f1)then 4 dx=dx+dx x2=x0-dx f2=f(x+x2*d,A,b) if(f2&lt;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&gt;=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)&lt;=tol)then x0=(a1+b1)/2 else if(f1&gt;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))&lt;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&gt;=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&lt;f1)then 4 dx=dx+dx x2=x0-dx f2=f(x+x2*d,A,b) if(f2&lt;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&gt;=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)&lt;=tol)then x0=(a1+b1)/2 else if(f1&gt;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