BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
3 R% N  ~0 F$ {9 m3 o    !!!iter整型变量，存放迭代次数；
0 c, Z6 J7 J5 |4 ?    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
8 M! O( e  y& \: G$ ]    !!!dir实型变量，存放搜索方向；
7 A' g) z6 `$ u1 f9 U3 d    program main
( q8 Y9 c$ c) e8 N2 b% O  c    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
    integer::n ,iter,i,j
    print*,'请输入变量的维数'
# N0 w: X! K" b4 L* D( f* K    read*,n
    allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),gradt1(n),x1(n))
* G, s' e' O' S( W* p) J    allocate(hessin(n,n),B1(n,n),G(n,n),G1(n,n))
, ~. S! K; N: m. {& x    print*,'请输入初始向量x'
' C4 ?3 H1 S5 W+ u    read*,x
    print*,'请输入hessin矩阵'
8 z8 [% s2 ^$ k3 X) j. p6 z    read*,hessin
    print*,'请输入矩阵b'
    read*,b
$ k6 J0 A' i$ L$ n; ?: I    iter=0
tol=0.00001</P>
* }: h$ L! a! k7 |* L/ j% a<> do i=1,n
    do j=1,n
       if (i==j)then
( p- D, _3 p2 x4 d" y# |1 Q9 y       B1(i,j)=1
- B! i) ^/ k6 h2 J    else
) \! |0 i# Y8 W2 ?: h0 w       B1(i,j)=0
    endif
; S0 d" ~& R6 [  G' }    enddo
# J" }0 l; C9 y' R enddo   
    gradt=matmul(hessin,x)+b
1 D% X8 C7 i/ S9 l( J' u, N; i100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
        !print*,'极小值点为：',x
1 ?; a9 _0 D" C8 y     !print*,'迭代次数:',iter
+ A- F! }6 J( H6 v0 d+ q     goto 101
    endif
% Z/ Q9 F* A$ r$ W3 V call gaussj(B1,n,(-1)*gradt)
dir=gradt
5 d7 `# a% M0 |7 q3 r    x0=golden(x,dir,hessin,b)
' d8 Y. A2 P) H5 H    x1=x+x0*dir
' M( j+ @- Z9 {/ Q gradt1=matmul(hessin,x1)+b
s=x1-x
# |; ~) x3 Z3 h: c0 w3 L. Z1 k# ` y=gradt1-gradt
, G' s3 D9 m+ ?# s* \4 N3 Q call vectorm(gradt,G)
G1=G
: `0 h  }( m/ s. U" @ call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
& }) Y7 T' t. u" \" Y" E' p gradt=gradt1
; W0 h' C9 h0 v    iter=iter+1
  if(iter&gt;10*n)then
; ^/ r! P  R6 A! Q" a: Y1 o0 s    print*,"out"
    goto 101
# t% v. i9 U: f5 a endif
# A$ e& N. k6 o& T4 s* K    print*,"第",iter,"次运行结果为",x
( \& [1 G9 H, m3 A; f6 z! g4 q print*,"方向为",dir  
    goto 100
! a+ W! O; m9 q6 i1 u, U$ [' E    contains</P>
<>    !!!子程序，返回函数值   
$ ?5 x: a1 q+ f/ n3 X0 }+ t( ?4 O    function f(x,A,b) result(f_result)
+ s7 }& A. ^, n2 p# r    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
) U  O/ s6 Y5 t- F- H    real::f_result
+ g7 K# s- j# z8 z1 l( P    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
    end function f
6 k) O* ]) c, `1 y/ a4 \ !!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
- \: N/ _" }* I) h1 p# j real,dimension(,intent(in)::p
# l6 W# Z, T8 @ real,dimension(:,,intent(out)::G
n=size(p)
do i=1,n
    !do j=1,n
       G(i,=p(i)*p
, k+ N+ X; J8 q' A    !enddo
. C( z0 {; }: D enddo
, \9 E) A& y* }+ O+ f6 _ end subroutine
3 g' Z' d( Y' _$ ^
- U, M& G3 h- Z: G& y6 E    !!!精确线搜索0.618法子程序 ，返回步长；
    function golden(x,d,A,b) result(golden_n)
    real::golden_n
9 F( _6 E6 ?* c2 N! Y    real::x0
    real,dimension(,intent(in)::x,d
    real,dimension(,intent(in)::b
    real,dimension(:,,intent(in)::A
) i2 S1 v  m* q$ U" Y& u) ]    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
    parameter(r=0.618)
    tol=0.0001
  G0 Y1 ]2 J: u; D7 |    dx=0.1
    x0=1
0 b: Y5 G+ T9 S" V    x1=x0+dx
4 _! l! F" F& o4 f    f0=f(x+x0*d,A,b)
/ @9 b; x/ O- t    f1=f(x+x1*d,A,b)
" G+ }2 G( L& m" z. X, P    if(f0&lt;f1)then
8 u( _1 V# t. h% Q4       dx=dx+dx
! F. m2 r- v# W. S* Z% ?        x2=x0-dx
, z2 m$ C6 S3 [" K+ P% f3 j" a        f2=f(x+x2*d,A,b)
; X1 ^1 {2 [9 Z' H+ \$ G        if(f2&lt;f0)then
. Y' x: p% w0 K$ C' a3 M' h- X$ l           x1=x0
. Z; O" n. }0 ~2 |+ p  I0 @/ F        x0=x2
        f1=f0
        f0=f2
        goto 4
        else
           a1=x2
        b1=x1
9 p5 k1 W3 g( V) w) k$ _( y        endif
3 p1 P: R4 p# D5 l( }    else
2       dx=dx+dx
( `3 A% ]8 E5 B& A$ Y( m        x2=x1+dx
        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
           b1=x2
. X$ c, C2 \! Z        a1=x0
& {' ~* Z- j5 v6 Q* q; Z# F9 ?        else
           x0=x1
# B8 T- S; a* d$ `        x1=x2
        f0=f1
        f1=f2
. F; o+ F) R. T! g, U# _3 \        goto 2
8 _0 l) V, G& }1 K8 }  R" m        endif
5 g: T/ {2 L; m; c* w2 Y    endif
' O6 O5 j) K. M4 x9 z, Q4 C    x1=a1+(1-r)*(b1-a1)
, Q5 W- F" Q) Z    x2=a1+r*(b1-a1)
    f1=f(x+x1*d,A,b)
$ ~' X" I2 @) ~( {' U6 c    f2=f(x+x2*d,A,b)
7 o  V- _/ W% k+ n7 }  B+ A3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
5 K! p- X) \3 f3 j, w0 P    else
        if(f1&gt;f2)then
$ U+ x; X# b/ t+ |5 E8 ~& |        a1=x1
. {; k2 e3 h3 q2 ~        x1=x2
" ]- [8 @* T/ q6 p: \        f1=f2
5 z* n+ O- K) E2 |        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
        goto 3
# M; v0 J' x! p$ |     else
- }# x0 W( u: i1 O% I9 f, N        b1=x2
7 w1 c; ]& T" G; A# Q        x2=x1
7 s- d" u: j) Z9 y" S+ x        f2=f1
7 h9 V# Y; N6 g. i# E        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
        goto 3
% {" i! I4 v+ t2 N% B- E# `7 L     endif
    endif
    golden_n=x0
    end  function golden</P>
( `# p: I9 S( P! A7 Q<>
. p5 S: \* u4 p  [" p6 ^    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
. |$ c' o9 o+ |0 M* I6 B    subroutine gaussj(a,n,b)
  D8 I& L/ ^" c1 Q- u    integer n,nmax
    real a(n,n),b(n)
8 T8 q. p6 C  K9 [. O    parameter(nmax=50)
    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
    real big,dum,pivinv  
8 |2 ^4 j: r0 j' M    do j=1,n
* k5 q' n; Y  Y% g) O, y       ipiv(j)=0
. j# k' h! X" b    enddo
7 h2 K: K- u( D  o+ }6 c    do i=1,n
: L" D( m% w& l0 f4 L' a       big=0.
8 B6 y0 h- A4 \9 P( k% ]       do j=1,n
* l6 s2 O& W5 l+ l1 i       if(ipiv(j)/=1)then
# L8 T8 p8 A$ Z6 m! q) _( k# f          do k=1,n
+ x/ [( K9 u. f# c; F& e          if(ipiv(k)==0)then
9 m" D! `* J1 L; y7 P$ j          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
9 @  Z4 J* a( M: C5 p# v) n6 @* y6 ^           irow=j
           icol=k
       endif
       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
       endif
       enddo
* C  i' h# {5 W8 q3 b    endif
    enddo
    ipiv(icol)=ipiv(icol)+1
9 J2 P1 m6 Q$ z6 Y, G6 `! L    if(irow/=icol)then
. |. z' s, e$ l1 O1 |9 a- f0 U       do l=1,n
          dum=a(irow,l)
, {$ m6 ]+ w8 V, O) V; H& S       a(irow,l)=a(icol,l)
& U5 B1 e! ^3 }! L/ e5 R* @; x7 o       a(icol,l)=dum
       enddo
* S2 F  [, c7 g. n5 m( U       dum=b(irow)
       b(irow)=b(icol)
0 \$ w6 {: r" h- ]. _       b(icol)=dum
    endif
    indxr(i)=irow
, ?2 ]8 X: }# o7 H8 B+ v    indxc(i)=icol
3 d$ S9 w& D0 n5 L8 v1 b    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
( q- C+ j+ W. @3 ?! s2 ]    a(icol,icol)=1.
+ X) N9 P5 ~* D. J8 u! V( u    do l=1,n
* Z. t% X/ I( _' o+ h+ D6 o        a(icol,l)=a(icol,l)*pivinv
; u/ J. a  n5 f    enddo
" V$ O4 A4 ]; I6 e) \* y    b(icol)=b(icol)*pivinv
: ]3 l& C: s) s2 [* A$ |" b' a    do ll=1,n
3 ~" M6 e( f1 f# i       if(ll/=icol)then
          dum=a(ll,icol)
       a(ll,icol)=0
7 e' R' Y# W' E6 r% T       do l=1,n
8 }/ k) P6 A7 w# f          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
9 o2 _; t  H9 a2 n' I" |1 m       b(ll)=b(ll)-b(icol)*dum
3 [; W" x6 F+ v5 q$ b4 k8 H4 {       endif
% o) M9 i1 n0 ^7 r: Y. W    enddo
    enddo
    do l=n,1,-1
. B: G+ q$ ?$ D. i6 L( F' g       if(indxr(l)/=indxc(l))then
% I1 _8 E5 I/ x3 K, }$ R       do k=1,n
          dum=a(k,indxr(l))
$ ?2 E5 S( f/ W$ d+ H- q3 i+ e6 p       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
       enddo
8 m. z7 H0 X; c% D    endif
# m8 k' a) L% F* h. [1 z    enddo
    end subroutine gaussj
. u& Z+ C, r( d, F101 end
</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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