有用c语言编共轭梯度算法的么?
偶是新手这里有用c语言编共轭梯度算法的么?
希望得到大家的帮忙
谢谢!!
lh_2777@sohu.com
[此贴子已经被作者于2005-11-22 20:14:10编辑过]
!!!本程序适用于求解形如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)<=estol)then
!print*,'函数的稳定点为：',x
!print*,'迭代次数为：',iter
goto 101
endif
dirf=(-1)*gradtf
10x0=golden(x,dirf,hessin,b)
x1=x+x0*dirf
k=k+1
iter=iter+1
if(iter>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)<=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)>0)then
x=x1
goto 100
else
goto 10
endif
endif
contains
!!!子程序，返回函数值
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
!!!精确线搜索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(f04 dx=dx+dx
x2=x0-dx
f2=f(x+x2*d,A,b)
if(f2 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>=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)<=tol)then
x0=(a1+b1)/2
else
if(f1>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
endfunction golden
101 end program main
我只有Fortran写的 不好意思，这个版主发过了
http://forum.vibunion.com/thread-5540-1-1.html C常用算法程序集上有 我不会Fortran
仍然谢谢你perl #include "math.h"
#include"iostream.h"
void grad(int n,double *a,double *b,double eps,double *x)
{
int i,j,k;
double *p,*r,*s,*q,alpha,beta,d,e;
p=new double;
r=new double;
s=new double;
q=new double;
for (i=0; i<=n-1; i++)
{
x=0.0; p=b; r=b;
}
i=0;
while (i<=n-1)
{
for (k=0; k<=n-1; k++)
{
s=0.0;
for (j=0; j<=n-1; j++)
s=s+a*p;
}
d=0.0; e=0.0;
for (k=0; k<=n-1; k++)
{
d=d+p*b; e=e+p*s;
}
alpha=d/e;
for (k=0; k<=n-1; k++)
x=x+alpha*p;
for (k=0; k<=n-1; k++)
{
q=0.0;
for (j=0; j<=n-1; j++)
q=q+a*x;
}
d=0.0;
for (k=0; k<=n-1; k++)
{
r=b-q; d=d+r*s;
}
beta=d/e; d=0.0;
for (k=0; k<=n-1; k++) d=d+r*r;
d=sqrt(d);
if (d<eps)
{
delete []p;
delete []r;
delete []s;
delete []q;
return;
}
for (k=0; k<=n-1; k++)
p=r-beta*p;
i=i+1;
}
delete []p;
delete []r;
delete []s;
delete []q;
return;
}
void main()
{
int i;
double eps,x;
static double a={
{5.0,7.0,6.0,5.0},
{7.0,10.0,8.0,7.0},
{6.0,8.0,10.0,9.0},
{5.0,7.0,9.0,10.0}};
static double b={23.0,32.0,33.0,31.0};
eps=0.000001;
grad(4,*a,b,eps,x);
for (i=0; i<=3; i++)
cout<<"x["<<i<<"]="<<x<<endl;
} 谢谢大家！！
谢谢风花雪月
[此贴子已经被作者于2005-11-22 20:41:47编辑过]
页:
[1]