fortran 90共轭梯度算法.c

用于科学计算的Fortran 90算法源程序

Fortran 90共轭梯度算法  
 

<P>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
    !!!输入函数信息，输出函数的稳定点及迭代次数；
    !!!iter整型变量，存放迭代次数；x0实型变量，开始存放进退法初始点；
    !!!x,x1为n维变量，分别存放函数在第k、k+1次迭代点
    !!!gradtf,gradts为n维实型变量，分别存放函数在第k、k+1次迭代点的梯度;
    !!!dirf,dirs为n维实型变量，分别存放第k、k+1次搜索方向；
    program main
    real,dimension(:),allocatable::x,x1,gradtf,gradts,dirf,dirs,b
    real,dimension(:,:),allocatable::hessin
    real::x0,c,estol
    integer::n,k,iter
    print*,'请输入变量的维数'
    read*,n
    allocate (x(n),x1(n),gradtf(n),gradts(n),dirf(n),dirs(n),b(n))
    allocate(hessin(n,n))
    print*,'请输入初始点x'
    read*,x
    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入向量b'     
    read*,b
    estol=0.000001
    iter=0
100 k=0
    gradtf=matmul(hessin,x)+b
    if(dot_product(gradtf,gradtf)&lt;=estol)then
        !print*,'函数的稳定点为：',x
  !print*,'迭代次数为：',iter
     goto 101
    endif
    dirf=(-1)*gradtf
10  x0=golden(x,dirf,hessin,b)   
    x1=x+x0*dirf
k=k+1
iter=iter+1
if(iter&gt;10*n)then
     print*,"out"
  goto 101
    endif
print*,"第",iter,"次运行结果为","方向为",dirf,"步长为",x0
print*,x1,"f(x)=",f(x1,hessin,b)
    gradts=matmul(hessin,x1)+b 
if(dot_product(gradts,gradts)&lt;=estol)then
    !print*,'函数的稳定点为：',x1
    !print*,'迭代次数为：',iter
    goto 101
endif
    if(k==n)then
    x=x1
    goto 100
else
    c=dot_product(gradts,gradts)/dot_product(gradtf,gradtf)
    dirs=(-1)*gradts+c*dirf
    dirf=dirs
    if(dot_product(dirf,gradts)&gt;0)then
       x=x1
    goto 100
    else
       goto 10
    endif
endif
       
   contains</P>
<P>    !!!子程序，返回函数值
    function f(x,A,b) result(f_result)
    real,dimension(:),intent(in)::x,b
    real,dimension(:,:),intent(in)::A
    real::f_result
       f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
    end function f</P>
<P>    !!!精确线搜索0.618法子程序，返回迭代步长
    function golden(x,d,A,b) result(golden_n)
    real::golden_n
    real::x0
    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
    dx=0.1
x0=1
    x1=x0+dx
    f0=f(x+x0*d,A,b)
    f1=f(x+x1*d,A,b)
    if(f0&lt;f1)then
4       dx=dx+dx
        x2=x0-dx
        f2=f(x+x2*d,A,b)
        if(f2&lt;f0)then
           x1=x0
        x0=x2
        f1=f0
        f0=f2
        goto 4
        else
           a1=x2
        b1=x1
        endif
    else
2       dx=dx+dx
        x2=x1+dx
        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
            b1=x2
         a1=x0
        else
            x0=x1
         x1=x2
         f0=f1
         f1=f2
         goto 2
        endif
    endif
    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)
3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
    else
        if(f1&gt;f2)then
        a1=x1
        x1=x2
        f1=f2
        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
        goto 3
     else
        b1=x2
        x2=x1
        f2=f1
        x1=a1+(1-r)*(b1-a1)
        f1=f(x+x1*d,A,b)
        goto 3
     endif
    endif
    golden_n=x0
    end  function golden
101 end program main</P>
<P>本程序由Fortran90编写，在Vistual Fortran 5上调试通过！希望大家批评指正！</P>