BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
2 }& `1 e# D9 T1 M% ^5 I    !!!dir实型变量，存放搜索方向；
2 X5 U8 l( i# A/ f7 @2 o# ~7 M/ ]    program main
  c7 Z9 {5 P9 b+ ]( _5 J    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
6 _& K3 _9 I$ D" F6 L; f: A4 E    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
- T0 J8 [1 o& i+ ]0 |, n0 G. ?" D    real::x0,tol
8 Z/ N, o( G! x1 ~2 d  {    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矩阵'
3 Y' v- {& p7 C% r0 K+ c    read*,hessin
8 G& E$ d# u: t1 o) U( k    print*,'请输入矩阵b'
: m7 {3 E. d! a; l! }) l    read*,b
    iter=0
$ q  C" q' y/ p+ o tol=0.00001</P>
6 j' o7 G7 C5 p% I/ ^<> do i=1,n
    do j=1,n
- t5 x2 |- |5 O) ]' G7 b       if (i==j)then
7 [5 g0 U6 n* j  _7 o, V4 ]       B1(i,j)=1
2 h+ Q7 ~5 p7 q) e    else
       B1(i,j)=0
1 D2 u* Z$ v; W: a% B7 v    endif
    enddo
enddo   
  Z$ F' v/ v0 D& _% a8 i    gradt=matmul(hessin,x)+b
% k( ^$ ?5 Y$ F8 G; z9 `100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
! {4 A% M& d) j        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
     goto 101
    endif
call gaussj(B1,n,(-1)*gradt)
; c8 v9 U7 A. _$ E) |' o" ^ dir=gradt
    x0=golden(x,dir,hessin,b)
% B) B8 X" `4 B. ?* y    x1=x+x0*dir
- L4 _$ c2 h" ?% d/ ?5 l gradt1=matmul(hessin,x1)+b
& R9 C: _& S/ L' w% I1 O 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
6 w0 F# e( s$ k) D8 t* [  W  ?    iter=iter+1
  if(iter&gt;10*n)then
    print*,"out"
) T1 N& X$ I# W: s    goto 101
endif
1 ^6 p, d# F( z) u! N    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
4 u. x6 f8 H! N- s. j: i6 e0 j    goto 100
9 j5 s4 b# s& c- ]+ R    contains</P>
<>    !!!子程序，返回函数值   
1 _  N3 k$ F* G) \& r    function f(x,A,b) result(f_result)
' Q' L8 m9 u5 `( s# w    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
" o" G* B# N) t- [( ~    real::f_result
7 @1 O7 f! U% W& h    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
    end function f
( U# T! a8 G6 S6 E" ^8 g# Z !!!子程序，矩阵与向量相乘
4 U0 e: h, v* }. f subroutine vectorm(p,G)
  E4 _5 Y4 t' G: [& O! j; E real,dimension(,intent(in)::p
real,dimension(:,,intent(out)::G
n=size(p)
do i=1,n
    !do j=1,n
& N; {( L7 g4 r+ f       G(i,=p(i)*p
    !enddo
enddo
) w( E6 z$ a3 |$ W, V! N& Z3 K end subroutine
4 I7 n  t, T* K/ U; D
    !!!精确线搜索0.618法子程序 ，返回步长；
) r/ P# h  i7 P+ H, X    function golden(x,d,A,b) result(golden_n)
7 l# H9 F" U: m" _% q4 l. Z    real::golden_n
& q; c( t$ p# c8 ~  t) D    real::x0
    real,dimension(,intent(in)::x,d
    real,dimension(,intent(in)::b
% L  i( t* D# C    real,dimension(:,,intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
    parameter(r=0.618)
2 e) \2 X- k  S3 p9 }    tol=0.0001
- P/ t5 X4 g2 I! D& K6 v8 J. F8 L5 h    dx=0.1
    x0=1
    x1=x0+dx
    f0=f(x+x0*d,A,b)
    f1=f(x+x1*d,A,b)
2 O- B: _* V) P# A    if(f0&lt;f1)then
# @( [& X; f" s  h5 h, x4       dx=dx+dx
        x2=x0-dx
+ u9 M$ ]( D$ L6 ]$ F6 r" V        f2=f(x+x2*d,A,b)
+ E- ~  ~4 d& \        if(f2&lt;f0)then
           x1=x0
" k1 b5 b6 `6 j) c1 T        x0=x2
        f1=f0
        f0=f2
& K! L5 E. M' w% ^+ l        goto 4
        else
           a1=x2
8 J9 T' K- P5 w) ^" A* L' J        b1=x1
: J3 r8 G* A3 q  T+ @0 ?7 s, [        endif
    else
2       dx=dx+dx
        x2=x1+dx
& c& Y( J$ U0 X# m+ U        f2=f(x+x2*d,A,b)
, c9 |: Z9 f4 F' I  X1 P$ F& r9 R        if(f2&gt;=f1)then
9 Q$ E& S( o/ t7 C           b1=x2
* N! k0 ~3 y& W/ q2 r        a1=x0
$ t$ j1 j, ]) e  b. a        else
           x0=x1
        x1=x2
: n5 D( m; t, U5 w& z        f0=f1
3 P$ ^4 _- D, C6 [) k$ l: G        f1=f2
4 K4 V: F  S# V9 R0 G3 O, L        goto 2
        endif
* A; P6 f" n. z6 g5 A0 x# J" X    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)
( v6 X2 N/ N. `( A8 y6 b3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
    else
        if(f1&gt;f2)then
7 a0 \& N2 R4 g" y! T/ r        a1=x1
1 r: t/ b9 Y2 P4 |        x1=x2
  J2 X3 |) Y/ M- |7 V3 {        f1=f2
0 ]1 _( m) o4 @        x2=a1+r*(b1-a1)
! S  s6 K0 {% l6 t( b        f2=f(x+x2*d,A,b)
        goto 3
2 }1 p; A: j  j' w  g3 j2 `     else
        b1=x2
        x2=x1
0 [. g- @$ ?0 g- }. O        f2=f1
        x1=a1+(1-r)*(b1-a1)
4 F4 g3 n& M/ c8 N' }        f1=f(x+x1*d,A,b)
- i* ~/ p1 n% s2 G8 f/ ^: J, y/ \        goto 3
     endif
5 I. w) r% Z' W0 k    endif
    golden_n=x0
) C! K- X8 }+ l3 e& a    end  function golden</P>
<>
2 @7 a2 g2 D9 N+ y( y    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
/ e2 o- u# N; Z    subroutine gaussj(a,n,b)
+ ?; `" s; i2 J! ^2 A    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
  |! f0 m! o( d* s% L    do i=1,n
; }6 T; K$ @% I9 a. \0 D3 d% f       big=0.
4 B( ?% i4 p" E1 N- l* ?5 H       do j=1,n
       if(ipiv(j)/=1)then
          do k=1,n
          if(ipiv(k)==0)then
( D( L2 E( m/ o5 R* [, ?2 s/ M          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
6 u. ^2 T* h, i3 [) J+ n+ j           irow=j
           icol=k
% z& B- g- k& `: I! U       endif
       else if(ipiv(k)&gt;1)then
, a/ ^' m1 Q- j5 e5 P) x% b          pause'singular matrix in gaussj'
3 N- C8 w2 ~% O9 ]! P       endif
( ], O/ E* r3 Y4 \: ?$ V. `8 Z       enddo
    endif
    enddo
) h2 k9 p9 C% P) c* C; {    ipiv(icol)=ipiv(icol)+1
0 H: ~. b1 ]$ n; q/ A: O0 n0 M    if(irow/=icol)then
       do l=1,n
          dum=a(irow,l)
       a(irow,l)=a(icol,l)
5 y) R1 g  O3 _( d       a(icol,l)=dum
" D( U# ^! @/ P, b       enddo
: X1 m! v- f8 [% v7 P; E, v+ Z# _       dum=b(irow)
* x/ z8 {/ [* ]6 u* K$ Z       b(irow)=b(icol)
9 M/ C6 N/ W( V1 @9 w6 Z0 M       b(icol)=dum
    endif
- ?1 @- W1 [  {3 P; c    indxr(i)=irow
    indxc(i)=icol
& A; D" J! Q. ^+ g1 M3 t/ h    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
$ `, ?% h7 n3 z5 o$ g& s+ n    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
% }! ~: l" ]. X8 |, y    enddo
" M# f* P' C9 g# m; J; K8 p    b(icol)=b(icol)*pivinv
    do ll=1,n
. D5 G8 k8 J* B; ]) B9 p; z0 y$ z       if(ll/=icol)then
# Y- K9 d$ q8 g# e          dum=a(ll,icol)
       a(ll,icol)=0
2 p! K4 Z$ Y& ?4 Y  @- D       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
5 i4 x7 v% R8 Q% U, F6 [       b(ll)=b(ll)-b(icol)*dum
       endif
% C" _! X; R; B: A$ z$ h    enddo
/ y1 i5 r* ^+ G: ]$ @" V    enddo
    do l=n,1,-1
& d6 g) K$ f9 A2 D$ K% }9 Q       if(indxr(l)/=indxc(l))then
2 _3 a& R' M& N! h- O# e       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
2 ^4 G9 R( N5 _/ o4 Y    end subroutine gaussj
$ u, z' M$ y  {101 end
</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

楼主，怎么这么多笑脸，晕了，有没有matlab的