BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
  t8 ?: b, q0 E; }    !!!输入函数信息，输出函数的稳定点及迭代次数；
2 z7 R/ g  u& `    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
0 F$ |" Y; M: J# N( O( N    !!!dir实型变量，存放搜索方向；
+ a2 O+ Y6 L! C: T% I    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
4 N( D6 m4 G0 S+ `    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
    integer::n ,iter,i,j
0 P, O; Z4 X# v  Z    print*,'请输入变量的维数'
7 _9 _4 ?- [1 c" }# g    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))
8 P) ~  ]2 P, U+ x) U7 s    print*,'请输入初始向量x'
$ j6 T' h2 j* q2 t" N6 y( N    read*,x
, y" h" T. q, n  e: O    print*,'请输入hessin矩阵'
7 ^. p7 I( o5 A0 q5 N6 E    read*,hessin
    print*,'请输入矩阵b'
    read*,b
    iter=0
" ]& R9 _/ c$ R6 @! h! ? tol=0.00001</P>
<> do i=1,n
8 C$ j, g; {' M% H    do j=1,n
       if (i==j)then
       B1(i,j)=1
9 D2 c" W6 n& x8 b    else
       B1(i,j)=0
% |& a; g" u+ @1 o; ?2 R    endif
    enddo
enddo   
    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
: {  a7 T4 n" z4 N  a        !print*,'极小值点为：',x
6 m" K- r' h. S$ _6 Q7 T     !print*,'迭代次数:',iter
     goto 101
    endif
! ^8 T! |# @; o. m) C call gaussj(B1,n,(-1)*gradt)
dir=gradt
    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir
. t) S! {( n( g0 S( X9 F; Q gradt1=matmul(hessin,x1)+b
  E  [- E  S- C, b8 L/ n s=x1-x
- T# P! w  s8 |0 n# ?' c- v: i8 J4 o y=gradt1-gradt
call vectorm(gradt,G)
  B* [$ m! D! \* G G1=G
9 f4 q( |- m/ D& k9 {8 K  } call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
6 K3 K0 M: d3 u  { x=x1
gradt=gradt1
, S* A7 p0 d. d9 _9 S. B6 }    iter=iter+1
( ]5 z8 O! [3 [  if(iter&gt;10*n)then
  o8 n; n4 o8 h/ n% p    print*,"out"
    goto 101
endif
0 I$ a3 e+ @) H1 w    print*,"第",iter,"次运行结果为",x
0 L( @+ Z# u- g% ^: J print*,"方向为",dir  
! m. H1 s5 I! q5 r* D7 P    goto 100
    contains</P>
" |/ V6 J: G% Q  ~* f3 c8 o<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
    real,dimension(,intent(in)::x,b
. Q* z- {# D0 f& d1 O) _    real,dimension(:,,intent(in)::A
8 \  T) ^& j5 ~& O2 ?/ x' t; `    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
1 s) q% a/ p( h& V    end function f
!!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
real,dimension(,intent(in)::p
$ q" _3 s" G1 w4 K. k4 Y' \+ b) P' f real,dimension(:,,intent(out)::G
1 a+ A/ F2 m! I: u$ W n=size(p)
do i=1,n
- [( V% I* e) t$ g    !do j=1,n
( @7 V( s) F' o4 |! f1 A! M       G(i,=p(i)*p
    !enddo
enddo
end subroutine

+ p; S7 {. s2 O    !!!精确线搜索0.618法子程序 ，返回步长；
    function golden(x,d,A,b) result(golden_n)
9 b6 m5 H) p* Y" |    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
- \  H8 V5 j! C" [! P  R    parameter(r=0.618)
0 p- X+ O) \4 [$ p5 v) _/ I4 Z9 o# [    tol=0.0001
/ ^$ i7 x4 s; f: r" X: D    dx=0.1
    x0=1
5 k# c% I  ~; {& _/ T) L2 B    x1=x0+dx
    f0=f(x+x0*d,A,b)
- k; b) l+ J4 |$ A( P5 J$ F& e  x    f1=f(x+x1*d,A,b)
. n% a% o# x6 ~; ?( _2 n: J    if(f0&lt;f1)then
$ n: \  o8 e' G: F5 k4       dx=dx+dx
        x2=x0-dx
4 c1 A( m  z; \6 I        f2=f(x+x2*d,A,b)
        if(f2&lt;f0)then
# ?$ z2 g/ }- p           x1=x0
  e2 M# {, ^4 c+ _) c        x0=x2
        f1=f0
        f0=f2
        goto 4
. P1 r8 X2 d* `8 V        else
+ S' U2 r6 Y- P9 _           a1=x2
) z  r: n% [/ S  _        b1=x1
& C. x( f: v: `  g3 q7 f        endif
5 u( K* L9 j# u7 l( G3 k  B; ?! s    else
2       dx=dx+dx
        x2=x1+dx
        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
! Y0 {% U8 q- [' H           b1=x2
# R) |) T! l) R1 Q( Z; ]4 s- ~        a1=x0
: R9 P, n3 S. ~# T        else
           x0=x1
' [# n. P/ ]8 t* h0 @8 p) m* Z        x1=x2
. F% z0 M) {) N9 S5 Y  S        f0=f1
        f1=f2
- n: i/ \8 c- e/ S- `        goto 2
' Z* T2 G( L8 k" x/ O        endif
    endif
    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
: x1 b4 T8 d. @2 b  `5 Y4 d, }, |3 N    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
* J9 o8 L' o$ x; R    else
- b6 `# [( S- J        if(f1&gt;f2)then
        a1=x1
' U; {" r$ d' [7 K# y0 W        x1=x2
        f1=f2
$ W, C% o. e* T5 n( w7 w        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
        goto 3
( x3 g& U% _1 U$ U' I- U6 b4 {) }     else
- J6 X% u5 B: W8 h: D8 Z        b1=x2
        x2=x1
        f2=f1
        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
        goto 3
5 a/ x( K8 f* m  B: q     endif
    endif
    golden_n=x0
6 t8 }6 Z$ k. S, t3 B/ g' {    end  function golden</P>
8 z' L3 B3 v  [& i<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
+ p% T1 I1 \* Q  r    subroutine gaussj(a,n,b)
    integer n,nmax
    real a(n,n),b(n)
6 R  m- U% u* \$ `    parameter(nmax=50)
" E( Z) ^* i+ |+ \* p    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
$ o5 x5 n* b3 U0 h/ T    real big,dum,pivinv  
    do j=1,n
+ q6 _% ^. D7 c/ ]: D) v/ h       ipiv(j)=0
    enddo
/ u+ y3 `& r. p; k0 X    do i=1,n
       big=0.
7 r5 v& i+ e  O9 U3 v) b2 W8 z1 |       do j=1,n
3 ^6 Y0 d2 L1 }       if(ipiv(j)/=1)then
          do k=1,n
          if(ipiv(k)==0)then
, h0 a0 D- `: E$ @          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
3 J% }3 R) L5 R# t0 w1 J! q; h           irow=j
( y% z$ f$ a8 ]5 O- }           icol=k
% c0 {8 Q3 I' y. l. D, ~( S       endif
       else if(ipiv(k)&gt;1)then
" D" o& U+ h1 G& a          pause'singular matrix in gaussj'
       endif
       enddo
    endif
    enddo
6 g$ M5 G( M$ y0 b$ O    ipiv(icol)=ipiv(icol)+1
- n0 p, \+ J9 ]3 B2 H! r: F    if(irow/=icol)then
       do l=1,n
          dum=a(irow,l)
) o- u$ p9 Z4 Y, d! ?; W       a(irow,l)=a(icol,l)
  D: V- Y9 o% J       a(icol,l)=dum
       enddo
4 c7 Z' {+ R7 q8 s: I6 n       dum=b(irow)
       b(irow)=b(icol)
0 ?% s7 ^& o" l* _# R1 p7 g1 A       b(icol)=dum
    endif
2 D# q+ w( [3 n! E    indxr(i)=irow
    indxc(i)=icol
    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
, [8 ~" c! v' U& ^    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
) d! B% b3 E1 O) C8 N    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
3 o8 z9 o# {1 M6 p( w. Q    enddo
    b(icol)=b(icol)*pivinv
4 V. w! s" \. P6 M1 h0 i    do ll=1,n
( w; z7 w/ B/ m  g/ \       if(ll/=icol)then
6 W3 K; z8 ^/ P          dum=a(ll,icol)
       a(ll,icol)=0
4 H: G, @5 a( t/ t       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
) O" b0 P6 |5 z3 X       enddo
       b(ll)=b(ll)-b(icol)*dum
+ H  h4 h1 N/ _       endif
    enddo
    enddo
    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
0 ?& K% f" P; j8 _       do k=1,n
          dum=a(k,indxr(l))
( a# G  H4 G6 E9 U5 Y; c# k; V) b       a(k,indxr(l))=a(k,indxc(l))
! e9 y. R8 l" R5 r       a(k,indxc(l))=dum
       enddo
) F$ Y1 e* |8 \- h4 r4 }' j3 O" q1 l    endif
    enddo
% s8 W/ z8 Y* o5 O8 p  D1 |    end subroutine gaussj
/ u4 ?& L* A4 T" @101 end
3 Y! _0 B2 t. P; d9 g; T/ o5 U/ P</P>
1 x" D6 |1 F" k* I) v9 P: M. u<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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