BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
" ?- f# u" \' l2 F/ Z4 {    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
: M% _" d* b+ d' c9 u6 D) \8 u! A    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
# ?* B! t: ]3 {/ R    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    real::x0,tol
    integer::n ,iter,i,j
" k3 R, d7 i. _  z: @4 T" l8 L    print*,'请输入变量的维数'
/ r7 s1 c/ K( K+ s1 S( `/ {" C    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))
" z* x6 q8 N: o    print*,'请输入初始向量x'
1 [: e- }2 p% X    read*,x
$ b" f' _1 Z  i. |    print*,'请输入hessin矩阵'
; x8 }5 j  a6 D: l6 L  V* _    read*,hessin
    print*,'请输入矩阵b'
2 P0 |$ z4 q$ [" g    read*,b
    iter=0
tol=0.00001</P>
<> do i=1,n
5 o" p7 ^2 [* y    do j=1,n
       if (i==j)then
       B1(i,j)=1
    else
       B1(i,j)=0
  |& B* H+ X5 p4 `  M7 A' R4 U    endif
0 k6 w+ q6 i$ H& {    enddo
enddo   
    gradt=matmul(hessin,x)+b
! M' M8 S% U5 j3 b0 S3 D" Y100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
, ^" m; {  \. B6 N5 Z        !print*,'极小值点为：',x
8 T5 R9 Z" p) R3 _     !print*,'迭代次数:',iter
4 Y# I' U. U# U  e     goto 101
& a# I' _( l' [' j- m# k    endif
. V6 A+ n- {2 b/ f. y5 N: v call gaussj(B1,n,(-1)*gradt)
1 X3 E) I  R3 ?4 S) n( B- r* r dir=gradt
    x0=golden(x,dir,hessin,b)
" g0 ^# I9 k9 f! i' [    x1=x+x0*dir
+ g6 q5 B) T9 n9 L. I gradt1=matmul(hessin,x1)+b
s=x1-x
y=gradt1-gradt
. |' i$ P# e, K call vectorm(gradt,G)
G1=G
: X2 V* x; ~" t4 q call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
5 I) o; X7 B& c9 Q7 M( p1 a x=x1
gradt=gradt1
    iter=iter+1
# h4 ]7 I- I  [' ]  if(iter&gt;10*n)then
    print*,"out"
    goto 101
endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
    goto 100
* q& }0 a& {* A1 y- \+ Q+ D    contains</P>
: g) E( q- ?) E0 w4 C<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
, D. L3 v: ?. y$ W    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
$ g$ ?3 O" J' G# ], {) F! e    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
    end function f
!!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
real,dimension(,intent(in)::p
; X* {) t- p) u) P real,dimension(:,,intent(out)::G
  v  C$ ]9 b- H7 ~: K n=size(p)
# ]+ P8 T1 V0 Z4 C  n) D do i=1,n
    !do j=1,n
. S: ^: K* r, h# t2 |0 r       G(i,=p(i)*p
1 [( P: h" W" w% n0 i- N* v    !enddo
: l& ]& w0 |5 Z) h$ p) f# H enddo
end subroutine
. v: r6 C0 s+ m% U
6 V9 a' r  i0 H8 \/ l. w    !!!精确线搜索0.618法子程序 ，返回步长；
0 }2 ~  R6 w( G2 B9 g. L    function golden(x,d,A,b) result(golden_n)
    real::golden_n
    real::x0
    real,dimension(,intent(in)::x,d
7 h4 I* v9 C4 a5 \, P* E    real,dimension(,intent(in)::b
    real,dimension(:,,intent(in)::A
8 p7 y* q1 G  k1 o/ C9 l    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
    parameter(r=0.618)
    tol=0.0001
* C! [: i" u2 i    dx=0.1
    x0=1
    x1=x0+dx
  d% b5 g6 f9 D; T    f0=f(x+x0*d,A,b)
0 Y' o7 \/ S2 N, ^- S+ |    f1=f(x+x1*d,A,b)
    if(f0&lt;f1)then
4       dx=dx+dx
# ?2 K# ~4 j4 X8 O7 {* s        x2=x0-dx
        f2=f(x+x2*d,A,b)
        if(f2&lt;f0)then
: Z2 B3 i  t* t- ]( A1 u           x1=x0
3 |9 _& T! E) t' W% C        x0=x2
, ?+ R7 n2 o. z% e        f1=f0
1 ^6 X: z6 n3 x9 y5 y. W( L3 I        f0=f2
! I8 y3 x+ S2 X* B6 C        goto 4
- P8 ?) [6 a; ?8 F. J# ^0 S        else
           a1=x2
7 |# i& P0 q+ \! q& ?, F6 a        b1=x1
        endif
! i! Y+ u% Y/ \8 d# t6 q& ^" D) s  A    else
2       dx=dx+dx
        x2=x1+dx
        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
6 I8 i5 a% o' E; G8 C           b1=x2
! [9 N& s3 k* u. E3 ~3 f. x        a1=x0
        else
           x0=x1
4 t4 p, E5 |( R        x1=x2
5 Q" l6 D1 L6 U        f0=f1
        f1=f2
4 S# r3 H4 d* x; ^% N8 u: \3 F        goto 2
        endif
    endif
+ ~; `# [# k- }& I: L) j  S    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
    f1=f(x+x1*d,A,b)
    f2=f(x+x2*d,A,b)
- J# t# t: R, D" A, l& o3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
    else
; m$ b6 O$ G+ o' y7 a$ _/ u        if(f1&gt;f2)then
% y9 @% d: j6 j2 |* H        a1=x1
        x1=x2
6 ^4 R1 n! V9 O$ ]2 h( h, v        f1=f2
        x2=a1+r*(b1-a1)
5 H3 }  g( R5 X1 T; V/ I        f2=f(x+x2*d,A,b)
        goto 3
     else
3 c/ a2 D0 R: H: z: T  ]$ Y( j        b1=x2
        x2=x1
1 T5 ?5 j- G3 g! z4 l+ M: v" |  E% j. [        f2=f1
        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
- T7 C9 |( k; V- {0 ?" _% T        goto 3
1 W" Q$ X1 {7 F% _4 C/ y     endif
8 T/ P/ m5 f2 Q+ l! W    endif
3 _% k; e8 g5 s5 H% e: J    golden_n=x0
# x3 g% s* G% H  n4 C    end  function golden</P>
) Q. m* D3 ^) I<>
. [0 c' I+ ~9 g) h8 X* F$ ^/ r    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
7 y2 O8 O$ h5 `" W5 ?' @- Y    subroutine gaussj(a,n,b)
    integer n,nmax
( o+ m5 Y0 p/ y6 x2 r9 V! }8 @    real a(n,n),b(n)
    parameter(nmax=50)
    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
" P- ~( b% c( u) \# w    real big,dum,pivinv  
    do j=1,n
8 A; ~% i9 M+ z4 Z) [       ipiv(j)=0
    enddo
/ W9 H- V) b" h* w( N2 W7 e    do i=1,n
- Q7 k# N$ ?) X; B; o) A       big=0.
       do j=1,n
       if(ipiv(j)/=1)then
2 y' H0 J! E6 S: g0 q1 O  ]+ y9 f6 v2 T          do k=1,n
          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
" A4 E0 I$ T$ s* @' N2 X+ m+ U* j           irow=j
9 |- F1 _" V, L           icol=k
2 @  z- h$ k6 x/ o  `# _2 E$ l       endif
& |+ u* _; k" x& u+ L       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
  ?& j0 @" `9 Y# E/ s  J$ X       endif
       enddo
    endif
+ w0 g/ t- C1 F) A    enddo
% S7 E, V& P# v- X& B    ipiv(icol)=ipiv(icol)+1
8 @! \- v! O; D; R8 A    if(irow/=icol)then
       do l=1,n
& H. T- X0 Z5 c* C          dum=a(irow,l)
' Q2 h4 i' G6 f: ?+ u       a(irow,l)=a(icol,l)
       a(icol,l)=dum
       enddo
0 j, {( n7 \5 p4 ^       dum=b(irow)
1 Q% Q: F2 J7 N9 d9 ]2 q       b(irow)=b(icol)
       b(icol)=dum
    endif
" K6 l5 C, ^! X2 n4 d" G    indxr(i)=irow
2 z: D3 t! m9 E* U    indxc(i)=icol
    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
; t$ C8 _, v+ i! e    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
    enddo
    b(icol)=b(icol)*pivinv
    do ll=1,n
       if(ll/=icol)then
          dum=a(ll,icol)
       a(ll,icol)=0
, K$ K- T* r5 \. w5 N& I( J       do l=1,n
" \. Q: i7 v* s4 D" m( o          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
       b(ll)=b(ll)-b(icol)*dum
  O) z  u$ L( {  C# m- U/ ^" ?2 }8 O       endif
' A( y" v' W: h) E- x4 y% u7 T    enddo
% W1 v6 ]3 ~  R- g6 m    enddo
    do l=n,1,-1
       if(indxr(l)/=indxc(l))then
       do k=1,n
          dum=a(k,indxr(l))
7 ]6 E- J8 x( _  ^# j( `8 s       a(k,indxr(l))=a(k,indxc(l))
3 ^) P' y) g/ H) l6 k  F; G       a(k,indxc(l))=dum
( f8 f% C; l  d" h, b$ H4 u2 N& }       enddo
    endif
    enddo
! Z! p. J% I2 Z0 q) `, K* s    end subroutine gaussj
/ X* i) J  ^9 S  K4 z6 y% k101 end
1 M; j9 {7 V2 z" O7 x: o" i</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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