BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
$ @6 w2 R8 k4 P6 O" p    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
' ~* B& I( A) P+ y' `5 m    !!!dir实型变量，存放搜索方向；
    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
7 D, ]) S5 R/ R! P9 \# T% K2 h    real::x0,tol
. D9 a" A2 W: X6 W; d& p    integer::n ,iter,i,j
    print*,'请输入变量的维数'
; F; q7 H" f5 Q    read*,n
; H/ _/ _' k) @" K    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'
    read*,x
    print*,'请输入hessin矩阵'
/ m0 q: b0 g: `5 u' L    read*,hessin
    print*,'请输入矩阵b'
2 w; w/ p- x& S! I/ Z    read*,b
0 u; u* ?( I) ^3 g* r    iter=0
: e" W, b8 G+ A tol=0.00001</P>
<> do i=1,n
    do j=1,n
       if (i==j)then
       B1(i,j)=1
' n  A- ~) C& u    else
+ e; P+ I5 w+ I3 f) r( [. h       B1(i,j)=0
    endif
    enddo
& n) n" k  g; H" M8 _  F5 O enddo   
    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
9 ?; b. ]. z7 g3 H6 k, r9 V4 X        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
$ M% L$ L) P0 V/ z3 v9 Y+ k  L     goto 101
+ G) e* S" w2 ~# [1 U: Z    endif
2 K6 w; c# x# v0 u, i- \ call gaussj(B1,n,(-1)*gradt)
dir=gradt
9 J8 l# u% N  \& [  M) t4 w    x0=golden(x,dir,hessin,b)
' q, ]) z7 `/ j  I2 |, R$ c5 m    x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
1 P) u4 E. M" C  ?0 e! R s=x1-x
5 @9 U# h6 T& m0 _5 S3 U$ I y=gradt1-gradt
; z8 [' O" M: Q0 T3 @4 Z call vectorm(gradt,G)
4 B/ Q! }* r9 _" `$ y2 M# L G1=G
2 S% @* t  e/ Q# D# P9 y; H& ~ call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
gradt=gradt1
& o+ B3 E& i; j/ S3 a$ [4 b    iter=iter+1
  if(iter&gt;10*n)then
) Y1 G+ L* ~9 I0 e    print*,"out"
    goto 101
, e; ~9 l. D- r, A) e! }4 B5 n endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
0 w! \! _* d  r. C    goto 100
6 v7 [# ~. I( }' M2 h+ t& K( t( {    contains</P>
<>    !!!子程序，返回函数值   
- D1 i3 n4 K0 `. Q$ B    function f(x,A,b) result(f_result)
* l2 N' i/ t( n" M+ V4 @0 X/ z3 d    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
5 Z5 m1 F. p$ P    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
' \6 F- r" K! \" P    end function f
!!!子程序，矩阵与向量相乘
8 Z# t3 E/ p3 g" @1 H& Z& Y% s subroutine vectorm(p,G)
7 W( F1 f0 m& q. @7 a& y3 v real,dimension(,intent(in)::p
real,dimension(:,,intent(out)::G
. F6 P$ {$ S+ y, k n=size(p)
do i=1,n
    !do j=1,n
       G(i,=p(i)*p
    !enddo
enddo
" e* W, C, p8 [% O( Q end subroutine

; R& k6 I! W' ?% [    !!!精确线搜索0.618法子程序 ，返回步长；
    function golden(x,d,A,b) result(golden_n)
    real::golden_n
' v4 ^4 e/ u+ k7 \: ?  j    real::x0
" ]) g) ]. v) r: `, E9 o+ j8 W    real,dimension(,intent(in)::x,d
) A, w; v2 E3 Z3 v! [% Y4 T$ j/ ~6 ]    real,dimension(,intent(in)::b
    real,dimension(:,,intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
# Q3 B  B# g' N6 J4 `    parameter(r=0.618)
    tol=0.0001
    dx=0.1
    x0=1
  H& o2 E$ e# k# R    x1=x0+dx
) n, O- t+ m0 U5 J& u% h8 N9 y    f0=f(x+x0*d,A,b)
' |: O0 K) o1 X: K0 [    f1=f(x+x1*d,A,b)
$ a, i/ F0 V+ P% c7 V    if(f0&lt;f1)then
4       dx=dx+dx
: X0 d+ ?" ~0 H% `" ^0 V0 i2 r        x2=x0-dx
. v( z# X- @' w        f2=f(x+x2*d,A,b)
- g4 P: c$ l9 a# m        if(f2&lt;f0)then
           x1=x0
7 O0 |- i  D. p# ~        x0=x2
! n( |3 ?1 \; j; g1 j! u        f1=f0
: A' Y6 u# }3 e        f0=f2
        goto 4
        else
# p2 e5 g+ n8 |7 m           a1=x2
9 c5 z0 P1 Q9 b1 p        b1=x1
        endif
    else
0 d  ?8 K/ N! ^  x2       dx=dx+dx
' x. l" Z8 d' J7 v        x2=x1+dx
; c' B7 a7 n, s5 [" d8 |! G        f2=f(x+x2*d,A,b)
4 `& w9 S1 M2 ]/ [! X/ q9 x- i  ?        if(f2&gt;=f1)then
           b1=x2
; ^' v+ O* e& q; c- ]        a1=x0
4 Z: N& G. D4 w& B# n/ I( L. b; ?        else
           x0=x1
        x1=x2
        f0=f1
' Q5 Q/ N, k+ f) A/ u/ y        f1=f2
        goto 2
        endif
    endif
    x1=a1+(1-r)*(b1-a1)
5 m- a5 w/ y3 ]! Z0 D% ]# U* X$ n    x2=a1+r*(b1-a1)
% c: n3 K% O) E+ {5 @    f1=f(x+x1*d,A,b)
    f2=f(x+x2*d,A,b)
' N- ^1 e( [/ M& S" u3   if(abs(b1-a1)&lt;=tol)then
- M% v8 k. K% C- Z' p7 l6 q        x0=(a1+b1)/2
    else
        if(f1&gt;f2)then
0 W/ w" ?4 V# A* w2 E* {        a1=x1
. U& n* Y; Z2 @        x1=x2
        f1=f2
        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
  M1 ~& Y& y0 n4 z6 }1 o        goto 3
     else
        b1=x2
        x2=x1
* C+ y* L( k: J$ ?        f2=f1
* ~) V3 i8 }0 Y3 n- ~( Q& m: O        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
  y- t5 W$ J- ^  d        goto 3
0 g1 l& v& d2 Y     endif
$ s" {# f% M! O! G9 S    endif
( R' d% f; d7 ]    golden_n=x0
    end  function golden</P>
' F/ D, {: t* a<>
4 X) u7 H4 T% l% @    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
    integer n,nmax
- x. l% I" t7 n* w    real a(n,n),b(n)
    parameter(nmax=50)
    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
. R6 B* `. J# D. W3 [    real big,dum,pivinv  
/ {/ {- H2 G, F; \+ J* N    do j=1,n
* V3 s! Q- Z* S" Y2 s& ]$ d+ ^       ipiv(j)=0
4 }$ f; u& m0 @$ E: ]3 \- M( s    enddo
    do i=1,n
5 ~5 Q0 n  ]9 P  d- {3 w       big=0.
$ l  g3 a3 V3 S5 e' E( y       do j=1,n
) W& f7 R- z  ?  k8 X* a       if(ipiv(j)/=1)then
          do k=1,n
* B# l! R- P" E: c8 Y          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
/ Z8 ^, a9 n8 g, v, f           irow=j
           icol=k
       endif
       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
2 \' Z$ s: W0 H2 s       endif
1 \9 t, e# o7 J% v- T$ }2 a       enddo
1 |. f8 R' |; M; D' C    endif
    enddo
    ipiv(icol)=ipiv(icol)+1
& N( }4 Q/ E+ }4 e+ s; J! P    if(irow/=icol)then
       do l=1,n
& s% h+ F8 @* V3 u/ E          dum=a(irow,l)
$ C# ?1 @6 _8 T7 l4 K       a(irow,l)=a(icol,l)
/ N2 E! m) x. z3 D       a(icol,l)=dum
       enddo
( z: _$ J  A+ Z       dum=b(irow)
       b(irow)=b(icol)
5 a, O3 p7 y) W! H6 \* z       b(icol)=dum
8 W' S  G0 w2 b  z    endif
1 n9 u3 B" {& e5 Y    indxr(i)=irow
    indxc(i)=icol
8 h& `$ i9 ^* c2 z3 X+ @    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
& [; ]# \% C1 z3 G. Q    do l=1,n
4 x% a/ D6 T) K) U        a(icol,l)=a(icol,l)*pivinv
" O$ E3 f: r) I+ O' N+ [    enddo
* e" a2 b6 ?4 X9 a9 q3 `8 @/ d/ |0 w    b(icol)=b(icol)*pivinv
. C' B6 x$ d" r% o- l6 G    do ll=1,n
7 c% `& J4 |. W8 L& r2 E       if(ll/=icol)then
          dum=a(ll,icol)
       a(ll,icol)=0
       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
0 }+ R: K+ M4 M$ B' A       b(ll)=b(ll)-b(icol)*dum
       endif
( h3 ^2 e. a! g% ?" h    enddo
  P& g( B$ W4 z, _6 o" E    enddo
1 A. H; A. D. w' U. B8 n& ]' d    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
4 T* h9 Q# P) _- }6 |       do k=1,n
" ^, e0 P% o1 {$ ^4 K/ c; U          dum=a(k,indxr(l))
       a(k,indxr(l))=a(k,indxc(l))
; q" d8 d/ f; Y# Z- \3 T( E. z       a(k,indxc(l))=dum
       enddo
    endif
    enddo
    end subroutine gaussj
101 end
</P>
2 v$ g8 R! U: f) ?<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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