有用c语言编共轭梯度算法的么?

lorin

偶是新手
这里有用c语言编共轭梯度算法的么?
希望得到大家的帮忙
谢谢!!
lh_2777@sohu.com

[此贴子已经被作者于2005-11-22 20:14:10编辑过]

MVH

1. !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
   
2. !!!输入函数信息，输出函数的稳定点及迭代次数；
   
3. !!!iter整型变量，存放迭代次数；x0实型变量，开始存放进退法初始点；
   
4. !!!x,x1为n维变量，分别存放函数在第k、k+1次迭代点
   
5. !!!gradtf,gradts为n维实型变量，分别存放函数在第k、k+1次迭代点的梯度;
   
6. !!!dirf,dirs为n维实型变量，分别存放第k、k+1次搜索方向；
   
7. program main
   
8. real,dimension(:),allocatable::x,x1,gradtf,gradts,dirf,dirs,b
   
9. real,dimension(:,:),allocatable::hessin
   
10. real::x0,c,estol
    
11. integer::n,k,iter
    
12. print*,'请输入变量的维数'
    
13. read*,n
    
14. allocate (x(n),x1(n),gradtf(n),gradts(n),dirf(n),dirs(n),b(n))
    
15. allocate(hessin(n,n))
    
16. print*,'请输入初始点x'
    
17. read*,x
    
18. print*,'请输入hessin矩阵'
    
19. read*,hessin
    
20. print*,'请输入向量b'
    
21. read*,b
    
22. estol=0.000001
    
23. iter=0
    
24. 100 k=0
    
25. gradtf=matmul(hessin,x)+b
    
26. if(dot_product(gradtf,gradtf)<=estol)then
    
27. !print*,'函数的稳定点为：',x
    
28. !print*,'迭代次数为：',iter
    
29. goto 101
    
30. endif
    
31. dirf=(-1)*gradtf
    
32. 10x0=golden(x,dirf,hessin,b)
    
33. x1=x+x0*dirf
    
34. k=k+1
    
35. iter=iter+1
    
36. if(iter>10*n)then
    
37. print*,"out"
    
38. goto 101
    
39. endif
    
40. print*,"第",iter,"次运行结果为","方向为",dirf,"步长为",x0
    
41. print*,x1,"f(x)=",f(x1,hessin,b)
    
42. gradts=matmul(hessin,x1)+b
    
43. if(dot_product(gradts,gradts)<=estol)then
    
44. !print*,'函数的稳定点为：',x1
    
45. !print*,'迭代次数为：',iter
    
46. goto 101
    
47. endif
    
48. if(k==n)then
    
49. x=x1
    
50. goto 100
    
51. else
    
52. c=dot_product(gradts,gradts)/dot_product(gradtf,gradtf)
    
53. dirs=(-1)*gradts+c*dirf
    
54. dirf=dirs
    
55. if(dot_product(dirf,gradts)>0)then
    
56. x=x1
    
57. goto 100
    
58. else
    
59. goto 10
    
60. endif
    
61. endif
    
62. 
    
63. contains

复制代码

1. !!!子程序，返回函数值
   
2. function f(x,A,b) result(f_result)
   
3. real,dimension(:),intent(in)::x,b
   
4. real,dimension(:,:),intent(in)::A
   
5. real::f_result
   
6. f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
   
7. end function f

复制代码

1. !!!精确线搜索0.618法子程序，返回迭代步长
   
2. function golden(x,d,A,b) result(golden_n)
   
3. real::golden_n
   
4. real::x0
   
5. real,dimension(:),intent(in)::x,d
   
6. real,dimension(:),intent(in)::b
   
7. real,dimension(:,:),intent(in)::A
   
8. real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
   
9. parameter(r=0.618)
   
10. tol=0.0001
    
11. dx=0.1
    
12. x0=1
    
13. x1=x0+dx
    
14. f0=f(x+x0*d,A,b)
    
15. f1=f(x+x1*d,A,b)
    
16. if(f04 dx=dx+dx
    
17. x2=x0-dx
    
18. f2=f(x+x2*d,A,b)
    
19. if(f2 x1=x0
    
20. x0=x2
    
21. f1=f0
    
22. f0=f2
    
23. goto 4
    
24. else
    
25. a1=x2
    
26. b1=x1
    
27. endif
    
28. else
    
29. 2 dx=dx+dx
    
30. x2=x1+dx
    
31. f2=f(x+x2*d,A,b)
    
32. if(f2>=f1)then
    
33. b1=x2
    
34. a1=x0
    
35. else
    
36. x0=x1
    
37. x1=x2
    
38. f0=f1
    
39. f1=f2
    
40. goto 2
    
41. endif
    
42. endif
    
43. x1=a1+(1-r)*(b1-a1)
    
44. x2=a1+r*(b1-a1)
    
45. f1=f(x+x1*d,A,b)
    
46. f2=f(x+x2*d,A,b)
    
47. 3 if(abs(b1-a1)<=tol)then
    
48. x0=(a1+b1)/2
    
49. else
    
50. if(f1>f2)then
    
51. a1=x1
    
52. x1=x2
    
53. f1=f2
    
54. x2=a1+r*(b1-a1)
    
55. f2=f(x+x2*d,A,b)
    
56. goto 3
    
57. else
    
58. b1=x2
    
59. x2=x1
    
60. f2=f1
    
61. x1=a1+(1-r)*(b1-a1)
    
62. f1=f(x+x1*d,A,b)
    
63. goto 3
    
64. endif
    
65. endif
    
66. golden_n=x0
    
67. endfunction golden
    
68. 101 end program main

复制代码

我只有Fortran写的

MVH

不好意思，这个版主发过了
http://forum.vibunion.com/thread-5540-1-1.html

frogfish

C常用算法程序集上有

lorin

我不会Fortran
仍然谢谢你perl

风花雪月

1. #include "math.h"
   
2. #include"iostream.h"
   
3. void grad(int n,double *a,double *b,double eps,double *x)
   
4. {
   
5. int i,j,k;
   
6. double *p,*r,*s,*q,alpha,beta,d,e;
   
7. p=new double[n];
   
8. r=new double[n];
   
9. s=new double[n];
   
10. q=new double[n];
    
11. for (i=0; i<=n-1; i++)
    
12. {
    
13. x=0.0; p=b; r=b;
    
14. }
    
15. i=0;
    
16. while (i<=n-1)
    
17. {
    
18. for (k=0; k<=n-1; k++)
    
19. {
    
20. s[k]=0.0;
    
21. for (j=0; j<=n-1; j++)
    
22. s[k]=s[k]+a[k*n+j]*p[j];
    
23. }
    
24. d=0.0; e=0.0;
    
25. for (k=0; k<=n-1; k++)
    
26. {
    
27. d=d+p[k]*b[k]; e=e+p[k]*s[k];
    
28. }
    
29. alpha=d/e;
    
30. for (k=0; k<=n-1; k++)
    
31. x[k]=x[k]+alpha*p[k];
    
32. for (k=0; k<=n-1; k++)
    
33. {
    
34. q[k]=0.0;
    
35. for (j=0; j<=n-1; j++)
    
36. q[k]=q[k]+a[k*n+j]*x[j];
    
37. }
    
38. d=0.0;
    
39. for (k=0; k<=n-1; k++)
    
40. {
    
41. r[k]=b[k]-q[k]; d=d+r[k]*s[k];
    
42. }
    
43. beta=d/e; d=0.0;
    
44. for (k=0; k<=n-1; k++) d=d+r[k]*r[k];
    
45. d=sqrt(d);
    
46. if (d<eps)
    
47. {
    
48. delete []p;
    
49. delete []r;
    
50. delete []s;
    
51. delete []q;
    
52. return;
    
53. }
    
54. for (k=0; k<=n-1; k++)
    
55. p[k]=r[k]-beta*p[k];
    
56. i=i+1;
    
57. }
    
58. delete []p;
    
59. delete []r;
    
60. delete []s;
    
61. delete []q;
    
62. return;
    
63. }

复制代码

1. void main()
   
2. {
   
3. int i;
   
4. double eps,x[4];
   
5. static double a[4][4]={
   
6. {5.0,7.0,6.0,5.0},
   
7. {7.0,10.0,8.0,7.0},
   
8. {6.0,8.0,10.0,9.0},
   
9. {5.0,7.0,9.0,10.0}};
   
10. static double b[4]={23.0,32.0,33.0,31.0};
    
11. eps=0.000001;
    
12. grad(4,*a,b,eps,x);
    
13. for (i=0; i<=3; i++)
    
14. cout<<"x["<<i<<"]="<<x<<endl;
    
15. }

复制代码

lorin

谢谢大家！！
谢谢风花雪月

[此贴子已经被作者于2005-11-22 20:41:47编辑过]