BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
" }4 o" b$ g1 b) z+ C* N9 l    integer::n ,iter,i,j
  j6 K; b  Q9 A  |0 _" g$ K2 `    print*,'请输入变量的维数'
% N2 n4 ~1 L$ U, F( W7 E, |) U% L    read*,n
0 E1 d( a" M% X& O8 ?+ x- D    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))
6 `. B5 k; j  u/ U6 B; z. j    print*,'请输入初始向量x'
7 G, c7 h# D5 Z, x    read*,x
    print*,'请输入hessin矩阵'
5 B* M8 b4 Y$ k    read*,hessin
    print*,'请输入矩阵b'
    read*,b
- i9 i# h  z0 Z3 g2 l0 J: g& n. {  W/ u    iter=0
tol=0.00001</P>
, a, h6 g& `* F. U1 U9 f* X' C<> do i=1,n
    do j=1,n
2 A6 Y) `4 C. v# V       if (i==j)then
       B1(i,j)=1
! a$ g0 I% ^. K( l. Z4 L    else
       B1(i,j)=0
    endif
    enddo
enddo   
; u) u8 e; O1 U% ]- y+ ?    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
0 O/ E# e1 N/ [6 q0 q     goto 101
# d" T6 m1 L+ ^4 T2 A- L) F3 n    endif
  O" y+ T2 {9 s1 f6 s7 X- c call gaussj(B1,n,(-1)*gradt)
/ f! G' V4 O; ?$ M9 P% H3 J  v dir=gradt
5 }* E/ m9 o% j5 M8 g& S! I  M    x0=golden(x,dir,hessin,b)
$ B7 B4 O- P; J    x1=x+x0*dir
4 g# y/ ]$ y& {* T- L gradt1=matmul(hessin,x1)+b
s=x1-x
9 d3 ]! Y- r" Q4 g y=gradt1-gradt
call vectorm(gradt,G)
1 W" L+ Y: b# r  |: }6 l5 t; _ G1=G
6 B: x) w6 ~/ A% J6 H; m call vectorm(y,G)
1 z. [& U, `' T6 p& X2 h) v# ~0 g B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
gradt=gradt1
" Z2 N2 p$ H5 ~# @' x9 X9 a- a    iter=iter+1
  if(iter&gt;10*n)then
' J* E" o( ?: e* K4 t7 r    print*,"out"
    goto 101
7 C& O. O7 ^  T. a/ Q/ W endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
    goto 100
    contains</P>
<>    !!!子程序，返回函数值   
# i2 D, M2 a9 \3 M    function f(x,A,b) result(f_result)
    real,dimension(,intent(in)::x,b
/ w  {- J0 f6 O  o# Y0 ~    real,dimension(:,,intent(in)::A
    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
: B  [. N5 y4 {% b+ B5 x& a3 `    end function f
0 |# O6 i4 P0 F8 s2 `7 q* u !!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
real,dimension(,intent(in)::p
7 t" \& M! i; l& t) L- c9 \" ? real,dimension(:,,intent(out)::G
n=size(p)
do i=1,n
    !do j=1,n
       G(i,=p(i)*p
, {, U, g  m' X) m- G    !enddo
* r4 ^" T$ `- n2 T: r enddo
end subroutine

    !!!精确线搜索0.618法子程序 ，返回步长；
    function golden(x,d,A,b) result(golden_n)
5 O+ ^6 N8 T; a$ G/ R* K* o* w    real::golden_n
! {1 d, m8 O& {" w- ^    real::x0
    real,dimension(,intent(in)::x,d
8 Z/ x5 m( Z' j) f( R    real,dimension(,intent(in)::b
    real,dimension(:,,intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
4 s% g) t1 M( I1 h' `( _    parameter(r=0.618)
3 n8 Q; ?" W$ a7 P2 g    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)
        if(f2&lt;f0)then
+ L" v0 p; W! |! C/ Q$ d" c+ z           x1=x0
* [2 T& A1 U3 N        x0=x2
* ?% [  x" w9 u) v( z        f1=f0
        f0=f2
        goto 4
+ f, R5 ^# F9 I  ?        else
           a1=x2
: e( o1 e1 H+ I! d4 d, z$ k, t+ F        b1=x1
        endif
5 C# O" k8 \% w9 a3 x    else
2       dx=dx+dx
- I9 S4 p& `; Y  ]2 f        x2=x1+dx
        f2=f(x+x2*d,A,b)
) r5 F5 d! b- H" _        if(f2&gt;=f1)then
           b1=x2
        a1=x0
        else
" _/ g. e0 v6 ?! l           x0=x1
8 D: N+ q  ?( A5 K0 i        x1=x2
0 T4 I2 a! ~- v& w        f0=f1
        f1=f2
0 B/ k' ?* X  x- M) E/ u        goto 2
        endif
    endif
- \% J$ H# q6 a% Y6 w    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
9 V! s- [' q* s( X$ X    f1=f(x+x1*d,A,b)
1 c+ D! s2 `/ _  X9 v    f2=f(x+x2*d,A,b)
* e& M& l8 p/ c! x; H3   if(abs(b1-a1)&lt;=tol)then
; _5 S7 a1 V3 D. }        x0=(a1+b1)/2
    else
        if(f1&gt;f2)then
! _# w. a5 N( s- C3 U        a1=x1
: Z( F4 Z' D1 s        x1=x2
0 l; z- {8 c) h( L        f1=f2
5 t: ?8 i( h: P; H3 Z0 S* A* X        x2=a1+r*(b1-a1)
& A; G! c( O  X8 R2 w        f2=f(x+x2*d,A,b)
        goto 3
     else
; @+ E6 R% w! h1 V9 Y, F, w5 [        b1=x2
6 N+ w) S. R) F5 t        x2=x1
        f2=f1
  Y* N  M# b* {  q$ }/ d& Y        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
& I! t. A, `" A% N        goto 3
     endif
, i* u7 L1 Y2 z5 \    endif
    golden_n=x0
    end  function golden</P>
<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
/ c- [( L2 k5 u# c, a% J5 ^. K- y    subroutine gaussj(a,n,b)
    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)
8 Y, Q' _& R) C3 \3 f) P    real big,dum,pivinv  
0 b, J0 c4 ]- o5 A, k    do j=1,n
       ipiv(j)=0
    enddo
0 I' v$ e/ L$ f6 `1 ?( `    do i=1,n
       big=0.
       do j=1,n
       if(ipiv(j)/=1)then
          do k=1,n
          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
& a% D; K5 S' M1 a           big=abs(a(j,k))
           irow=j
           icol=k
. w7 r% \* @8 }: m: K       endif
       else if(ipiv(k)&gt;1)then
! t' ?. @, t8 g6 N  d" ?) q$ C' h          pause'singular matrix in gaussj'
       endif
/ i2 {+ A- r4 M7 r) b       enddo
& {, h( G9 M5 o6 r" J    endif
, g- I  e  d+ v% D" H# H    enddo
    ipiv(icol)=ipiv(icol)+1
2 v3 s. k% G/ L8 _, y# D+ G    if(irow/=icol)then
0 ^( p5 Q  @/ S4 {2 F5 V+ W& z2 s       do l=1,n
          dum=a(irow,l)
8 v1 c$ R2 B& |; T- ~1 t       a(irow,l)=a(icol,l)
       a(icol,l)=dum
       enddo
# Z7 i6 O: Z) o4 i$ X% R) e6 g       dum=b(irow)
       b(irow)=b(icol)
       b(icol)=dum
6 E# C8 i) b" I1 P% B$ s& I    endif
7 F9 E3 o& u, H* X    indxr(i)=irow
    indxc(i)=icol
    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
8 V8 f  Y6 o) I4 T  u  M    enddo
# d% M0 c, H8 g1 E    b(icol)=b(icol)*pivinv
    do ll=1,n
       if(ll/=icol)then
1 ~$ ]( d+ b& {' h. ?8 d# C3 r: B9 t          dum=a(ll,icol)
       a(ll,icol)=0
       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
7 }4 P1 G7 o7 @% y. n* [+ H       enddo
5 ]1 ^" {  d/ r: n: R+ h       b(ll)=b(ll)-b(icol)*dum
9 w/ _. a% P( J       endif
' Y8 b1 ~7 M# |    enddo
    enddo
/ p  O0 A) v! n+ f3 s  t    do l=n,1,-1
$ @  X. M0 x3 p# @6 k- q. Q       if(indxr(l)/=indxc(l))then
       do k=1,n
: x( h8 ]$ n7 t/ V, l' @/ M          dum=a(k,indxr(l))
) ^0 c% I" ?; d4 ^* f3 |, ^3 X7 c( `       a(k,indxr(l))=a(k,indxc(l))
# i7 G2 n, w7 I# M/ N8 Y       a(k,indxc(l))=dum
: B3 b& ~8 `5 d% `       enddo
" j  v5 U- `8 J+ O  n    endif
& H: Z2 h8 ^. y( i    enddo
    end subroutine gaussj
1 p1 |$ G' {, R2 `101 end
6 ]# w9 a/ Y3 j( P</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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