BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
# S4 T# |& y' T/ T2 {) Z8 @    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
: F1 q& |( ?. a1 C% F# U0 ^    !!!dir实型变量，存放搜索方向；
% i! P7 F2 A# y0 x* e( e    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
9 ]) N! b! c' |8 H4 ]    real::x0,tol
7 F$ W8 v  c# [    integer::n ,iter,i,j
" e$ ?7 R& M, W    print*,'请输入变量的维数'
! @0 v% ?" I& t6 j+ S3 e% B    read*,n
) U4 o2 f* W: R. E9 E7 u    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))
5 O9 ]2 K- h8 C0 o    print*,'请输入初始向量x'
    read*,x
- P* r' N" `) T( G    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
) u3 ]6 k/ p; _, C5 s% W    read*,b
    iter=0
3 u5 F# u& G( h& }1 Y, j tol=0.00001</P>
<> do i=1,n
    do j=1,n
% I. Q- J4 ~; @: I5 x4 {/ |       if (i==j)then
& Q" v6 X6 y" d3 o# I! |$ l       B1(i,j)=1
0 u' S5 y* M: a. d    else
       B1(i,j)=0
3 K5 t* p. ]$ [, W  ^, f    endif
+ T  A3 T* _1 _( \    enddo
7 S  V. A6 w) o" K9 _7 ?- C enddo   
) \' D* ^) d0 f/ C( z    gradt=matmul(hessin,x)+b
# q& f2 Y) e2 i! d* ?100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
9 G6 F) x+ h4 r( c        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
     goto 101
6 n: z. a* n# U9 [6 u    endif
call gaussj(B1,n,(-1)*gradt)
! |) j9 z/ p  _/ U3 F4 @ dir=gradt
    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir
; y2 }& Q* i: f  w6 [- G! u gradt1=matmul(hessin,x1)+b
s=x1-x
8 @5 _( L+ N  x( R& I2 a y=gradt1-gradt
3 H. q5 l: Z7 X$ M/ N6 h call vectorm(gradt,G)
G1=G
call vectorm(y,G)
) k- o* ?4 X5 c. C B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
5 `& q- J% D" X5 G3 m3 O. L gradt=gradt1
    iter=iter+1
  if(iter&gt;10*n)then
( \; @5 j* B# B9 N# e    print*,"out"
* z$ G, ]4 N2 O+ j% n; V* e    goto 101
endif
6 J. n2 O+ x; q    print*,"第",iter,"次运行结果为",x
8 X1 S7 K3 u' [) d print*,"方向为",dir  
$ q+ M/ Q( p( n9 W! B    goto 100
2 [% g: k. c* a7 F" x5 |; r6 o" v    contains</P>
<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
, M' B" s% Q" D- z- z( Z# M    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
    end function f
' F1 N! y! ?  B$ A! e6 e' U2 ^: u* p( H !!!子程序，矩阵与向量相乘
5 t: j3 V' i4 W1 V subroutine vectorm(p,G)
" M. O& \+ H0 g7 ? real,dimension(,intent(in)::p
; _3 u$ Y* _$ s, l$ `$ l real,dimension(:,,intent(out)::G
n=size(p)
do i=1,n
- w3 Y. n6 l  [- C+ k! ~) q    !do j=1,n
$ v5 l; S) I. Q. y) }2 p       G(i,=p(i)*p
& T# X9 h* r, T1 l+ ]    !enddo
' E+ t* l1 L0 H% s3 @! V enddo
' c  o) X0 X1 V9 \. [. V end subroutine

" s. Q1 D% J0 T  B4 w3 b$ t" `- x    !!!精确线搜索0.618法子程序 ，返回步长；
    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
/ p5 N6 W) i% Y9 ?. j5 n    real,dimension(:,,intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
% ~; Z! ^+ w1 U( d    parameter(r=0.618)
    tol=0.0001
    dx=0.1
; @6 u# k4 N. T$ Z. ]4 N* w    x0=1
4 ^/ p9 C1 u+ r1 t. ^8 n8 p    x1=x0+dx
    f0=f(x+x0*d,A,b)
8 J1 o: @1 i! S  M  @4 W9 n    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
; q: l4 w( `/ u7 n) s        f0=f2
        goto 4
        else
, ^- ~$ y& W( W+ e           a1=x2
1 i$ u8 W* ]) q# h        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
( z4 e8 c1 m) P0 a        else
, x: a% s# I$ l" F$ H7 M1 h) X           x0=x1
7 }5 |: ]- E4 c; J8 a6 `        x1=x2
        f0=f1
5 o5 B5 p5 u! Z  w# H        f1=f2
        goto 2
        endif
' F+ Z! V$ C/ C0 `3 u; m( v3 U# @& P5 v    endif
    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
4 ]! h; q7 b+ b/ I1 E    f1=f(x+x1*d,A,b)
    f2=f(x+x2*d,A,b)
3   if(abs(b1-a1)&lt;=tol)then
; k) d, @* e& t, n3 [& t, @        x0=(a1+b1)/2
1 @  n' O' C5 g. }1 W2 E8 f! ^    else
        if(f1&gt;f2)then
. t3 j, W( ~5 p: J4 R; z. ]        a1=x1
6 X/ |9 [& I9 e# P. p7 ]8 L6 E        x1=x2
        f1=f2
5 }2 U6 ^! x: C: r0 g        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
3 M+ C3 _, t2 N* j9 t        goto 3
8 m/ `" f1 l' U     else
        b1=x2
  x# q  r& N1 I# |( D1 Y# w( O        x2=x1
        f2=f1
. F' Z/ ^; m/ |& o        x1=a1+(1-r)*(b1-a1)
6 Y4 @% T4 c. P$ M: R        f1=f(x+x1*d,A,b)
        goto 3
     endif
, L% o5 H% k' x! G' s. e( x( a9 E+ f    endif
    golden_n=x0
! G. c1 m+ p3 M4 y8 y1 g" ~/ M    end  function golden</P>
<>
) _6 D) O) Y+ N    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
; K* x5 K, T6 @' R5 ]) n    subroutine gaussj(a,n,b)
& @0 |; U- x8 u& t( H" q    integer n,nmax
    real a(n,n),b(n)
" V- R4 Q. I% \7 q  o0 ~+ O$ ]+ F    parameter(nmax=50)
1 y* N# ^* L8 ~0 _5 E: @: ?    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
3 k, C2 v. ]5 |) M4 h4 d4 n  T    real big,dum,pivinv  
6 @+ T, P! u% p2 H" [$ B0 g9 |$ ^    do j=1,n
       ipiv(j)=0
( K  A, V, S( F' N9 d' h    enddo
    do i=1,n
       big=0.
. x: T* B' F* A; u7 J+ n) I! q' M$ y       do j=1,n
2 C  i1 {+ R) w3 I5 F       if(ipiv(j)/=1)then
1 \; E; O) @) I$ f4 }7 K& T& H. ~+ g7 `          do k=1,n
& n4 z3 ]5 Y" g6 ]1 k$ z; s+ d          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
3 y' a' J" E7 M3 b. V2 T           big=abs(a(j,k))
           irow=j
' k, s" n& U( l% m. @           icol=k
2 x# ?0 g  S! M4 S# a9 x2 `# x       endif
5 E0 S0 _: i# W* a       else if(ipiv(k)&gt;1)then
, L  ?/ N5 u2 S" U, r8 T          pause'singular matrix in gaussj'
       endif
: |' v  [) w5 b9 @% O( R       enddo
    endif
    enddo
( u: I& ^4 R/ L1 H6 u3 b( `8 H    ipiv(icol)=ipiv(icol)+1
    if(irow/=icol)then
/ g' c* A* `: i1 ~; ]- z1 Y       do l=1,n
) L, ]) t0 b1 v6 u          dum=a(irow,l)
" {+ \/ v( P! R2 S1 O" t( |       a(irow,l)=a(icol,l)
       a(icol,l)=dum
       enddo
       dum=b(irow)
3 G) x5 n8 \7 I5 A. I* I       b(irow)=b(icol)
       b(icol)=dum
    endif
    indxr(i)=irow
# e/ I) i2 v3 o3 h; \    indxc(i)=icol
( H1 Q- M" w# [! A8 g    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
1 Y! l$ Q/ ~' U' j    a(icol,icol)=1.
    do l=1,n
! t8 i1 u/ O' a2 E  b        a(icol,l)=a(icol,l)*pivinv
    enddo
0 G; ~5 V: U/ }0 D    b(icol)=b(icol)*pivinv
& ?) A' I$ B0 a7 G" [    do ll=1,n
       if(ll/=icol)then
, Y3 K2 X+ L8 Z' X& Z6 a/ r: @# l          dum=a(ll,icol)
       a(ll,icol)=0
+ b% S- F$ w& K# `# d       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
& w# e( [& s5 r, h6 Z' h       enddo
       b(ll)=b(ll)-b(icol)*dum
       endif
2 b  t6 s) F2 }0 A    enddo
. i2 l7 M$ m. w    enddo
* w, W6 }- t$ x- B2 G    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
) v# l. s5 y; C9 l# z  B- s) s       do k=1,n
: p, `9 [+ V/ H5 z) b          dum=a(k,indxr(l))
       a(k,indxr(l))=a(k,indxc(l))
; A' }* A$ r9 G& D       a(k,indxc(l))=dum
       enddo
& g7 i) e6 ^) H) D3 ]% g    endif
% U: C/ ?6 b2 T; N    enddo
2 l& |& e) _% ?6 S    end subroutine gaussj
0 E2 G2 ^/ M" \8 B5 R" o) {101 end
</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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