fortran 90稳定的牛顿法.c

Fortran 90

Fortran 90稳定的牛顿法  
 

<P>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
    !!!x,x1为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；x0实型变量，存放步长。
    program main
    real,dimension(:),allocatable::x,gradt,dir,b,x1
    real,dimension(:,:),allocatable::hessin
    real::x0
    integer::n ,iter
    print*,'请输入变量的维数'
    read*,n
    allocate(x(n),gradt(n),dir(n),b(n),x1(n))
    allocate(hessin(n,n))
    print*,'请输入初始向量x'
    read*,x
    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
    read*,b
    iter=0
tol=0.000001
  2 gradt=matmul(hessin,x)+b
    dir=(-1)*gradt
    call cholesky(hessin,dir,n)
    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir 
    iter=iter+1
if(iter&gt;10*n)then
   print*,"out"
   goto 101
endif
print*,"第",iter,"次运行结果为","direction",dir ,"step length",x0
print*,x1,"f(x)=",f(x1,hessin,b)
    if(dot_product(x1-x,x1-x)&lt;=tol)then 
       !print*,x1 ,iter 
       goto 101
    else
       x=x1
    goto 2
    endif
    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
  !!!利用矩阵的cholesky分解；求得迭代方向
    subroutine cholesky(G,d,n)
    real,dimension(:),intent(out)::d
    real,dimension(:,:),intent(in)::G
    real::mb,sm,c,tol ,beita,sita
    integer::i,j,r,n
    real::b(n),s(n),y(n),p(n)
    real::L(n,n)
    do i=1,n-1
       b(i)=0
       do j=i+1,n
          b(i)=max(b(i),abs(G(i,j)))
    enddo
    enddo
mb=0
do i=1,n-1
    mb=max(mb,b(i))
enddo
c=0
do i=1,n
    c=max(c,abs(G(i,i)))
enddo
beita=max(sqrt(c),sqrt(mb/n))
tol=0.1
i=1
s(1)=max(tol,sqrt(abs(G(1,1))))
do j=2,n
    s(j)=G(j,1)/s(1)
    enddo
sita=0
do j=2,n
    sita=max(sita,abs(s(j)))
enddo
if(sita&lt;=beita)then
    do j=1,n
       L(j,1)=s(j)
    enddo
else
    L(1,1)=sita*s(1)/beita
    do j=2,n
       L(j,1)=beita*s(j)/sita
    enddo
endif
    i=i+1
100 if(i==n)then
       sm=0
       do r=1,n-1
       sm=sm+L(n,r)*L(n,r)
    enddo
    L(n,n)=max(tol,sqrt(abs(G(n,n)-sm)))
    do j=1,n
       do r=j+1,n
       L(j,r)=0
    enddo
    enddo
    goto 10
    
else
     sm=0
  do r=1,i-1
     sm=sm+L(i,r)*L(i,r)
  enddo
  s(i)=max(tol,sqrt(abs(G(i,i)-sm)))
        do j=i+1,n
     sm=0
        do r=1,i-1
        sm=sm+L(j,r)*L(i,r)
     enddo
     s(j)=(G(j,i)-sm)/s(i)
   enddo  
   sita=0
   do j=i+1,n
      sita=max(sita,abs(s(j)))
   enddo
endif
if(sita&lt;=beita)then
      do j=i,n
      L(j,i)=s(j)
   enddo
else
      L(i,i)=sita*s(i)/beita
   do j=i+1,n
      L(j,i)=beita*s(j)/sita
   enddo
endif
i=i+1
goto 100
!!!解方LL'P=g,令L'p=y,先求得y,然后解出p
10  y(1)=-d(1)/L(1,1)
do i=1,n
    sm=0
    do r=1,i-1
       sm=sm+L(i,r)*y(r)
    enddo
    y(i)=-(d(i)+sm)/L(i,i)
enddo
p(n)=y(n)/L(n,n)
do i=1,n-1
    sm=0
    do r=1,i-1
       sm=sm+L(n-r,n-i)*p(n-r)
    enddo
    p(n-i)=(y(n-i)-sm)/L(n-i,n-i)
enddo</P>
<P>    d=p</P>
<P> end subroutine
   
    !!!精确线搜索0.618法子程序,返回值为迭代步长
    function golden(x,d,A,b) result(golden_n)
    real::golden_n
    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,x0
    parameter(r=0.618)
x0=1
    tol=0.0001
    dx=0.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 program main
</P>
<P>本算法由Fortran 90语言编写，在Vistrual Fortran 5上编译通过，本程序由沙沙提供！</P>