BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
- m9 Y! A7 W5 q( R: }; W& m5 K    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
9 P: v0 G, z1 l; b/ R    !!!dir实型变量，存放搜索方向；
' a) |5 j! i2 T( z; o  ~# }" q    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
+ S: K- x- M7 S5 k1 y: F( c    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
6 R9 j6 Z' o3 ]' |    real::x0,tol
  \$ ^9 \6 t' g: \  y! r9 w    integer::n ,iter,i,j
- W! d. ]  L* m3 E+ Y6 a1 _    print*,'请输入变量的维数'
    read*,n
; {: ?" h& T& Z" b  ~% b5 L    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
! v& [* G$ `1 p9 [# h    print*,'请输入hessin矩阵'
# [' a) X- \" n. p    read*,hessin
$ e3 D" R' O$ n/ ^    print*,'请输入矩阵b'
$ c, c! v8 I4 O# @% [    read*,b
    iter=0
tol=0.00001</P>
<> do i=1,n
1 `) x$ D. ~7 n2 _0 d+ i) u    do j=1,n
+ c2 K. g: @3 T$ o) q" L       if (i==j)then
       B1(i,j)=1
    else
       B1(i,j)=0
9 l$ g3 n, C( Q1 k7 n8 k* Z2 @7 |    endif
    enddo
' P! c6 x9 n* }/ t5 ~) J enddo   
* m' f9 b. U; g" L9 {! j! }7 _    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
9 z( k  Q: {9 W7 D3 d2 P        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
     goto 101
    endif
0 k$ F" x4 [0 X( w1 N4 i  k call gaussj(B1,n,(-1)*gradt)
* C' Z5 _, O. V- r dir=gradt
    x0=golden(x,dir,hessin,b)
7 L  G4 L: ?' v2 H9 ^8 w. f    x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
s=x1-x
; K: e- L$ O0 R. n; c2 _ y=gradt1-gradt
call vectorm(gradt,G)
( ^1 `5 h' k# G" A G1=G
2 P9 X) }+ o* u( f call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
- m: W0 u, {; X$ w! `% a x=x1
. T$ i* ]0 S) V3 P) M  z- j7 T gradt=gradt1
    iter=iter+1
  if(iter&gt;10*n)then
  {8 G, v5 ]9 [# J6 H3 [) C    print*,"out"
    goto 101
endif
    print*,"第",iter,"次运行结果为",x
2 t" f7 Q% v! K8 c& h% N7 N print*,"方向为",dir  
    goto 100
    contains</P>
% Y. C" O% W( t<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
+ \8 v; l( Y7 C3 A- z" x: F    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
$ D4 m1 a. B% f8 c8 ]% w! o    real::f_result
, Z; r$ i& a+ {* }3 _    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
' ?) `8 [" c/ w$ _    end function f
+ k4 t, i0 Q3 M6 X6 [$ q !!!子程序，矩阵与向量相乘
1 _2 i% g6 b+ s1 x subroutine vectorm(p,G)
real,dimension(,intent(in)::p
real,dimension(:,,intent(out)::G
n=size(p)
9 w5 n5 p% l$ k4 C do i=1,n
4 F" Y" z, d8 {& p    !do j=1,n
       G(i,=p(i)*p
7 m& d: Q$ D% S  v2 y: y# `) g    !enddo
) [( W  N4 Z- A+ H! s enddo
end subroutine

    !!!精确线搜索0.618法子程序 ，返回步长；
+ U! e3 V3 u/ Q4 {- E" C2 z    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
1 }0 Y- J/ V0 n7 ~1 ~    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
9 W6 [, }4 p( R+ l. }    f0=f(x+x0*d,A,b)
9 Q! C3 O8 d0 K# [. R7 S    f1=f(x+x1*d,A,b)
5 [# `0 w; \# G/ W, B    if(f0&lt;f1)then
4       dx=dx+dx
        x2=x0-dx
        f2=f(x+x2*d,A,b)
5 z% T8 |9 k8 }/ w, o0 H        if(f2&lt;f0)then
           x1=x0
        x0=x2
/ l5 V+ V7 @# S- m        f1=f0
        f0=f2
& n% H/ K7 i: K: U        goto 4
/ ?/ o. [4 C: U# k. f' M$ i5 Y, a        else
; [; w& p# @1 j           a1=x2
! z7 i4 U; e( K  u        b1=x1
1 h: [0 Z+ K( W& F+ R6 y        endif
    else
2       dx=dx+dx
; r1 S5 S& `2 f- k& {' m% @        x2=x1+dx
        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
, ~* z4 X$ B6 A6 _# k* J$ ~           b1=x2
        a1=x0
2 V7 _  N! ?3 w7 s- Z$ a        else
: g! v% z# ~2 A$ Q- V- `9 X/ K8 X4 U, G           x0=x1
        x1=x2
- h( D, I8 h' {7 z; |2 Y        f0=f1
$ U, V( a! P: ^0 N4 }9 c- ]: ?        f1=f2
; X) g; ^( S- R6 e4 [        goto 2
        endif
    endif
) V3 w+ j0 J* [5 u: D3 f/ e+ `    x1=a1+(1-r)*(b1-a1)
( R& n3 m* y- }! ^& w/ o& S    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
; `: T7 \' L$ X        x0=(a1+b1)/2
    else
        if(f1&gt;f2)then
        a1=x1
7 V5 v6 u* K4 w$ r$ T4 v5 g        x1=x2
        f1=f2
        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
        goto 3
+ ?8 z6 d$ T/ V* y     else
        b1=x2
7 u3 g# g7 u& B7 C* |! J5 M! f! V        x2=x1
: z* d- N: d, V" E, |        f2=f1
8 s; A, g, E" j( c6 k3 o        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
; y/ S" _7 f3 |" Q" ^) i& U        goto 3
! R# m) I0 `3 v# N" n. r" q     endif
    endif
3 d5 f6 w' U$ Q    golden_n=x0
# S/ S! \4 f( I- d  r+ S* t    end  function golden</P>
6 q0 \. {" H- k9 X<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
- j: j. k% p) [# O+ q: x    subroutine gaussj(a,n,b)
    integer n,nmax
% u, |0 w$ x( t& |: Q% Y: p" ?; K    real a(n,n),b(n)
. N- s. r0 ?  w) k0 \! [    parameter(nmax=50)
( Y! \' Y  @* o& d2 C6 n    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
4 }  G. j4 w* K0 F1 y+ t) g    real big,dum,pivinv  
; u1 O% O, M' C: M5 C    do j=1,n
       ipiv(j)=0
8 \+ @* I4 f8 ]2 \1 Z    enddo
    do i=1,n
7 R- t2 V( ~! U* k       big=0.
       do j=1,n
7 T% k2 I* m, b  }3 L/ X       if(ipiv(j)/=1)then
- }& L/ Q( \/ ~8 b          do k=1,n
/ U$ g$ @$ K  k" t          if(ipiv(k)==0)then
) i9 Q9 M% B5 E. r  Z; d          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
           irow=j
+ ]4 S* ?  b' Z3 o           icol=k
       endif
       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
       endif
: V$ z! z/ s% k, L9 z. v$ u       enddo
, H" V, U: s' x, V1 _    endif
    enddo
6 E6 A6 Y, A0 N/ k    ipiv(icol)=ipiv(icol)+1
$ w% F! V! i$ w' x% U3 `    if(irow/=icol)then
       do l=1,n
          dum=a(irow,l)
       a(irow,l)=a(icol,l)
8 ]- v/ {3 v4 [0 X3 K9 c       a(icol,l)=dum
# U% _% @7 y5 ?, w" T       enddo
       dum=b(irow)
, w0 N% k' v" C9 a  J5 T2 a4 E       b(irow)=b(icol)
       b(icol)=dum
* x9 M- C' G5 P; I% f7 z4 J' K    endif
" Y: f( \$ v% e& o. c1 q. Q6 G    indxr(i)=irow
    indxc(i)=icol
    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
* Y0 A6 `( q% {! ~) l    pivinv=1./a(icol,icol)
2 B  A" a- h; j/ M. s# \    a(icol,icol)=1.
6 K5 H! t+ {! O: f; }9 [    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
    enddo
3 t5 f  W+ ~  B6 b# M9 y    b(icol)=b(icol)*pivinv
    do ll=1,n
3 Z! U1 x, y' L" A3 D" E& O6 b       if(ll/=icol)then
: V0 ?, K3 Y1 B( Z          dum=a(ll,icol)
( w, R+ I* f! D& }, K' g       a(ll,icol)=0
       do l=1,n
9 v# j7 Z! |9 H          a(ll,l)=a(ll,l)-a(icol,l)*dum
: q: b. r/ Q  E( n: y1 ?       enddo
  O% p/ w, T* J/ G       b(ll)=b(ll)-b(icol)*dum
       endif
& ?4 ]4 m. F" [0 e8 B    enddo
    enddo
    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
  Z* n  N' ~% q3 O; R2 L( ~       do k=1,n
; P' q5 D- T) s          dum=a(k,indxr(l))
       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
) C, d- l' O; \       enddo
* U" ^& `# J2 C+ Z! g# J8 c, p    endif
    enddo
    end subroutine gaussj
101 end
  \2 S5 k! E8 z) `# b9 T5 y</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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