BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
: @; X3 x- D, F( ?, b9 G    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
* C  v0 o4 ~4 C/ g5 J3 N    !!!dir实型变量，存放搜索方向；
+ X& _# g/ R) o! W9 D    program main
8 V& M/ N8 M7 D    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
, ~/ z8 E* y; |  G5 y    integer::n ,iter,i,j
2 r7 h+ U  e/ ^# Z* z    print*,'请输入变量的维数'
0 D0 G6 b, d! ^( T* o- q6 u# h    read*,n
" |; f+ z  P: a5 |/ L( R/ J    allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),gradt1(n),x1(n))
+ X( ?) s7 y+ p* ]2 z& j    allocate(hessin(n,n),B1(n,n),G(n,n),G1(n,n))
5 U0 N! x" t* O8 t    print*,'请输入初始向量x'
    read*,x
$ g. B5 A: g) I6 y3 n2 g6 h    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
    read*,b
    iter=0
4 A- z8 h7 Y% M- f# \$ s tol=0.00001</P>
' E$ z) b2 W; d<> do i=1,n
; x: j/ X9 H' \) Y6 b+ F    do j=1,n
       if (i==j)then
% N$ }2 w! y# k  f+ i       B1(i,j)=1
& R3 X. G  {& H' ^  Q    else
       B1(i,j)=0
    endif
    enddo
enddo   
" W8 O0 \/ W9 [    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
0 s3 N9 Z5 l8 z: E        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
! F9 N- H9 g/ H& i     goto 101
4 q. p# D9 `. z: ~9 Q    endif
5 f- {0 `- d9 U, e3 ?7 O# h6 }' C call gaussj(B1,n,(-1)*gradt)
$ |6 m! f  E; U) x/ J1 ^) L dir=gradt
3 S- ]9 v& u5 g( }: u    x0=golden(x,dir,hessin,b)
: L% h2 G4 Q: j    x1=x+x0*dir
6 }8 L5 k# X( w' p3 A; m/ f7 K/ R gradt1=matmul(hessin,x1)+b
$ @4 I" P" t$ o s=x1-x
; J2 e1 b& X( P, ^: _! ^- t y=gradt1-gradt
call vectorm(gradt,G)
G1=G
call vectorm(y,G)
4 o# U$ L# l) A! t6 B B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
gradt=gradt1
    iter=iter+1
  if(iter&gt;10*n)then
- ~/ d( ^9 f5 _    print*,"out"
    goto 101
endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
9 Y, f2 O6 R( }2 e    goto 100
7 M# f9 t% V, a/ E+ i- }  V/ o    contains</P>
<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
; c# u% n% j' u    real,dimension(,intent(in)::x,b
2 |/ |$ G  T$ e* L/ t' X    real,dimension(:,,intent(in)::A
    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
/ H* Y" l3 L3 _+ l) X" d    end function f
+ X! V; ?: U% [! L3 w !!!子程序，矩阵与向量相乘
+ G0 T! H# k; t subroutine vectorm(p,G)
real,dimension(,intent(in)::p
real,dimension(:,,intent(out)::G
n=size(p)
do i=1,n
    !do j=1,n
5 \, Z: H. o# H: k       G(i,=p(i)*p
    !enddo
enddo
; W$ k0 Z5 \; @* T; `6 d end subroutine
" v. q. A4 s2 g' ?- W1 a8 J
, J; @2 E6 [' n, {7 U2 r    !!!精确线搜索0.618法子程序 ，返回步长；
' a$ Z2 y/ q0 [5 i) ]" _    function golden(x,d,A,b) result(golden_n)
, E1 s# ?* A! \6 u" \    real::golden_n
6 I, M/ ?4 K# j0 Q) A, B6 A    real::x0
4 x; R" x8 z) u3 R: A1 C" ]% n5 @    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
( b0 J0 f( `1 g- L( D    parameter(r=0.618)
    tol=0.0001
    dx=0.1
" I5 k- p. G6 w9 ^6 I6 `    x0=1
# m* ]6 f9 b" W* T9 Z" e( U& j" d    x1=x0+dx
8 W0 a& ^7 B4 Q/ K    f0=f(x+x0*d,A,b)
+ O: s+ g+ ?8 @% [/ x( d2 [    f1=f(x+x1*d,A,b)
    if(f0&lt;f1)then
4       dx=dx+dx
+ r& M& K# z( i+ V        x2=x0-dx
        f2=f(x+x2*d,A,b)
        if(f2&lt;f0)then
' T* P& D# O0 w: ~& X           x1=x0
        x0=x2
        f1=f0
        f0=f2
7 {8 X7 n' @5 o        goto 4
        else
           a1=x2
; o6 e5 d0 q% f* f        b1=x1
        endif
    else
2       dx=dx+dx
        x2=x1+dx
        f2=f(x+x2*d,A,b)
" @( A  o/ [  n" ^        if(f2&gt;=f1)then
0 x- j8 w; c8 J; m6 T) s           b1=x2
        a1=x0
" v. s3 w$ S) }        else
& D. q4 T! d1 f           x0=x1
        x1=x2
9 o+ T. X3 Q  |* P; k5 O        f0=f1
: n* Q" c- S+ p! S        f1=f2
# @' Q- g9 c9 ~        goto 2
        endif
! [. }7 j" o& @4 k# p& _    endif
" h& P6 ]  ?8 c$ y3 e7 L& q  M    x1=a1+(1-r)*(b1-a1)
    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
        x0=(a1+b1)/2
    else
        if(f1&gt;f2)then
        a1=x1
& I+ ^: R- U, e# ~* B( ^        x1=x2
        f1=f2
        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
        goto 3
, \" S3 C* H9 I9 T     else
        b1=x2
        x2=x1
        f2=f1
! O; b9 c: U- b( j& F* U: t        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
        goto 3
     endif
    endif
    golden_n=x0
    end  function golden</P>
1 `4 v$ D2 H3 [+ m5 Z<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
3 ~4 I/ P- h, C    subroutine gaussj(a,n,b)
, z  y# ~# L  M4 G0 \! n    integer n,nmax
    real a(n,n),b(n)
    parameter(nmax=50)
    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
* A' P- x  ~: I; l$ G& a+ C    real big,dum,pivinv  
    do j=1,n
& f' M7 }8 J: K       ipiv(j)=0
& V+ `# c3 N& E" j; r$ h3 P+ s% M    enddo
    do i=1,n
       big=0.
0 Z0 ^7 P5 }- T  i+ Y2 H9 V& a, a       do j=1,n
       if(ipiv(j)/=1)then
          do k=1,n
2 _. B+ P8 t2 |          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
           irow=j
           icol=k
       endif
       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
& ?* |! B3 p. I       endif
       enddo
    endif
    enddo
    ipiv(icol)=ipiv(icol)+1
- M3 {3 w7 Z' U$ v    if(irow/=icol)then
' D6 x9 O( ?+ U3 h       do l=1,n
          dum=a(irow,l)
       a(irow,l)=a(icol,l)
5 Q# M! ?& Z  `1 w' i2 N6 u7 Q       a(icol,l)=dum
( F$ V0 Q/ q1 N0 A1 Z. C       enddo
       dum=b(irow)
' D' J0 p, y2 k( ~8 S       b(irow)=b(icol)
( g" ^! C; P  v% w0 @: F       b(icol)=dum
    endif
* s; D* Y; h8 c7 a    indxr(i)=irow
6 ~8 B# `. [2 v0 B; x    indxc(i)=icol
3 }' |0 o+ W( g* O    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
    enddo
    b(icol)=b(icol)*pivinv
7 U+ t7 T8 v) A    do ll=1,n
       if(ll/=icol)then
1 F4 a& ^' c3 ?5 K2 _          dum=a(ll,icol)
/ N$ f2 l0 {1 U0 Q6 u* H       a(ll,icol)=0
; {, D0 R6 ~& Z' z8 J  N* ?       do l=1,n
% |# a6 H2 O) H, o          a(ll,l)=a(ll,l)-a(icol,l)*dum
8 i5 \  S1 Y( y$ ?- ~       enddo
       b(ll)=b(ll)-b(icol)*dum
9 A# B" E$ D, _. `, T6 d/ x- H       endif
    enddo
$ K, Y; `& m9 e1 n( v    enddo
    do l=n,1,-1
- v; j/ i2 e5 A+ Z8 k2 Z       if(indxr(l)/=indxc(l))then
       do k=1,n
          dum=a(k,indxr(l))
1 R- u! F- f" O# ~       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
       enddo
    endif
' I9 f! C3 E% ?- ?. \4 I) F/ }6 n    enddo
    end subroutine gaussj
101 end
</P>
2 N  e4 K# @3 Q& x5 O<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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