fortran 90最速下降法.c

Fortran 90

Fortran 90最速下降法！  
 

<P> !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
    program main
    real,dimension(:),allocatable::x,gradt,dir,b
    real,dimension(:,:),allocatable::hessin
    real::x0,tol
    integer::n ,iter
    print*,'请输入变量的维数'
    read*,n
    allocate (x(n),gradt(n),dir(n),b(n))
    allocate(hessin(n,n))
    print*,'请输入初始点x'
    read*,x
    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入向量b'
    read*,b
    tol=0.000001
    iter=0
100 gradt=matmul(hessin,x)+b
    dir=(-1)*gradt
    if(dot_product(gradt,gradt)&lt;=tol)then
        !print*,'输出函数稳定点',x
     !print*,'输出迭代次数',iter
     goto 101
    else
        x0=golden(x,dir,hessin,b)
        x=x+x0*dir
     iter=iter+1
  if(iter&gt;10*n)then
     print*,"out"
  endif
  print*,"第",iter,"次运行结果为","direction",dir,"step length" ,x0
  print*,x,"f(x)=",f(x,hessin,b)
        goto 100
    endif
    contains</P>
<P>    !!!子程序，返回函数f(x)在x点的函数值; 
    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(A,x),x)+dot_product(b,x)
    end function f</P>
<P>    !!!精确线搜索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</P>
<P>101 end
</P>
<P>本算法由Fortran 90语言编写，在Vistrual Fortran 5上编译通过，本程序由沙沙提供！</P>