BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
  L0 ]* o- ~! T; A$ l" z( V    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
3 B- |! n1 N. f3 m) a    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
/ i+ a9 x/ ~+ G  \( D( j; g+ S: J    program main
$ E/ k. C2 q" }% T$ a+ R    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
: l) |6 j* b4 C    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
    integer::n ,iter,i,j
/ U1 D1 m' i- G3 R& E+ u- B8 F, {    print*,'请输入变量的维数'
    read*,n
, U( S' p0 n5 M2 w- a% Z2 O    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矩阵'
    read*,hessin
0 j' p* }# K" ^1 a0 n; G  J; `  d    print*,'请输入矩阵b'
' x6 B: I  D$ n- c    read*,b
6 L% |6 V3 `% t9 l  G2 h/ E    iter=0
/ A& M0 i1 v4 M0 ~, H tol=0.00001</P>
<> do i=1,n
    do j=1,n
       if (i==j)then
       B1(i,j)=1
    else
+ q  u& x. C: b4 n' k       B1(i,j)=0
2 V. C: [: m( E) c( ]$ D    endif
    enddo
7 \! Y& G) D7 G) j1 K enddo   
    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
6 M2 p* q3 ?4 o9 e& l        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
: A. E) C* y8 m2 P  p+ ~$ {     goto 101
0 a& N5 X- b0 o# W: ]    endif
+ i# f" ?( N# A: c4 C  g" s call gaussj(B1,n,(-1)*gradt)
- N! U# X7 [5 f; z& @  q6 x  C dir=gradt
    x0=golden(x,dir,hessin,b)
" I! C! y  k% U4 u    x1=x+x0*dir
# D& O4 m& X; j; c$ ?) C" L gradt1=matmul(hessin,x1)+b
, V6 U8 U1 Q% q9 E) n6 p s=x1-x
y=gradt1-gradt
call vectorm(gradt,G)
! f1 y5 e" R5 B G1=G
/ @. v1 G8 M" {4 ?. {8 V4 ? call vectorm(y,G)
+ }) f3 f, Z. p3 w, k- e7 b  B" ~9 I* w B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
7 Y) h9 k9 ]' W0 ~1 F" C+ v x=x1
gradt=gradt1
5 n& z: O' C( q  @# E" `% q4 m& P    iter=iter+1
8 }+ s2 s% S7 \! ]  if(iter&gt;10*n)then
    print*,"out"
    goto 101
endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
    goto 100
    contains</P>
4 Z2 s3 l2 ]0 q) M) |3 M<>    !!!子程序，返回函数值   
: Y+ r  P7 J0 k; F( ?( f    function f(x,A,b) result(f_result)
6 q/ d' h- k3 g$ {9 Y+ U0 E. f6 Q    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
    real::f_result
6 N+ k1 ]  s$ t5 e8 m% L0 n    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
    end function f
!!!子程序，矩阵与向量相乘
6 Y7 a. s2 E, w# Y$ ]5 q subroutine vectorm(p,G)
$ d' P6 u3 w, H, ` real,dimension(,intent(in)::p
real,dimension(:,,intent(out)::G
' _) r/ F% K- }6 o n=size(p)
; {" e, U2 e+ A do i=1,n
% ]% k* o7 k; D" }; o+ F    !do j=1,n
7 e" @% q  G# {* u* g. H       G(i,=p(i)*p
    !enddo
+ I5 X* ]* {2 ?. H enddo
end subroutine
  U+ B" u- h* u5 u6 O, a  p
1 R( M! S8 t* e4 S' j+ Y    !!!精确线搜索0.618法子程序 ，返回步长；
) E( S  k5 Y" r( V2 M! r  a    function golden(x,d,A,b) result(golden_n)
    real::golden_n
$ T1 }" ?9 X! D' `8 f    real::x0
  |2 y0 h, E' W; W# T- x    real,dimension(,intent(in)::x,d
    real,dimension(,intent(in)::b
# p+ N. n8 C5 @% u2 g7 ?    real,dimension(:,,intent(in)::A
. }0 O6 O% q) q+ C( }. ^0 z" V" i    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
    f0=f(x+x0*d,A,b)
    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)
: }" x. N5 _, I% ]% w# f' F$ _        if(f2&lt;f0)then
           x1=x0
$ z- {! S6 ~; L        x0=x2
/ f5 R; R& x9 g$ u! u+ E  E8 g        f1=f0
        f0=f2
2 g) O1 x2 r! b* o7 |1 m& U5 k! {1 K        goto 4
* |2 Q% V0 T. U' g1 h        else
           a1=x2
        b1=x1
1 s5 h' X' [% ]( Q: {1 `5 S        endif
    else
2 D& E' y& D( O. R) Z2       dx=dx+dx
        x2=x1+dx
- f  T1 B6 B, ?, }+ [        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
: ~( M0 ^% g0 ?0 {           b1=x2
        a1=x0
        else
           x0=x1
        x1=x2
        f0=f1
. o8 F4 R! F! U- A( [% s. [( Q' Q        f1=f2
        goto 2
9 j6 z: ?* c0 ?# }* w        endif
/ x$ `1 G' i( h    endif
+ j5 e- |+ d; Z" P, R+ F    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
    f1=f(x+x1*d,A,b)
$ q/ v6 ~7 G7 n. y; T' h3 b7 x    f2=f(x+x2*d,A,b)
3   if(abs(b1-a1)&lt;=tol)then
3 Y9 D$ m  h# q1 ?$ p& \, r        x0=(a1+b1)/2
. V$ r+ s& Q, j3 g, Q    else
        if(f1&gt;f2)then
- D! T0 ]; |, f6 b/ a! ^& ]! p7 U* {        a1=x1
        x1=x2
3 B6 j  L. V% G4 c, _9 N        f1=f2
        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
+ F( f' ^+ @5 {" r" F        goto 3
+ F( F) d  S3 N' h% S7 N     else
! E2 A. k/ u5 a4 u+ n        b1=x2
6 T8 i$ ^, O& J' D4 a( h3 O' I& \: A& ]        x2=x1
% @6 `9 C" Y0 |' Q0 d        f2=f1
' f+ ]! z6 T" {6 o8 j        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
8 A. c6 Q1 Z4 t: q; Q        goto 3
     endif
    endif
! `+ w3 X) l$ n9 b! Z9 A7 d    golden_n=x0
    end  function golden</P>
2 {3 V% |+ ?3 G$ o<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
    integer n,nmax
6 t6 U, p9 w& P1 b2 ~    real a(n,n),b(n)
2 ~: e8 z) ]# r; @! h5 i    parameter(nmax=50)
% E2 I( J& [; h6 W2 [  t  k7 t* r1 U    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
$ c* Z: N" F! ^; j1 }6 X& N! Y    real big,dum,pivinv  
    do j=1,n
       ipiv(j)=0
+ @- h' g6 w. ]/ @! w# u1 M    enddo
$ M: z9 o* h, n. y) K    do i=1,n
  k! h% M. s% q1 X3 j0 `8 |/ l7 A       big=0.
* c+ l( K1 u4 o) y. P$ I% }$ l) o       do j=1,n
& R/ T$ U2 o4 q# W       if(ipiv(j)/=1)then
          do k=1,n
! T+ d# g! R5 `% b/ e( Q          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
8 }  [  L0 d7 I           irow=j
) i1 H* K# N' }4 Z# g+ U8 P7 z  ]           icol=k
* c8 j0 H) m( x' n3 d3 i- A! l" `       endif
- w: b0 a/ m& s3 |* d       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
0 j3 G/ r6 ^1 B- ?8 \       endif
       enddo
! c2 v6 ?; `4 ^2 Q9 B* E    endif
. e, C* d2 i$ I+ ]    enddo
- O  c- r7 H+ _& \3 T    ipiv(icol)=ipiv(icol)+1
    if(irow/=icol)then
       do l=1,n
          dum=a(irow,l)
. t  X# E2 Y) H$ R       a(irow,l)=a(icol,l)
       a(icol,l)=dum
       enddo
       dum=b(irow)
       b(irow)=b(icol)
  n$ P8 f$ A0 H3 Y5 }8 b       b(icol)=dum
    endif
$ V2 S! ^2 a' b8 f    indxr(i)=irow
# W5 Y1 u7 j/ u# n' {1 X    indxc(i)=icol
    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
8 u/ k3 \# N/ P5 R3 z0 A    pivinv=1./a(icol,icol)
4 d9 v! `) g3 z1 X7 X    a(icol,icol)=1.
    do l=1,n
* w3 p9 I" g2 _5 y$ u. |  C- e        a(icol,l)=a(icol,l)*pivinv
    enddo
" K3 L$ N7 Q+ I+ T+ U) f1 ]    b(icol)=b(icol)*pivinv
    do ll=1,n
       if(ll/=icol)then
          dum=a(ll,icol)
2 u8 F( a7 i8 ?       a(ll,icol)=0
       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
       b(ll)=b(ll)-b(icol)*dum
       endif
  Z, O  C7 X' V% Y8 X, y. C    enddo
$ m5 Z6 f7 y+ u    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))
- x1 y( o/ p% h) ^* Y% ?       a(k,indxc(l))=dum
6 l! J- I: s& `  X       enddo
' P, c. X2 U% S    endif
% m4 n4 K, Z9 ]  n    enddo
+ `( x% v! [. c5 d, x    end subroutine gaussj
101 end
  [) h6 V2 p: i8 ~4 }+ w  r. a</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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