共轭梯度法C语言（西安交大）
#include
#include
#define N 10 /*定义矩阵阶数*/
void main()
{
int i,j,m,A[N][N],B[N];
double X[N],akv[N],dka[N],rk[N],dk[N],pk,pkk,ak,bk;
for(i=0;i<="">
for(j=0;j<n;j++)< p="">
{
if(i==j)
A[i][j]=-2;
else if(abs(i-j)==1)
A[i][j]=1;
else
A[i][j]=0;
}
}
for(i=0;i<="">
if(i==0||i==N-1)
B[i]=-1;
else
B[i]=0;
}
printf("\n"); /*输出系数矩阵A*/
printf("the Maxtr A is\n");
for(i=0;i<n;i++)< p="">
{
printf("\n");
for(j=0;j<n;j++)< p="">
printf("%3d",A[i][j]);
}
printf("\n"); /*输出矩阵B*/
printf("the Maxtr b is\n");
for(i=0;i<n;i++)< p="">
printf("%3d",B[i]);
printf("\n");
printf("input the Maxtr X\n"); /*给X0输入一个初始向量*/ for(i=0;i<n;i++)< p="">
X[i]=0;
printf("X chushi:\n");
for(i=0;i<="">
/*开始计算*/
for(i=0;i<="" p="">
{
pk=0.0;
for(j=0;j<n;j++)< p="">
pk+=A[i][j]*X[j];
正则化共轭梯度法
akv[i]=pk;
}
for(i=0;i<n;i++)< p="">
rk[i]=B[i]-akv[i];
for(i=0;i<="" p="">
dk[i]=rk[i];
for(m=0;m<="" p="">
{
for(i=0;i<="" p="">
{
pk=0.0;
for(j=0;j<n;j++)< p="">
pk+=A[i][j]*dk[j];
dka[i]=pk;
}
pk=0.0;pkk=0.0;
for(i=0;i<n;i++)< p="">
{
pk+=dka[i]*dk[i];
pkk+=rk[i]*rk[i];
}
if(pkk==0) /*如果残差pkk=0,计算结束*/ break;
ak=pkk/pk;
for(i=0;i<="" p="">
X[i]=X[i]+ak*dk[i];
for(i=0;i<="" p="">
{
pk=0.0;
for(j=0;j<n;j++)< p="">
pk+=A[i][j]*X[j];
akv[i]=pk;
}
for(i=0;i<n;i++)< p="">
rk[i]=B[i]-akv[i];
pk=0.0;
for(i=0;i<="" p="">
pk+=rk[i]*rk[i];
bk=pk/pkk;
for(i=0;i<="" p="">
dk[i]=rk[i]+bk*dk[i];
}
printf("\n");
printf("X cacualtate value is:\n"); /*输出运算结果X*/
for(i=0;i<n;i++)< p="">
printf("%f\n",X[i]);
}
</n;i++)<>
</n;i++)<>
</n;j++)<>
</n;i++)<>
</n;j++)<>
</n;i++)<>
</n;j++)<>
</n;i++)<>
</n;i++)<>
</n;j++)<>
</n;i++)<>
</n;j++)<>

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。