共轭梯度法C语言(西安交大)

1、精选优质文档-倾情为你奉上#include<stdio.h>#include<math.h>#define N 10 /*定义矩阵阶数*/void main()int i,j,m,ANN,BN;double XN,akvN,dkaN,rkN,dkN,pk,pkk,ak,bk;for(i=0;i<N;i+) /*输入系数矩阵A*/for(j=0;j<N;j+)if(i=j)Aij=-2;else if(abs(i-j)=1) Aij=1;else Aij=0; for(i=0;i<N;i+) /*输入矩阵B*/ if(i=0|i=N-1) Bi=-1;

2、else Bi=0; printf("n"); /*输出系数矩阵A*/ printf("the Maxtr A isn"); for(i=0;i<N;i+) printf("n"); for(j=0;j<N;j+) printf("%3d",Aij); printf("n"); /*输出矩阵B*/printf("the Maxtr b isn");for(i=0;i<N;i+) printf("%3d",Bi);printf("

3、n"); printf("input the Maxtr Xn"); /*给X0输入一个初始向量*/ for(i=0;i<N;i+) Xi=0; printf("X chushi:n"); for(i=0;i<N;i+) /*显示给定的初始向量X0*/ printf("%3.0f",Xi); /*开始计算*/ for(i=0;i<N;i+) /*计算rk*/ pk=0.0; for(j=0;j<N;j+) pk+=Aij*Xj; akvi=pk; for(i=0;i<N;i+) rki=Bi-a

4、kvi; for(i=0;i<N;i+) /*给rO赋值*/ dki=rki; for(m=0;m<N;m+) /*开始迭代循环*/ for(i=0;i<N;i+) /*计算ak*/ pk=0.0; for(j=0;j<N;j+) pk+=Aij*dkj; dkai=pk; pk=0.0;pkk=0.0; for(i=0;i<N;i+) pk+=dkai*dki; pkk+=rki*rki;if(pkk=0) /*如果残差pkk=0,计算结束*/ break; ak=pkk/pk; for(i=0;i<N;i+) /*计算X(k+1)*/ Xi=Xi+ak*dki; for(i=0;i<N;i+) /*计算r(k+1)*/ pk=0.0; for(j=0;j<N;j+) pk+=Aij*Xj; akvi=pk; for(i=0;i<N;i+) rki=Bi-akvi; pk=0.0; for(i=0;i<N;i+) /*计算bk*/ pk+=rki*rki; bk=pk/pkk; for(i=0;i<N;i+) /*计算d(k+1)*/ dki=rki+bk*dki; printf("n"); printf("X cacualtate value is