BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
1 U( f! x, W9 H- C! m, x    !!!输入函数信息，输出函数的稳定点及迭代次数；
8 {' e" ^  o% j" [& j4 }    !!!iter整型变量，存放迭代次数；
" s' I4 t5 e+ _! W    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
0 _* v) ?! A/ F2 M  U+ Y3 B: [    !!!dir实型变量，存放搜索方向；
6 ^3 `6 b: L( h) t, C/ p    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
: y1 r7 `4 v' V% @7 {' U$ J4 X    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
' I/ k+ I/ ]% D6 s% F1 k    real::x0,tol
    integer::n ,iter,i,j
    print*,'请输入变量的维数'
    read*,n
    allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),gradt1(n),x1(n))
% a7 c! t5 k$ B* w    allocate(hessin(n,n),B1(n,n),G(n,n),G1(n,n))
    print*,'请输入初始向量x'
    read*,x
! S8 w& d9 j& E4 {# s3 E; g( p    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
. c  \) z  B5 J3 L$ x: j2 u( X( o    read*,b
    iter=0
. C+ e) a0 B' m8 b& K  j7 V tol=0.00001</P>
<> do i=1,n
8 e, b+ b& i! T$ X6 Y    do j=1,n
       if (i==j)then
       B1(i,j)=1
9 t9 q! j1 i1 t( O3 a    else
' G% P6 [0 K, {, `4 d       B1(i,j)=0
3 A7 U" |0 S* o4 h) b    endif
    enddo
enddo   
5 g$ R8 W7 t* A; s; \* N    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
* h4 C8 k- z$ a* T2 ?  G4 l        !print*,'极小值点为：',x
# ?- C; E4 P- R. K" [: `) B     !print*,'迭代次数:',iter
     goto 101
    endif
7 U1 S$ T& D3 R* w/ M' ]. v% V. K call gaussj(B1,n,(-1)*gradt)
dir=gradt
8 I' O* g% l9 G1 ?  }8 |    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir
$ z4 Q/ w/ I8 g) s) I$ H& J3 c gradt1=matmul(hessin,x1)+b
3 c/ v) R4 W& e9 m# k% U4 h s=x1-x
y=gradt1-gradt
call vectorm(gradt,G)
G1=G
6 B; I! S# t  ]% k$ r2 d0 @3 A  _1 p call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
' z3 Y' F4 ^! J, u7 f3 ` gradt=gradt1
- j- p; ~6 m) d; K    iter=iter+1
  if(iter&gt;10*n)then
    print*,"out"
    goto 101
) m) J! d# q% B3 [- K( N endif
    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
    goto 100
7 \* i* j9 `& f1 K5 }4 J, }    contains</P>
<>    !!!子程序，返回函数值   
1 h1 e, j# B$ r    function f(x,A,b) result(f_result)
& r2 p2 c) g8 U% K% s    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
; v2 n. S8 o- Z: c+ q    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
1 w6 G7 g! ~/ f# C    end function f
!!!子程序，矩阵与向量相乘
; A* B8 ?5 g9 ^! W+ t- h0 l# C subroutine vectorm(p,G)
- H- S4 \7 d' o5 r6 _+ h+ | real,dimension(,intent(in)::p
real,dimension(:,,intent(out)::G
2 N7 [4 `; @/ j& F' r) A1 E6 j  V/ T n=size(p)
do i=1,n
    !do j=1,n
: G0 h1 u. }0 {$ d, W* s( K       G(i,=p(i)*p
! T+ V% s8 R* G1 F* ?9 \9 [    !enddo
0 t3 n/ G) N/ ^ enddo
# G% Z1 |2 b3 y: x; i end subroutine
$ H. A8 I) P, g3 ~2 I( \
, ]9 A1 n# m3 H2 Z8 }    !!!精确线搜索0.618法子程序 ，返回步长；
5 A3 s+ {6 U/ T1 p    function golden(x,d,A,b) result(golden_n)
    real::golden_n
    real::x0
    real,dimension(,intent(in)::x,d
& p7 V" n$ P0 c) `7 e; l    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
    dx=0.1
    x0=1
    x1=x0+dx
! w7 V+ f* u7 Q; s2 C5 R    f0=f(x+x0*d,A,b)
    f1=f(x+x1*d,A,b)
/ E( |+ c: V, Q2 F2 z    if(f0&lt;f1)then
" ^1 z0 T3 I( K% j  H4       dx=dx+dx
        x2=x0-dx
1 a- o  T; k- B* C# E. _        f2=f(x+x2*d,A,b)
/ p9 e! ~1 ^$ k; O7 Q        if(f2&lt;f0)then
' h$ c$ Z8 X7 C: k' D           x1=x0
1 L$ ~% H. l+ G( ^        x0=x2
        f1=f0
        f0=f2
        goto 4
        else
- u+ c. h4 F8 l4 m           a1=x2
8 I  o1 e9 h- z$ m! b. c; i& s        b1=x1
        endif
    else
. M7 w9 `9 y! ?" [' u$ W2       dx=dx+dx
8 O" x: I+ p3 {! s; Q% r# |9 _        x2=x1+dx
        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
           b1=x2
3 Y2 w% n  r6 {* G5 Y- ^        a1=x0
) l" j% v) y5 g$ p. D+ ?& i- _        else
           x0=x1
' v  u( L% q! v5 k4 S        x1=x2
        f0=f1
* `" s: R: b( k& ^% I        f1=f2
" l% N, O9 L; w! r/ w  f% t        goto 2
) \! ]* V3 k4 [$ {" R) R$ U$ P        endif
5 g: t, y. V" j& X    endif
! Z1 I: H1 \) I+ M- _  ?9 {4 o    x1=a1+(1-r)*(b1-a1)
4 W7 ~9 h9 ]1 ^  F8 e    x2=a1+r*(b1-a1)
& y; [5 H$ p- H" r4 F$ D    f1=f(x+x1*d,A,b)
    f2=f(x+x2*d,A,b)
3   if(abs(b1-a1)&lt;=tol)then
7 }, d$ ?, H1 A/ Y9 g- w. O        x0=(a1+b1)/2
    else
* e6 ]; n- U+ [4 l6 T        if(f1&gt;f2)then
        a1=x1
) K! V% q7 V6 X+ M+ b0 i7 i        x1=x2
        f1=f2
7 R' P. `, s2 D1 _  j        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
9 P1 u4 ~* |1 Q' D* b  u) l        goto 3
: O" A* O  Z# w3 G     else
. K. j9 s0 k( }9 p# n        b1=x2
        x2=x1
        f2=f1
2 I- s3 E. a; u8 q) J        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
        goto 3
     endif
  Y- C$ `4 S5 t1 o. r    endif
    golden_n=x0
; J/ l# c, g: F; x    end  function golden</P>
% @3 `* p4 N. @' S<>
7 q2 ^' ~5 N, K$ M& \    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
    integer n,nmax
! E- S$ S: R' L/ g    real a(n,n),b(n)
    parameter(nmax=50)
9 Q, `( u+ [( p. s3 X( T8 C    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
    real big,dum,pivinv  
    do j=1,n
       ipiv(j)=0
    enddo
+ ?" N! r. {' c' j    do i=1,n
! j( Y6 W2 _' T0 i7 w& W       big=0.
       do j=1,n
1 l: e  g0 ^3 w" L; j       if(ipiv(j)/=1)then
          do k=1,n
: |: O4 I2 N7 e2 F% j4 U" u  \          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
4 ^: k6 r* o9 ~2 z( g! A* e8 X           big=abs(a(j,k))
           irow=j
+ w) J# `* F* }! G3 o           icol=k
9 x! T# W# Y! Z       endif
, U% _  a& V: H5 j- u7 L' c" E       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
. M% a2 m: c1 O* N% k+ y" i$ }       endif
       enddo
    endif
    enddo
    ipiv(icol)=ipiv(icol)+1
) ]* v& n5 {* }" C- H) G1 G( Z  O: V    if(irow/=icol)then
       do l=1,n
          dum=a(irow,l)
/ N2 J" {) c1 W1 W! N/ q  \3 F       a(irow,l)=a(icol,l)
: _  m' p7 @* Y, Z/ y$ E; \* u7 C  o       a(icol,l)=dum
       enddo
       dum=b(irow)
       b(irow)=b(icol)
       b(icol)=dum
    endif
; S$ X' o* A! Q6 n/ P* j. c" Q    indxr(i)=irow
$ c6 c; X) p3 b# X$ y3 A6 D    indxc(i)=icol
    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
! E$ `9 Z0 ]- }2 [0 W7 k$ \7 S, r    a(icol,icol)=1.
% G* x0 C, M  J  O2 M4 T5 s    do l=1,n
) e9 w- z/ L9 M# x7 [! i1 R! {7 ~        a(icol,l)=a(icol,l)*pivinv
    enddo
5 ]) R" ^$ ?/ {5 C; m    b(icol)=b(icol)*pivinv
: e9 S$ e& b3 F' ~% y    do ll=1,n
       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
; ~  |3 Q" Y3 @" s) Q& e       enddo
' N; B; M1 Y, K8 P$ D) Z) w+ ^& z9 |       b(ll)=b(ll)-b(icol)*dum
2 V' Y: e4 E# a8 j       endif
    enddo
# X) Q6 X: D0 O/ t) R    enddo
    do l=n,1,-1
+ ]$ y( s) g3 w: z       if(indxr(l)/=indxc(l))then
1 Z# S$ L5 c) y: t  b/ }+ a/ c9 u       do k=1,n
+ C6 s# T8 \- S! R% V0 b          dum=a(k,indxr(l))
3 J7 B, e+ j( K6 \' w       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
       enddo
    endif
" x8 R$ {1 Q& \5 f9 u( j+ E  V    enddo
    end subroutine gaussj
101 end
( V. I3 s+ ?; g9 I0 o: R0 D</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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