Cholesky平方根算法C语言

#include
#include
#include
extern int flag=1; //0时正确 1时错误
int Cholesky(double **a , double *b , int H)
{
double *y = (double *) malloc(sizeof(double)*H);
double s; int i , j ,k ;
/*********** 求解L矩阵，L下三角 *************/
a[0][0] = sqrt(a[0][0]); //求解a11
for (i = 1; i < H; i++) //第一列
a[i][0] = a[i][0] / a[0][0]; //求解ai1
for (k = 1; k < H; k++) //第k列
{
for (i = k; i < H; i++) //下三角，i>k
{ //求aik
for (s = 0, j = 0; j < k; j++)
s += a[i][j] * a[k][j]; //求中间和
if (i > k) //下三角元素
a[i][k] = (a[i][k] - s) / a[k][k];
else if (i == k) //akk，主对角线元素
a[i][k] = sqrt(a[i][k] - s);
}
}
/*********** 求解L矩阵转置,L上三角 *************/
for(j = 0 ; j < H ; j++)
for(i = j+1 ; i < H ; i++)
a[j][i] = a[i][j]; //转置
/************** 检验方程是否错误 ***************/
for(i = 0 ,flag = 0 ; i < H ; j++)
{
if(a[i][i] == 0)
{
printf("方程系数输入有误\n");
flag = 1 ; //flag置1
}break;
}
if(flag != 1) // flag标志位不为0时继续计算
{/******************* 求解y[i] ******************/
y[0]= b[0] / a[0][0]; //先求y1
for (k = 0; k < H; k++)
{
for (s = 0, j = 0; j < k; j++)
s += a[k][j] * y[j];
y[k] = (b[k] - s) / a[k][k]; //再求yk...n
}
/******************* 求解x[i] ******************/
b[H-1]= y[H-1] / a[H-1][H-1]; //先求xn
for (k = H-2; k >= 0; k--)
{
for (s = 0, j = k+1; j < H; j++)
s += a[j][k] * b[j];
b[k] = (y[k] - s) / a[k][k]; //再求xi...1
}
/************** 打印L矩阵，L转置 ****************/
printf("\nL的下三角矩阵为：\n");
for(i = 0 ; i < H ; i++)
{
for(k = 0 ; k <= i ; k++)
printf("%lf ",a[i][k]);
printf("\n");
}
printf("\nL的转置矩阵为：\n");
for(j = 0 ; j< H ; j++)
{
for(i = 0 ; i < H ; i++)
{
if(i < j)
printf(" ");
else
printf("%lf ",a[i][j]);
}
printf("\n");
}
/***************** 打印Y的值 *******************/
printf("\nY的值为：\n");
for(i = 0 ; i < H ; i++)
printf("y%d=%lf\n",i+1,y[i]);
/***************** 打印X的值 *******************/
printf("\nX的值为：\n");
for(i = 0 ; i < H ; i++)
printf("x%d=%lf\n",i+1,b[i]);
free(y); //释放动态开辟的空间
return 0;
}
else if(flag == 1)
flag = 1; //flag标志位置1，退出
return 0;
}
int main()
{ double **a , *b;
int i, j , H;
while(flag == 1)
{ printf("请输入未知数个数n:\n");
scanf("%d", &H);
a = (double**)malloc(sizeof(double*)*H); //为二维数组分配H行
b = (double *) malloc(sizeof(dou

ble)*H);
for (i = 0; i < H; i++) //为每列分配H个大小空间
a[i] = (double*)malloc(sizeof(double)*H);
printf("请输入系数矩阵：\n"); //初始化
for (i = 0; i < H; i++)
for (j = 0; j < H; j++)
scanf("%lf",&a[i][j]);
printf("请输入向量b：\n");
for (i = 0; i < H; i++)
scanf("%lf",&b[i]);
Cholesky(a,b,H);
} free(b); //释放动态开辟的空间
for (i = 0; i < H; ++i)
free(a[i]); //a行
free(a); //a列
}