svd.c

由c语言编写的解非线性方程

#include<stdio.h>
#include"svdqj.c"
#define M   5
#define N   5
#define ka  6 // ka=max(M,N)+1
void main()
{
    FILE *fp1,*fp2;FILE *fp3;
	int i,j;
	double a[M][N],b[M],x[N],aa[N][M],u[M][M],v[N][N],eps=0.0001,s;
    char ch1[20],ch2[20],flag1,flag2;
	printf("\n");
	printf("*******************************************************************\n\n");
    printf("                  SVD奇异值分解求解线性方程组\n\n");
	printf("*******************************************************************\n\n");
 
	do
    {
		flag1='a';
	printf("  1.请输入存放线性方程组系数矩阵A的文件名：");
    scanf("%s",ch1);
    printf("  2.请输入存放方程组右端向量b的文件名：");
    scanf("%s",ch2);
    printf("\n");
    if((fp1=fopen(ch1,"r"))==NULL)
	{
      printf("没有找到目标文件%s,是否重新输入文件名<Y/N>?",ch1);
      getchar();
      flag1=getchar();
      
	}
    if((fp2=fopen(ch2,"r"))==NULL)
	{
      printf("没有找到目标文件%s,是否重新输入文件名<Y/N>?",ch2);
      getchar();
      flag1=getchar();
	  
	}
	}while((flag1 == 'y')||(flag1 == 'Y'));
	if((flag1 == 'n')||(flag1 == 'N')) return;
    flag1='y';
    for(i=0;i<M;i++)
	{
	   for(j=0;j<N;j++)
	   {
         fscanf(fp1,"%lf",&a[i][j]);
	   }
        fscanf(fp2,"%lf",&b[i]);
	}
   fclose(fp1);fclose(fp2);

    printf("MAT A(%d*%d)=UWV IS:\n",M,N);
    for(i=0;i<M;i++)
	{
		for(j=0;j<N;j++)
			printf("%13.4f",a[i][j]);
		printf("\n");
	}
	printf("vector b is:\n");
	for(i=0;i<M;i++)
		printf("%13.4f",b[i]);
    printf("\n\n");
    printf("******文件数据读入是否正确<Y/N>?******\n");
	printf("  Y:  文件数据正确，程序继续运行\n");
	printf("  N:  文件数据错误，返回修改数据\n");
    printf("<command>：");
    getchar();
    flag1=getchar();
	
	if((flag1 != 'Y')&&(flag1 != 'y')) return;
    printf("\n");
	printf("对A进行奇异值分解...显示分解结果<Y/N>?");
	getchar();
    flag2=getchar();
	i=ginv(a,M,N,aa,eps,u,v,ka);
	
	if((flag2 == 'y')||(flag2 == 'Y'))
    {
	printf("MAT U(%d*%d) IS:\n",M,M);
	for(i=0;i<M;i++)
	{
		for(j=0;j<M;j++)
		   printf("%13.4e",u[i][j]);
		printf("\n");
	}
	printf("\n");
    printf("MAT V(%d*%d) IS:\n",N,N);
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		   printf("%13.4e",v[i][j]);
		printf("\n");
	}
	printf("\n");
    printf("MAT W(%d*%d) IS:\n",M,N);
	for(i=0;i<M;i++)
	{
		for(j=0;j<N;j++)
		   printf("%13.4e",a[i][j]);
		printf("\n");
	}
    printf("\n");
    printf("MAT A+(%d*%d) is:\n",N,M);
	if(i>0)
	{
        for(i=0;i<N;i++)
		{
	 	  for(j=0;j<M;j++)
			printf("%13.4f",aa[i][j]);
		 printf("\n");
		}
		
	}
	}
     printf("\n");
	for(i=0;i<N;i++)
	{
		s=0.0;
		for(j=0;j<M;j++)
		{ s=s+aa[i][j]*b[j];}
		x[i]=s;
	}
	printf("方程组AX=b的解向量X为:\n\n");
    fp3=fopen("svd_x.txt","w");
	for(i=0;i<N;i++)
	{
		printf("%13.4f",x[i]);fprintf(fp3,"%8.2f",x[i]);
	}
	printf("\n");fclose(fp3);
}