lsqr.c

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

/*  LSQR程序   */
#include<stdio.h>
#include<math.h>
#define M  5
#define N  5
void LSQR(float X[N])
{
  FILE *fp1,*fp2;
  int  i,j,k,n=0;
  char ch1[10],ch2[10];
  float A[M][N],B[N][M],D[M];
  float S[M],W[M],R[N],P[N],a,b,c,d,e,f,g,h;
  
/*参量说明：A[M][N]线性方程组系数矩阵；B为A的转置；X为解向量；D[M]为各个线性方程的常数项矩阵；
其余均为中间变量*/
     
  printf("请输入系数矩阵文件名：");
  scanf("%s",ch1);
  printf("请输入方程常数项矩阵文件名：");
  scanf("%s",ch2);
  printf("\n");
 if((fp1=fopen(ch1,"r"))==NULL)
  {
    printf("can not open file!\n");
    return;
  }
/*打开增广矩阵数据文件*/

if((fp2=fopen(ch2,"r"))==NULL)
  {
    printf("can not open file!\n\n");
    return;
  }
/*打开方程常数项矩阵数据文件*/
printf("  LSQX正在求解此方程组...\n");
printf("\n");

  for(i=0;i<M;i++)
  {
	  for(j=0;j<N;j++)
     {
       fscanf(fp1,"%f",&A[i][j]);B[j][i]=A[i][j];
     }
       fscanf(fp2,"%f",&D[i]);S[i]=D[i];
  }
   fclose(fp1);fclose(fp2);
 
 for(j=0;j<N;j++)
  {
    a=0;
   for(i=0;i<M;i++)
    {
     a=a+B[j][i]*D[i];
    }
    R[j] = a;
    P[j] = R[j];
    X[j] = 0;    /*解向量X赋初值为0*/
  }/*赋迭代初值*/

 do 
  {
    k=0;
    c=0;
    d=0;
     for(i=0;i<M;i++)
     {
       b=0;
       for(j=0;j<N;j++)
        {
           b=b+A[i][j]*P[j];
           
        }
        W[i] = b;
        d =(float) (d+pow(W[i],2));
      } 
	 for(j=0;j<N;j++)
	 {
      c=(float)(c+pow(R[j],2));
	 }
 
     e=(float)(c*pow(d,-1));
 
     for(j=0;j<N;j++)
     {
       X[j] = X[j]+e*P[j];
       if(fabs(e*P[j])>=0.00001) k=1; 
     }
	 for(i=0;i<M;i++)
	 {
      S[i]=S[i]-e*W[i];
	 }
     
     g=0;
     for(j=0;j<N;j++)
     {
       f=0;
       for(i=0;i<M;i++)
       {
         
         f = f+B[j][i]*S[i];
       }
       R[j] = f;
       g=(float)(g+pow(R[j],2));
     }
     h=(float)(g*pow(c,-1)); 
     for(j=0;j<N;j++)
     {
       P[j] = R[j]+h*P[j];  
     }n++;
  }while(k);
  printf("运算结束，运算结果为：\n");
  printf("迭代次数为：%d\n",n);
}

void main()
{
  FILE *fp;
  float X[N];
  int j;

  LSQR(X);	
  fp=fopen("LSQR_X.txt","w");
   printf("解向量为：\n");
   for(j=0;j<N;j++)
   {
     printf("%8.2f",X[j]);
     fprintf(fp,"%8.2f",X[j]);
   }
   printf("\n");
   fclose(fp);
 }