BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
- U; X* K/ E) b, c! @9 d    !!!iter整型变量，存放迭代次数；
3 a( H3 `% u- ?0 [# y. q! i8 b    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
    program main
' i, [" o! N( H# d" x* M    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
7 P/ y# v3 U5 }. g& L# p9 Q    integer::n ,iter,i,j
. u1 l' U: s. P; M1 U    print*,'请输入变量的维数'
  u( v9 y" m" M: m4 E    read*,n
    allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),gradt1(n),x1(n))
3 J, s9 g' S& w* q    allocate(hessin(n,n),B1(n,n),G(n,n),G1(n,n))
    print*,'请输入初始向量x'
    read*,x
    print*,'请输入hessin矩阵'
" _8 d$ y2 S& I# H    read*,hessin
# [% K5 [: ]( }6 z    print*,'请输入矩阵b'
5 B1 E- n# i5 _; j6 Y6 S    read*,b
    iter=0
tol=0.00001</P>
5 |" e+ g$ ^& k2 r) q<> do i=1,n
- O, V6 D# e$ i% [    do j=1,n
       if (i==j)then
: R7 F. ~. Y: G9 X* u       B1(i,j)=1
9 n! b- W8 \, _! T7 ~9 d    else
       B1(i,j)=0
    endif
    enddo
enddo   
$ l  b! q5 v- r. P* F% m    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
: P7 \7 r7 G8 W; ~        !print*,'极小值点为：',x
' a3 d# H/ b- [; c% z     !print*,'迭代次数:',iter
; }. a, Y! z3 S' r; r     goto 101
    endif
, K: j+ B# P9 i+ P1 h) h2 r4 b call gaussj(B1,n,(-1)*gradt)
dir=gradt
    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
1 `4 u! L% I  l/ a6 Y s=x1-x
$ m( R  P( R" ]( k! t$ ]3 V* { y=gradt1-gradt
4 x: `2 i: r) f3 f9 | call vectorm(gradt,G)
) G) q9 E: W+ Q3 u# v! i G1=G
! I$ s. x; O% x* L9 _ call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
7 C* x+ S" D! ^3 Q# [ x=x1
0 ~/ f/ y+ U, Y: O3 A3 v; _( w gradt=gradt1
) L% ]; Y' d; d6 Q( Z( G+ O    iter=iter+1
. H# l9 z- G/ j4 B3 l% i  if(iter&gt;10*n)then
    print*,"out"
    goto 101
( l( E- @+ S9 f: p endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
7 u6 K* i$ n# L6 O) B; w    goto 100
( j: Z2 X) h& }/ p& o8 V    contains</P>
9 v, x8 H2 K+ u' K<>    !!!子程序，返回函数值   
. D8 {+ Z2 k1 |) Z* @, [* m    function f(x,A,b) result(f_result)
; r0 e/ \  Y- V" \    real,dimension(,intent(in)::x,b
; ]: g% X1 }6 p7 |' ]7 e    real,dimension(:,,intent(in)::A
    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
" V2 A2 B: i* H% A0 U0 w- e+ o    end function f
: o  x- }! q: J/ m" C# ~9 B !!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
real,dimension(,intent(in)::p
: d- B/ h" E& `4 w! m, N: X( T0 h0 W real,dimension(:,,intent(out)::G
n=size(p)
do i=1,n
* B6 K3 w) f5 [8 m    !do j=1,n
2 j* Y( N, F! B1 s6 S1 t       G(i,=p(i)*p
    !enddo
enddo
* G' {) }7 z9 }% q6 ]: x end subroutine
/ d3 H8 j: q# S6 r, j
    !!!精确线搜索0.618法子程序 ，返回步长；
2 X8 k0 [3 J8 p: z    function golden(x,d,A,b) result(golden_n)
    real::golden_n
# m1 {" o6 M. E7 h    real::x0
5 H$ Q9 }0 }/ o4 U: Y    real,dimension(,intent(in)::x,d
    real,dimension(,intent(in)::b
    real,dimension(:,,intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
    parameter(r=0.618)
    tol=0.0001
, b! N6 T9 ]3 Q% F: K& ^5 I    dx=0.1
+ U. V2 f$ G1 y! `% y, T    x0=1
% T" r; V- Q/ O3 b. l' n* F    x1=x0+dx
    f0=f(x+x0*d,A,b)
    f1=f(x+x1*d,A,b)
$ e% J5 |# R( E' ^) J    if(f0&lt;f1)then
4       dx=dx+dx
5 l% K  d/ u" W8 D+ l        x2=x0-dx
+ t" G- [6 @, W        f2=f(x+x2*d,A,b)
6 f) K2 M* Y2 z( b4 r- j        if(f2&lt;f0)then
; {* a: d4 n- I! s( u           x1=x0
6 h! ]* G/ c0 }) c+ i9 D7 Y        x0=x2
5 s0 U" u; L6 C% ^8 H8 e( o        f1=f0
' n$ g+ T5 f/ r6 M7 G, G) u        f0=f2
        goto 4
/ ]: |4 ~2 ?8 a  i1 P        else
           a1=x2
+ ~$ A, ]8 n! N/ T        b1=x1
        endif
/ A8 W: E0 ^% T$ ]    else
2       dx=dx+dx
        x2=x1+dx
4 y3 ?9 J2 q! T3 w0 @        f2=f(x+x2*d,A,b)
/ p  h% m- T( m5 @2 Z        if(f2&gt;=f1)then
           b1=x2
% d+ [( j  C3 W) f% s5 p        a1=x0
; n4 Q% T+ a' S0 ^" U8 j        else
           x0=x1
        x1=x2
% _0 o6 a- t' z0 P% @9 e. `" J        f0=f1
        f1=f2
        goto 2
! e/ u& I) q: j/ l0 \+ N; X        endif
  ~# K. `% ?3 l5 o, h& W    endif
0 J3 O0 H' i( [* i) h9 O    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
' P" h) |1 F% r5 e: S    f1=f(x+x1*d,A,b)
    f2=f(x+x2*d,A,b)
2 }" F) A  O. n( S5 E3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
, B! x2 ?* k+ {' S8 S    else
        if(f1&gt;f2)then
! Q6 v8 B# [' `9 A) [- }/ N        a1=x1
        x1=x2
( Y/ z& |) a! n" R        f1=f2
( D6 C8 [; h; z, r; A        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
6 T1 f# J0 c! }7 ~, z& T        goto 3
     else
' ~( D$ |, m4 K- ^) T: s        b1=x2
        x2=x1
! P8 M! e( d9 n* {& r1 t        f2=f1
        x1=a1+(1-r)*(b1-a1)
$ N2 a! T+ @. r" |6 L8 z( R: D        f1=f(x+x1*d,A,b)
$ o/ i* Z) y" C! Z2 p( n        goto 3
     endif
8 ~# A/ z5 t5 ~5 s) E! h    endif
    golden_n=x0
4 Y, ?' z6 B; X6 T    end  function golden</P>
<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
    integer n,nmax
    real a(n,n),b(n)
    parameter(nmax=50)
2 {4 y0 B; |  F$ e& f: c5 `8 P    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
    real big,dum,pivinv  
5 P3 o! {. L" H, T' |/ P" p    do j=1,n
       ipiv(j)=0
: X: \4 }0 b( X, S3 r: Y' b2 k$ B    enddo
    do i=1,n
       big=0.
3 M1 {" |9 o2 O+ X+ `/ r       do j=1,n
) N6 n- Z6 Q3 V& k) K       if(ipiv(j)/=1)then
9 F+ X$ |% t' q3 G* X          do k=1,n
$ o- w* w5 `9 ?4 `          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
           irow=j
7 ~5 y5 A3 [) @) s           icol=k
; o+ y) c$ {  I- `5 \       endif
       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
7 R0 X) \5 q7 T$ Z" `       endif
; ], p* n/ s) i% o$ e! W       enddo
    endif
    enddo
    ipiv(icol)=ipiv(icol)+1
0 G$ k9 @8 v0 ~8 }& S    if(irow/=icol)then
       do l=1,n
          dum=a(irow,l)
       a(irow,l)=a(icol,l)
       a(icol,l)=dum
# t5 ?0 k" E* @       enddo
- j* V$ {; Z% A; v- e2 t- D       dum=b(irow)
       b(irow)=b(icol)
3 a# Y+ z* O& O) R. Y4 e& D; H       b(icol)=dum
& ?- x1 E  E2 {; y2 n! K' {2 h    endif
/ p4 a1 q) F& M1 p    indxr(i)=irow
    indxc(i)=icol
, p  U* Z% C0 ?5 f" u    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
* m9 r* [0 s$ \/ P6 s1 r" K5 o4 ~8 k    do l=1,n
; z- _* t, z2 `8 v, }) ^        a(icol,l)=a(icol,l)*pivinv
    enddo
    b(icol)=b(icol)*pivinv
    do ll=1,n
( V2 y# P* f, H! D: \0 U) T6 x( k       if(ll/=icol)then
6 A% y, V2 F$ O0 I/ f  \          dum=a(ll,icol)
       a(ll,icol)=0
1 `8 Y1 c0 P1 e9 x) G       do l=1,n
3 ^$ a/ H2 I& B+ H9 ^          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
       b(ll)=b(ll)-b(icol)*dum
       endif
) ~  U7 Z1 ?$ A1 [    enddo
  _& ?7 b- F9 t    enddo
2 t( \, U1 \. t9 c/ B, `7 @8 h/ T6 Q    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
       do k=1,n
* d( ]9 Q9 u# t" s0 _1 l          dum=a(k,indxr(l))
' t! i) ?, k: T" Y       a(k,indxr(l))=a(k,indxc(l))
: d' t% @, \1 C! t* Q+ Z' k5 E       a(k,indxc(l))=dum
4 U: ?( m& l. e9 [! A0 b       enddo
8 K6 P% W- {0 c& m1 ?( Z    endif
5 I  W  s- {8 R8 P/ D    enddo
    end subroutine gaussj
101 end
</P>
9 ^9 k, [4 n2 d% R% |1 P# _<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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