BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
/ V% n, B  Q  ?$ U    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
4 q7 P0 h3 E9 f, p    !!!dir实型变量，存放搜索方向；
    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
    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))
    allocate(hessin(n,n),B1(n,n),G(n,n),G1(n,n))
2 y7 A- I8 F8 {( ], [    print*,'请输入初始向量x'
9 J8 B  A  K9 Z5 z    read*,x
( t( G7 e' i( [- I/ N9 @    print*,'请输入hessin矩阵'
    read*,hessin
  ~- h$ K! N, [" u; M    print*,'请输入矩阵b'
    read*,b
    iter=0
/ A$ ]  Z6 \! D# w* K5 J+ Z tol=0.00001</P>
<> do i=1,n
7 L: b- K$ c" T8 a( N! s1 T    do j=1,n
       if (i==j)then
       B1(i,j)=1
, ]3 A& G: x+ [8 b    else
2 E1 {9 U. p+ ?# s# c3 U% f# w+ c       B1(i,j)=0
    endif
    enddo
/ E  ~* w5 }2 f enddo   
    gradt=matmul(hessin,x)+b
! i5 ?6 B$ l* E  o$ U9 J( G100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
     goto 101
! J# Y3 _/ ]# x3 e% O) Q    endif
! I% x6 e) U& o: P call gaussj(B1,n,(-1)*gradt)
dir=gradt
    x0=golden(x,dir,hessin,b)
) R) a0 j: A& J; s# J    x1=x+x0*dir
* d. e1 M' h/ ~7 R- z gradt1=matmul(hessin,x1)+b
- o  c- k/ d! c" p s=x1-x
4 Z' [7 g+ \+ P0 I* T y=gradt1-gradt
call vectorm(gradt,G)
' e; w  ~" m# M) k4 o G1=G
% Z& m! P: K( D' Q" U7 S call vectorm(y,G)
- B: x% X1 W: ^ B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
7 Q! D  j% Z* v' g2 i# n x=x1
gradt=gradt1
7 m6 Y3 l4 i/ R0 C% I6 L    iter=iter+1
  if(iter&gt;10*n)then
    print*,"out"
    goto 101
endif
* R6 `6 {9 ~! Y- i& F  d* a+ W0 K    print*,"第",iter,"次运行结果为",x
print*,"方向为",dir  
- r2 |9 p: V' ^0 U( _) |    goto 100
    contains</P>
8 Y6 F( d$ Q" U& M! v9 y<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
/ u, f+ N4 D) W8 f4 A    real,dimension(,intent(in)::x,b
    real,dimension(:,,intent(in)::A
5 [* m1 C* s! j- k/ ]    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
    end function f
2 s4 W8 N. W" J# Z3 a% L7 d, [* Q !!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
real,dimension(,intent(in)::p
9 r& @  S  g1 d" Z9 g# t real,dimension(:,,intent(out)::G
+ i9 n& F5 i$ |" h4 G n=size(p)
do i=1,n
    !do j=1,n
       G(i,=p(i)*p
% }& o3 ~) G" l2 {" s3 f: }    !enddo
3 f8 a3 M6 ~, { enddo
& Q9 T( x% S4 ?' G/ u7 P) @6 Q$ W end subroutine
& p! I. c$ R$ p) o  |1 _+ M4 t) b* R
2 s  A. W5 J$ N( t; x) I7 P" G: O    !!!精确线搜索0.618法子程序 ，返回步长；
$ V& ^& N5 I% d    function golden(x,d,A,b) result(golden_n)
    real::golden_n
7 B2 ^! C- }) x% q% F: C* s    real::x0
    real,dimension(,intent(in)::x,d
8 ]3 @6 g2 f+ ]0 S% e5 L0 ^    real,dimension(,intent(in)::b
1 R/ t5 Q7 P$ j4 b. H( L    real,dimension(:,,intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
1 [. J. O+ h+ L3 H    parameter(r=0.618)
0 b0 i( r8 G: I; G! U    tol=0.0001
! V! M- O; N0 p' t, M" ?% A7 ?    dx=0.1
2 q) V+ i8 {  w5 ~# T) O    x0=1
    x1=x0+dx
% C9 M/ G) D6 K7 T0 L    f0=f(x+x0*d,A,b)
    f1=f(x+x1*d,A,b)
" S  o" T! |) c) w    if(f0&lt;f1)then
4       dx=dx+dx
8 d1 }) b% E8 g% L7 d. \* _        x2=x0-dx
' ~. q, P2 E' @7 |1 N: O        f2=f(x+x2*d,A,b)
; D; x' o- S5 @$ H        if(f2&lt;f0)then
$ T5 ~8 W  i: r% m' y0 r3 P) N           x1=x0
8 I; t& Y, m9 l        x0=x2
! V. v2 h7 C' W) u' y2 L% ~  T        f1=f0
- q7 B$ N9 o: S5 E2 z" r6 b) R        f0=f2
        goto 4
: c. G6 U4 V8 {9 u2 r9 e        else
( U: C- G# v- W( N           a1=x2
        b1=x1
0 x% `9 R! f5 n        endif
# r! M! q- _+ L* L/ y% Q# t    else
2       dx=dx+dx
5 \/ J5 w( s( B4 u: V        x2=x1+dx
) u- H; k$ \" O& Q! C0 l4 H        f2=f(x+x2*d,A,b)
! Q+ u7 I1 Z% @- }* J        if(f2&gt;=f1)then
           b1=x2
# [# m  r: w. [' {  z& X- _        a1=x0
        else
           x0=x1
: f, X, S. w7 D( r0 O  i6 a        x1=x2
0 v2 V$ j3 w4 s3 N4 ]8 H        f0=f1
( v1 f' E$ K+ x* o3 K# E        f1=f2
        goto 2
        endif
    endif
    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
    f1=f(x+x1*d,A,b)
+ u# N2 p1 e4 u    f2=f(x+x2*d,A,b)
: y9 D6 j2 E# J9 v0 j3   if(abs(b1-a1)&lt;=tol)then
5 m4 Q0 F& t9 X4 z) f2 X" S        x0=(a1+b1)/2
    else
4 D& z$ P4 |( z1 F/ Q3 Y5 u        if(f1&gt;f2)then
7 ~7 @9 ?9 g: `7 S, g; B        a1=x1
        x1=x2
( v: U. p( A* J* _- G& _        f1=f2
        x2=a1+r*(b1-a1)
; j) N$ T; L9 t2 \        f2=f(x+x2*d,A,b)
& {9 g7 }7 Y7 y: C& r        goto 3
     else
: B  b) R! q  M2 X- p1 C        b1=x2
        x2=x1
        f2=f1
        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
' b0 l  [  s( j8 M7 i8 ^1 R        goto 3
- H2 ]0 L* R0 U3 p+ {& y. t; q- l6 g     endif
+ \1 y' ~7 j2 ?* m! u& J  x    endif
' a& E: w5 @1 U6 t% B& v    golden_n=x0
+ R3 Z/ X" o9 n) f3 \    end  function golden</P>
. \. B, A9 d) q; ]2 r, I- }<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
    subroutine gaussj(a,n,b)
* w1 s- x# H% K' {% F    integer n,nmax
    real a(n,n),b(n)
) h" h* E  G+ U    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
: t- z# I$ m4 e       ipiv(j)=0
    enddo
    do i=1,n
9 e. X4 P) L/ X0 n& J/ L8 `1 o       big=0.
       do j=1,n
       if(ipiv(j)/=1)then
          do k=1,n
0 `! r8 v5 b8 ^          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
& @) {; v$ F/ B8 Q, `6 g. q           big=abs(a(j,k))
5 U7 M4 ^. I9 f+ m           irow=j
; o; X% Q9 T1 N3 C7 @           icol=k
       endif
       else if(ipiv(k)&gt;1)then
9 N$ u6 o" o1 i$ s7 J          pause'singular matrix in gaussj'
       endif
       enddo
    endif
7 e5 V- o& u! J* v/ n    enddo
    ipiv(icol)=ipiv(icol)+1
- a9 X7 L4 _7 b  G  l    if(irow/=icol)then
       do l=1,n
          dum=a(irow,l)
4 f0 T6 x: q2 k; q/ T       a(irow,l)=a(icol,l)
% b  @0 ^$ i- @8 V, ^  A' w       a(icol,l)=dum
6 b' [7 n& a1 d7 G& X, S9 o       enddo
       dum=b(irow)
       b(irow)=b(icol)
       b(icol)=dum
% A9 B% \5 J9 {- g$ w1 ?$ l    endif
    indxr(i)=irow
  J  W6 s% O$ F( x2 ~1 N    indxc(i)=icol
    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
  L' W* U9 r5 f8 z) `    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
    enddo
* m7 ]5 B+ ]+ g$ x6 A    b(icol)=b(icol)*pivinv
" \6 J& g9 S& f7 H5 ]2 J7 \% P    do ll=1,n
  b$ I: Y  T+ f" n       if(ll/=icol)then
) a6 d/ q9 C' L          dum=a(ll,icol)
       a(ll,icol)=0
       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
1 r4 z' v  I& T       enddo
       b(ll)=b(ll)-b(icol)*dum
" s% ^$ p+ g5 |9 v       endif
    enddo
    enddo
    do l=n,1,-1
" G' ?6 I) \/ ]( E8 E7 R       if(indxr(l)/=indxc(l))then
8 b& h& x3 D$ T% o7 O# R3 i2 \       do k=1,n
6 ~) X4 M5 l8 T          dum=a(k,indxr(l))
' P. t. G+ \4 t' j  ^/ h       a(k,indxr(l))=a(k,indxc(l))
" J$ y, q$ e' ~: J  ^8 i7 x# d2 N       a(k,indxc(l))=dum
% c: o2 ~* r* u5 j; c1 V       enddo
# |- O) }, i0 a9 o# Q* B7 T/ G    endif
6 C- N8 t4 u, B0 x$ \6 S6 J; V    enddo
5 m/ Z7 T4 p  I2 A  x, |    end subroutine gaussj
2 d6 o) d1 |- R+ O! I101 end
</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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