BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
+ U$ ~' X' U4 b. h; m  @    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
' l  |9 \9 B) [+ F( `8 t2 F    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
    integer::n ,iter,i,j
    print*,'请输入变量的维数'
0 C' B  L( \. M$ W7 r. a* K/ p. b. t( ~    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))
    print*,'请输入初始向量x'
, h7 G6 J- v6 Z    read*,x
    print*,'请输入hessin矩阵'
6 b- [* ?, I! t0 a$ z, p! F9 V- `7 v+ `    read*,hessin
    print*,'请输入矩阵b'
    read*,b
    iter=0
tol=0.00001</P>
' z: h: {* s8 R( Y: t3 X4 G  n<> do i=1,n
8 T6 W: w6 \% g& s# `" P    do j=1,n
) _* ~4 O' ^6 W5 V  l- h       if (i==j)then
       B1(i,j)=1
  P' ^+ S8 t1 l/ }/ q    else
% _8 o  A8 }+ [5 W6 _       B1(i,j)=0
    endif
3 ~2 n. ]/ |5 N! [    enddo
enddo   
    gradt=matmul(hessin,x)+b
4 L) t2 C! f: v/ N  K. C100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
6 H0 D& G4 s, |9 |% Q' ]1 [        !print*,'极小值点为：',x
# O* S) K2 P4 f     !print*,'迭代次数:',iter
     goto 101
    endif
6 z" ^& T+ S  U% n# q call gaussj(B1,n,(-1)*gradt)
dir=gradt
    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
8 R1 g" {: X2 z4 d2 k% X8 G s=x1-x
. p! ^  t% u; _. h( O y=gradt1-gradt
call vectorm(gradt,G)
" ?% [- ~) J% A* Q  E G1=G
* Y5 p4 z$ i8 ^2 }3 z3 K* p call vectorm(y,G)
' y; {  T- }! {4 z6 f6 D B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
% r6 {+ O2 [- L7 ]' d x=x1
gradt=gradt1
    iter=iter+1
  if(iter&gt;10*n)then
    print*,"out"
6 h" F6 w0 T+ @6 f. B5 v    goto 101
; f1 S) r- N1 p" D! E' S& o endif
( p0 \) ?- `4 @. V+ x  o    print*,"第",iter,"次运行结果为",x
9 K9 f, z: F' Y2 y print*,"方向为",dir  
3 s6 a2 y5 t  n# O6 k" ~: F0 Z    goto 100
    contains</P>
) P! ~) h7 [9 z9 H) j; _8 [<>    !!!子程序，返回函数值   
/ w' E+ P7 @# L! |8 }/ p9 C, t    function f(x,A,b) result(f_result)
    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
    real::f_result
0 D0 N( P" Q, j6 x) n    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
2 D$ Q% p' H- T7 v7 ]$ g    end function f
!!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
6 ]  S# U1 m5 e8 {# J1 h4 U$ E real,dimension(,intent(in)::p
4 b$ g/ P6 g. K! v) [; c1 s) C8 K real,dimension(:,,intent(out)::G
n=size(p)
) ~$ I8 _* U" P6 r* T do i=1,n
    !do j=1,n
       G(i,=p(i)*p
$ U' V9 h" }) o  s* W  r& h8 X    !enddo
enddo
; i0 L, {0 _) k0 t2 l end subroutine
# D, x5 R. _# f* I# I
    !!!精确线搜索0.618法子程序 ，返回步长；
0 y3 p2 o7 X1 j    function golden(x,d,A,b) result(golden_n)
    real::golden_n
0 C8 k; w" N/ |; M) M% A* C% W    real::x0
- w0 \+ P5 \, o" s" ]3 a    real,dimension(,intent(in)::x,d
  [1 P, Q0 G2 c3 u    real,dimension(,intent(in)::b
7 W3 C* S8 ~- }/ F+ L$ O  r    real,dimension(:,,intent(in)::A
9 u- l. o# w" ?0 l    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
    parameter(r=0.618)
5 ~1 }2 z+ s5 Z    tol=0.0001
    dx=0.1
    x0=1
2 w0 o8 T% K& t- @% h  F/ x/ M    x1=x0+dx
    f0=f(x+x0*d,A,b)
' \2 Z( d( ?0 R" \" {& S    f1=f(x+x1*d,A,b)
    if(f0&lt;f1)then
4       dx=dx+dx
        x2=x0-dx
4 Z: N- M& z  E! i        f2=f(x+x2*d,A,b)
        if(f2&lt;f0)then
           x1=x0
        x0=x2
        f1=f0
6 R+ e- i/ C4 I- z$ ~& `9 g        f0=f2
4 P! [, E' u( L        goto 4
1 y. O4 T2 g, S1 f4 ]. P        else
           a1=x2
2 M- I0 w# \2 s6 V) I        b1=x1
5 C8 J4 z( s; c! m9 W9 `$ v$ M        endif
    else
2       dx=dx+dx
' J( P" P6 u% e, h4 N* y% b) c; G: [        x2=x1+dx
        f2=f(x+x2*d,A,b)
: [8 _. A. p) X1 Y  }        if(f2&gt;=f1)then
- }$ E' j/ I0 z6 Q$ g$ X           b1=x2
9 f4 p* ~9 Q9 j7 w1 R        a1=x0
        else
  R: Z) J! m) C" s0 F# b           x0=x1
/ a" t# z( g0 w2 w/ v        x1=x2
        f0=f1
        f1=f2
& w$ X, X/ K0 U; g' F9 Q0 ]        goto 2
        endif
' }$ C3 i4 g5 X* l) V9 B    endif
" ~2 h: `& b( X' F    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
    f1=f(x+x1*d,A,b)
' K# i2 X- N+ p& ^' t    f2=f(x+x2*d,A,b)
3   if(abs(b1-a1)&lt;=tol)then
0 @+ V) c% [8 F+ K" ~        x0=(a1+b1)/2
    else
) J  ?* k, d% H        if(f1&gt;f2)then
- j( E# B& l* C! E$ A        a1=x1
        x1=x2
        f1=f2
7 ^4 q4 V2 Z& b  {& |1 U        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
        goto 3
8 W  m) |3 e2 A0 G; p# T5 ]     else
/ c$ N- o0 X" F7 L) r5 W- z        b1=x2
! q0 H" X8 [: J+ H1 f( t. F+ {        x2=x1
6 z" O6 H1 q: |5 R$ A        f2=f1
        x1=a1+(1-r)*(b1-a1)
% c, G7 q1 M$ g, C; t; M/ Z        f1=f(x+x1*d,A,b)
( h! a% M2 {4 s& E, m4 W! Q9 }        goto 3
7 o8 o; B/ `4 J+ h     endif
9 K9 K9 d# o* O* _/ G    endif
    golden_n=x0
    end  function golden</P>
5 N1 d2 r9 b& @: e" R' Z/ D- \<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
! w; W% Q9 L! S* P    integer n,nmax
! x; I$ y; ^3 }( }! t; N7 Z$ m    real a(n,n),b(n)
' R8 j9 l8 R) D# j    parameter(nmax=50)
    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
    real big,dum,pivinv  
    do j=1,n
$ G  d, L9 R* ?1 f6 M4 @4 f       ipiv(j)=0
    enddo
, V# s: x: Q$ p9 j    do i=1,n
( V1 M7 {( |# k" x# F5 m" Y; x7 R       big=0.
       do j=1,n
* p1 l' c1 U1 h; R& q) ?7 Z. I       if(ipiv(j)/=1)then
3 r6 e2 f9 O" Y% q! t          do k=1,n
9 L' e  ^: [& x& s; }9 p# N& n          if(ipiv(k)==0)then
$ k3 K- v' m8 @+ t( y  ~% k          if(abs(a(j,k))&gt;=big)then
+ n" f1 d5 m" p           big=abs(a(j,k))
, ^: y! X8 y/ t  L# F3 Z3 a           irow=j
  m# z6 s# z/ ^6 j, ~           icol=k
       endif
' l' u- r" L% V; S       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
       endif
       enddo
    endif
    enddo
    ipiv(icol)=ipiv(icol)+1
/ b' C' i  E& }    if(irow/=icol)then
: A' e% t6 Q% A# B$ S       do l=1,n
% V3 U- {' n$ o- n, u, `          dum=a(irow,l)
       a(irow,l)=a(icol,l)
       a(icol,l)=dum
) X* {/ q- P% H  _) X       enddo
       dum=b(irow)
       b(irow)=b(icol)
2 X  r9 A; Y  b& m2 {       b(icol)=dum
% l( E( l$ p9 m/ |4 _4 |    endif
    indxr(i)=irow
% M: f4 U4 W5 G9 H1 G6 A    indxc(i)=icol
5 v! W0 x- H# @    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
& l5 E8 p( Q; @5 H5 x    a(icol,icol)=1.
    do l=1,n
9 e- B+ G3 x# L6 J, S        a(icol,l)=a(icol,l)*pivinv
  D9 F9 i2 Z6 f: i' x    enddo
    b(icol)=b(icol)*pivinv
    do ll=1,n
! v9 \! f+ k9 n$ [+ b! a       if(ll/=icol)then
8 V' S+ y* z  @, |' `- b1 h% q, r          dum=a(ll,icol)
: c9 g% F2 M' T- h# L       a(ll,icol)=0
5 z5 s4 Z4 g' p7 E" Y; H       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
2 {% f5 j# e! g& A       b(ll)=b(ll)-b(icol)*dum
       endif
% h' m$ U2 o/ C4 C    enddo
: ?+ }6 h# b& x3 l4 V# r" t    enddo
8 i, y! x, `0 ]& u# R6 e) y. I- }+ n    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
       do k=1,n
          dum=a(k,indxr(l))
- K. r0 b, P. v, A       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
       enddo
2 z" F! f9 p! ^0 q* ^: q) |8 P    endif
' B  W. m6 q& v. e    enddo
    end subroutine gaussj
& f8 ~) V7 A8 h101 end
! l  F1 e. D: T</P>
8 F  X0 ]/ y! Q7 n<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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