BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
, ^0 L' g2 ]* q/ Q  B    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
1 H8 S0 y* @7 [+ n: |    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
    integer::n ,iter,i,j
    print*,'请输入变量的维数'
    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))
; c8 w* T& ^8 v2 L* h    print*,'请输入初始向量x'
- M. Z3 L8 q& [- i2 n    read*,x
* N7 x1 v6 e6 X    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
    read*,b
    iter=0
; I2 p' p' l( n/ @5 R# A/ {; V tol=0.00001</P>
5 W; J8 m: g$ X0 I6 b$ p4 T<> do i=1,n
    do j=1,n
$ p+ s( X1 R3 K6 U       if (i==j)then
       B1(i,j)=1
6 `& x; l9 G" `! M, H! t& ^. Q    else
5 P; y; [) C. y3 o( E2 Z: T       B1(i,j)=0
    endif
; f5 f! ]: M2 ]$ a1 `    enddo
enddo   
0 \8 B# d+ |7 G( Z: s8 W& T    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
- K- W4 v, Y0 f8 T     goto 101
4 l9 M! |/ l9 i    endif
4 l4 S4 I8 A; }, W7 o5 _" z call gaussj(B1,n,(-1)*gradt)
4 M0 Y/ T. d7 j. J2 | dir=gradt
    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
+ h% f$ C: g' \- a" M, b/ N# N s=x1-x
% W7 s) v- ?: ~; U8 v7 n2 P y=gradt1-gradt
7 U" U% @7 r+ a9 d call vectorm(gradt,G)
G1=G
6 U7 `3 W. S. _  X  l2 H& s call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
9 i3 `, t* D2 |3 g8 ^/ | x=x1
gradt=gradt1
    iter=iter+1
3 w3 j3 j8 j" ~% m" r  if(iter&gt;10*n)then
/ ?* m' }# p: M& Z+ e3 i9 Z    print*,"out"
" |, [; W8 A. R& q8 ?    goto 101
  e  u1 }0 V) e5 C endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
    goto 100
    contains</P>
<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
" j* M, A7 v5 Q, W9 j    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
1 W# G( y1 u# t' x& Y3 E    end function f
$ e$ l, i: J) S5 o% G  w( N' P !!!子程序，矩阵与向量相乘
- y4 B$ y( o7 D1 ? 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
       G(i,=p(i)*p
3 @) @5 A! N& H- I8 u# X* B/ D    !enddo
& f: L2 f! f! I8 q! A9 n enddo
end subroutine

. `0 M2 W' S9 g    !!!精确线搜索0.618法子程序 ，返回步长；
' Y" X$ z7 Y9 ~. R# z    function golden(x,d,A,b) result(golden_n)
    real::golden_n
! V% U- T$ O$ F9 K9 P% [; _    real::x0
    real,dimension(,intent(in)::x,d
: Q; b: }! E( f- f( }' I- p    real,dimension(,intent(in)::b
' K9 w" N+ D( r6 G' [; V    real,dimension(:,,intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
1 f" B$ d$ N* L. _; s    parameter(r=0.618)
" y; u+ Z0 E; T$ k$ ^% }+ g, \    tol=0.0001
    dx=0.1
    x0=1
    x1=x0+dx
* Z8 R2 Z1 X9 V8 e    f0=f(x+x0*d,A,b)
! \$ {! h5 i" k8 F8 o* S; k    f1=f(x+x1*d,A,b)
0 s$ @; y6 T/ `" z9 @6 ^! k    if(f0&lt;f1)then
4       dx=dx+dx
" S( i% z- V7 b- O        x2=x0-dx
        f2=f(x+x2*d,A,b)
        if(f2&lt;f0)then
           x1=x0
        x0=x2
# U) x, a$ p4 E3 N4 d; L$ `6 A/ i        f1=f0
% Y' }3 Z4 R' I' v9 s/ r# B        f0=f2
        goto 4
* v  |4 c- {1 s) I# G# o        else
           a1=x2
  I! [! s: o0 E# s. q3 D9 ?        b1=x1
        endif
+ W" e% ^/ w0 \* ]. T( {    else
9 {& T1 D) |/ Q  j! F/ k6 b2       dx=dx+dx
+ C* H7 D6 s! `7 l        x2=x1+dx
$ }1 j9 \+ Z8 @! s- K        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
           b1=x2
$ h* W2 z# T+ x- b        a1=x0
7 _# y2 t4 x% v: d/ V. m5 x        else
2 R6 B; e( K3 V6 @  g; }' K' p$ {           x0=x1
        x1=x2
        f0=f1
# Z& ?2 t7 m! Z) P        f1=f2
' y1 @9 m$ R8 |9 M2 X) G        goto 2
% P+ \. }2 T! K3 x1 S        endif
    endif
    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
- S% `$ v' v0 I( V! v4 L( ^. i9 w    f1=f(x+x1*d,A,b)
+ c0 M3 S8 M" ?    f2=f(x+x2*d,A,b)
" Z$ J/ f* f1 M$ a- S' s3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
    else
- i4 L6 x- N# c. T) t- ~' J        if(f1&gt;f2)then
        a1=x1
        x1=x2
& S# [9 e# W* ]  ?& \, n5 d( o        f1=f2
+ h( \  m' J8 S9 [. y) }3 S" D* K        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
        goto 3
     else
$ e. Y! e& f# L. k% h+ Z* J        b1=x2
        x2=x1
( ~5 A$ K( m8 R* I0 C1 T        f2=f1
3 x: U2 }+ J& @) |# ^, H% `. y7 {, c        x1=a1+(1-r)*(b1-a1)
" L# A" I9 J3 V        f1=f(x+x1*d,A,b)
        goto 3
     endif
# _) q* D: `# o# b+ a    endif
2 h+ D& x( m4 \" r$ o    golden_n=x0
    end  function golden</P>
) _  U0 V8 C+ g0 t<>
- a. a2 z/ R. ~; w- l0 I    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
" e- u( S- s- S9 a" {' S    integer n,nmax
9 o1 b' E; o! @4 c" r    real a(n,n),b(n)
; J$ Z0 X6 Q1 I' q% J& J$ [1 F' ]    parameter(nmax=50)
+ b6 L' G4 m6 U# }1 k% Z9 Z$ c    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
    real big,dum,pivinv  
& F$ S. L$ W1 R, o2 A4 k7 u    do j=1,n
       ipiv(j)=0
    enddo
    do i=1,n
       big=0.
8 }1 [$ r2 T6 j       do j=1,n
       if(ipiv(j)/=1)then
          do k=1,n
4 i! c2 \5 i, a  m9 L, y          if(ipiv(k)==0)then
: G. Y0 o6 i( S% u( r          if(abs(a(j,k))&gt;=big)then
% M, T& d' D% R           big=abs(a(j,k))
! M7 s$ i; x! R8 L           irow=j
           icol=k
       endif
       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
       endif
7 G6 k" O, b0 }* s  k* L( z       enddo
6 k1 c/ ?9 \/ M! Z( c2 f    endif
% u' c4 X; \% u$ J) H  x    enddo
- U+ N, y( [3 Z    ipiv(icol)=ipiv(icol)+1
    if(irow/=icol)then
; ]" }  y3 [; z       do l=1,n
. B, J. z: t5 J, i          dum=a(irow,l)
2 H5 ^( G1 i  s) U1 G% ~       a(irow,l)=a(icol,l)
' v$ W# E" @$ y% _) U" ^! y, p5 ]4 |       a(icol,l)=dum
       enddo
. y- n" c/ P6 Z. {8 ^9 A       dum=b(irow)
' U- i4 r0 v, e" ?+ T  Y       b(irow)=b(icol)
       b(icol)=dum
# W, D8 G# a  r; d    endif
    indxr(i)=irow
4 h+ p1 M) _% V3 Z4 D    indxc(i)=icol
' r/ W" B; ?) G2 }2 h  @( {. M4 v    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
  b# ]9 H% b/ X  |    do l=1,n
0 S/ a1 G( P( ?9 f        a(icol,l)=a(icol,l)*pivinv
# j" p- i* p% e9 l# g" T' w    enddo
' n( O. i1 v6 @* X" l    b(icol)=b(icol)*pivinv
0 [: ?8 Y7 K9 x& e    do ll=1,n
: R6 f. N! H- Q6 k       if(ll/=icol)then
: C1 l! g* H5 O- @0 Q          dum=a(ll,icol)
       a(ll,icol)=0
. i: G7 b7 f' ]       do l=1,n
5 [# {9 ?! ^/ P: I          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
. V; H) d! P0 H0 ~3 ^  ^       b(ll)=b(ll)-b(icol)*dum
       endif
. u+ e7 ^& o7 S3 U+ z    enddo
. Q- R* l; _1 {% L) o6 u. K    enddo
1 v: C" q0 C  ~- X* M7 @    do l=n,1,-1
. x5 L. U5 j6 X" Q       if(indxr(l)/=indxc(l))then
       do k=1,n
4 o. B  J! m" r1 m7 v: r( ?/ A          dum=a(k,indxr(l))
. {# {* m: m$ L. _! [       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
7 U: ]& G. ~' I% h* c5 d& r: D       enddo
) u* [! ~- C+ t3 e0 Y% P. G1 T, v! C    endif
) F( O/ R, {! K8 [5 t7 ?    enddo
8 t( `) ^6 `5 r! X; @    end subroutine gaussj
8 z/ I+ i7 d' _7 Z101 end
</P>
) b" E. Z3 @5 n' g# Y8 x<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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