BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
5 o' a1 b5 v6 H6 }- K0 C2 S    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
" H. }+ ?) i, Z7 M* }2 S% r2 j    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
1 w$ J; d( Z$ e6 }    !!!dir实型变量，存放搜索方向；
- G5 X8 Z4 G: i4 j6 A    program main
' f" Y' p0 o% d: a8 ^8 b8 l, w2 F    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
( S7 b* a0 h0 q  i- {% S& v    real::x0,tol
    integer::n ,iter,i,j
/ v( F5 g* ]: G, _* w    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))
% V  @9 x1 Z$ Z: V5 v4 m    print*,'请输入初始向量x'
  A  n3 z, y* J8 A" j, Y    read*,x
5 l' }( P. E) I    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
/ j( ~0 ?" o$ u6 h! D( k  a( U; R8 k    read*,b
    iter=0
tol=0.00001</P>
% [! z' F" l+ G6 R2 Q( c<> do i=1,n
8 V; j  `' V' p    do j=1,n
       if (i==j)then
       B1(i,j)=1
    else
       B1(i,j)=0
6 w% g6 B3 n' O$ ]: F    endif
$ c6 L' u! r, t2 F9 a! _4 O  O7 {    enddo
/ K% [  J$ s5 C  J3 P) i enddo   
: @; M/ U3 g( r* G8 V8 `4 o% @; ]7 ^    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
5 o5 r0 L) j9 w8 p+ \, a2 ~8 s5 r        !print*,'极小值点为：',x
5 O& X( X8 n- C( q: C/ x9 v/ j3 e     !print*,'迭代次数:',iter
     goto 101
    endif
call gaussj(B1,n,(-1)*gradt)
9 ~, I: u) b* c1 d! z. W8 H' I dir=gradt
    x0=golden(x,dir,hessin,b)
* A8 M! e9 W8 y# T; M" j( I, n    x1=x+x0*dir
+ J8 R; t* d& W# E- } gradt1=matmul(hessin,x1)+b
s=x1-x
; N. ]% K3 ]9 V  `* W5 S 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
# t' v; n; P) n- m. d8 { x=x1
gradt=gradt1
    iter=iter+1
, H6 L9 B5 D6 T. R7 k& |  if(iter&gt;10*n)then
& f8 l5 t+ z' K% W& T: \- X- `    print*,"out"
    goto 101
8 n3 i5 [4 R+ f' h endif
    print*,"第",iter,"次运行结果为",x
  c8 r# d7 y1 x1 S5 p# P print*,"方向为",dir  
, s+ h8 v$ P0 W9 |/ G2 g0 `    goto 100
* q* t7 I+ H% {! I6 {& N9 B& S    contains</P>
% t, D. o& Q& \! V' @  T<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
1 q8 B* G: U- l! Z! l; V    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)
& b7 o0 _+ u- j: C: R9 c    end function f
+ i6 g5 L0 J; \9 w$ d( ?/ v' B !!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
real,dimension(,intent(in)::p
real,dimension(:,,intent(out)::G
  H, W5 c& E% l$ w1 H( ] n=size(p)
do i=1,n
    !do j=1,n
       G(i,=p(i)*p
    !enddo
enddo
end subroutine

    !!!精确线搜索0.618法子程序 ，返回步长；
0 t. C; z& ^7 Z+ f1 K    function golden(x,d,A,b) result(golden_n)
5 k9 }& o  x# ^2 D    real::golden_n
    real::x0
) `# a! t4 I: Q+ J* y# P    real,dimension(,intent(in)::x,d
    real,dimension(,intent(in)::b
    real,dimension(:,,intent(in)::A
* {6 I+ v% e; O, P6 Y    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
    parameter(r=0.618)
    tol=0.0001
    dx=0.1
2 a) [* E# o1 L  ~# U& @    x0=1
  h1 {. |0 I/ f    x1=x0+dx
    f0=f(x+x0*d,A,b)
/ o$ P2 p, Y( }" |    f1=f(x+x1*d,A,b)
1 R5 N' D5 `3 J$ K- |. p    if(f0&lt;f1)then
4       dx=dx+dx
        x2=x0-dx
1 B7 [+ o& z) q- w% T8 |$ _        f2=f(x+x2*d,A,b)
- a- b: I1 f  C% K        if(f2&lt;f0)then
           x1=x0
' y' o$ P% f$ j9 E! M  `        x0=x2
        f1=f0
8 J  o" E# Y6 S8 n        f0=f2
        goto 4
        else
* P, G" A7 N+ e1 b           a1=x2
( X" ~% \; j! v2 J6 c        b1=x1
        endif
    else
9 F2 P2 E4 {2 G# t0 h2       dx=dx+dx
; P: n# e5 i% I" v- D        x2=x1+dx
        f2=f(x+x2*d,A,b)
; B; M& \4 F% }/ b! S  N" e' g* o        if(f2&gt;=f1)then
           b1=x2
        a1=x0
        else
           x0=x1
        x1=x2
        f0=f1
        f1=f2
        goto 2
( D) L/ n$ D' t9 D) v' f9 z        endif
    endif
    x1=a1+(1-r)*(b1-a1)
& }% q: _+ A7 y5 Z0 W    x2=a1+r*(b1-a1)
0 e" o$ f! D, o    f1=f(x+x1*d,A,b)
; V+ _6 _5 U2 x5 g  }& O    f2=f(x+x2*d,A,b)
3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
7 q- `/ D- p1 X' V    else
  G: Z; P! t: e  X& N  F. j        if(f1&gt;f2)then
  [5 X: J* z. S+ J/ ]- z+ V0 s        a1=x1
        x1=x2
        f1=f2
        x2=a1+r*(b1-a1)
" }9 Q4 o2 X, u! G4 e+ I1 C        f2=f(x+x2*d,A,b)
" n  m% ]3 i7 ~. i5 j8 Y: T        goto 3
     else
, b! ]8 A5 t  X+ M/ _1 W        b1=x2
        x2=x1
4 z* D; S9 J7 C8 @0 `' {0 R        f2=f1
        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
        goto 3
     endif
( y  T: }* r2 A; A    endif
    golden_n=x0
    end  function golden</P>
<>
% n# e8 K5 u: |4 {1 `, Y0 @$ Y    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
    integer n,nmax
) D2 t. j  u( Y' Y1 l- W0 h' x    real a(n,n),b(n)
    parameter(nmax=50)
    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
! t/ E: y! F. z3 z( C5 U; z. I    real big,dum,pivinv  
  [# A* ^* I: b( t3 @    do j=1,n
       ipiv(j)=0
, y- A! M* O/ g+ G- o2 C    enddo
& S4 b5 U2 G6 k0 s' g' g    do i=1,n
       big=0.
& y4 @; @4 ?8 _( n# C+ h       do j=1,n
       if(ipiv(j)/=1)then
          do k=1,n
          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
1 J' ]& [$ [7 ~% r* {( s! f* n           big=abs(a(j,k))
+ j. T& h5 @1 V' ]3 P# \/ ], M0 o# Q' Y           irow=j
           icol=k
3 b# x0 \# l" l1 U       endif
       else if(ipiv(k)&gt;1)then
3 u6 m! l  P$ ^1 ]3 h0 x* h          pause'singular matrix in gaussj'
       endif
       enddo
) N  ?. G1 I/ H- F6 q$ [' X    endif
! T+ \# \2 S! X+ j" O$ q/ [9 f    enddo
0 N" i- x' a- D/ i: m9 T& z7 B    ipiv(icol)=ipiv(icol)+1
& n4 y& ~$ u6 m    if(irow/=icol)then
/ }9 _* d% P* A4 D: x. r$ [       do l=1,n
          dum=a(irow,l)
       a(irow,l)=a(icol,l)
. [) r: ]5 c& c4 F. J       a(icol,l)=dum
; W7 S- g; J) B( t- W6 T* x" j       enddo
$ `3 a# j7 Q0 _) H1 b3 [       dum=b(irow)
4 `" _  t$ }0 N7 D* P( _       b(irow)=b(icol)
       b(icol)=dum
- p& C2 L- @+ A, y    endif
    indxr(i)=irow
" ]: s  v* u' E9 Z3 v    indxc(i)=icol
1 F1 g2 S' ]0 k2 C    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
. x8 X1 L0 W, M$ f" v; d    a(icol,icol)=1.
    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
3 y) v& ^# y* f9 Q  b2 s3 f0 m/ i    enddo
    b(icol)=b(icol)*pivinv
7 m( [* {, S- [/ T: M5 [    do ll=1,n
- i/ ~8 [- ~" V. p% _$ a- C# Y       if(ll/=icol)then
/ S2 N* O3 \4 O$ Z: Z0 U5 R- b          dum=a(ll,icol)
3 l1 A) a" Z: t: z5 N& z       a(ll,icol)=0
8 u9 v' I( m. Z1 g       do l=1,n
& i- [6 g9 J% J# c. Y1 O$ Q! |8 L          a(ll,l)=a(ll,l)-a(icol,l)*dum
, L( d& ~5 W0 `+ j       enddo
7 I' F/ @6 \- q: y" }4 \; [       b(ll)=b(ll)-b(icol)*dum
       endif
+ B/ X" U! H+ x, r    enddo
4 K7 O# r  ]" ~    enddo
    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
       do k=1,n
          dum=a(k,indxr(l))
       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
       enddo
- R9 f) R) H$ |5 I- a  D    endif
4 Q$ e+ v" K) p+ A! ]8 _    enddo
    end subroutine gaussj
101 end
</P>
! s/ `5 \  b2 v3 l! r; ^9 T<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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