brinv0.c

优化计算的新方法 可以直接进行数值计算

  #include "stdio.h" 
  #include "io.h"
  #include "stdlib.h"
  #include "math.h"
  #include <sys/types.h>
  #include <sys/stat.h>
  #include <iostream.h>


FILE *fp1,*fp2,*fp3;  

void main()
  {
	  float *um,*pm,*qm;
	  long arraysize;
	  long filelengh;
	  float z1=0,z2=0,z3=0,z4=0,z5=0,q1=0,q2=0,q3=0,b=0;
	  int fh;
	  
	 int i,k,ys;
	 float w[5],z[5];
	 int j,jj;
	 float qs;
  
    static float a[4][10];
	
	if((fp3=fopen("fp3.xls","w"))==NULL)
	      { printf("w3?\n");   exit(0); }
	  if((fp2=fopen("fp2.xls","w"))==NULL)
	      { printf("w2 ?\n");  exit(0); }
	  if((fh = _open( "tt41.xls",   _S_IREAD))  != -1 )
	  {
         filelengh=_filelength(fh);
		 _close( fh );
	  }
	  arraysize=filelengh/28;

	  if((fp1=fopen("tt41.xls","r"))==NULL)
	      { printf("r1?\n");   exit(0); }
	
   	if((um=(float *)malloc(arraysize*sizeof(float)))==NULL)
	{
		printf("can't allocate more memory,terminating.\n");
		exit(1);
	}
	if((pm=(float *)malloc(arraysize*sizeof(float)))==NULL)
	{
		printf("can't allocate more memory,terminating.\n");
		exit(1);
	}
	if((qm=(float *)malloc(arraysize*sizeof(float)))==NULL)
	{
		printf("can't allocate more memory,terminating.\n");
		exit(1);
	}
	for(i=0;i<arraysize;i++)
		fscanf(fp1,"%f\t%f\t%f\n",&um[i],&pm[i],&qm[i]);
	for(i=0;i<arraysize;i++)
	{
		z1=z1+um[i]*um[i]*um[i]*um[i];
		z2=z2+um[i]*um[i]*um[i];
		z3=z3+um[i]*um[i];
		z4=z4+um[i];
		z5=z5+1;
		q1=q1+um[i]*um[i]*pm[i];
		q2=q2+um[i]*pm[i];
		q3=q3+pm[i];
	}
	z1=z1*(-2);
	z2=z2*(-2);
	z3=z3*(-2);
	z4=z4*(-2);
	z5=z5*(-2);
	q1=q1*(-2);
	q2=q2*(-2);
	q3=q3*(-2);

	a[0][0]=1;a[0][1]=0;a[0][2]=0;a[0][3]=0;a[0][4]=z1;a[0][5]=z2;a[0][6]=z3;a[0][7]=1;a[0][8]=-1;a[0][9]=q1;
	a[1][0]=0;a[1][1]=1;a[1][2]=0;a[1][3]=0;a[1][4]=z2;a[1][5]=z3;a[1][6]=z4;a[1][7]=1;a[1][8]=-1;a[1][9]=q2;
	a[2][0]=0;a[2][1]=0;a[2][2]=1;a[2][3]=0;a[2][4]=z3;a[2][5]=z4;a[2][6]=z5;a[2][7]=1;a[2][8]=-1;a[2][9]=q3;
	a[3][0]=0;a[3][1]=0;a[3][2]=0;a[3][3]=1;a[3][4]=-1;a[3][5]=-1;a[3][6]=-1;a[3][7]=0;a[3][8]=-1;a[3][9]=-1;

	for(i=0;i<4;i++)
	{
		if(a[i][9]<0)
			goto loop1;
	}
	w[1]=a[0][9];w[2]=a[1][9];w[3]=a[2][9];w[4]=a[3][9];z[1]=0;z[2]=0;z[3]=0,z[4]=0;
	goto end;
loop1:
	for(i=0;i<4;i++)//找最大的-q
	  	if(-a[i][9]>qs)
		{
			qs=-a[i][9];
			j=i;//j记录主元所在的行
		}
	for(i=0;i<10;i++)//对应的j行为主行，z0对应的列为主列，进行主元消去
		a[j][i]=-a[j][i];
	for(i=0;i<j;i++)
		for(k=0;k<10;k++)
			a[i][k]=a[i][k]+a[j][k];
    for(i=j+1;i<4;i++)
		for(k=0;k<10;k++)
			a[i][k]=a[i][k]+a[j][k];
		ys=4+j;
loop2://第二步
	for(i=0;i<4;i++)
		if(a[i][ys]>0)
			goto loop22;
//loop21:
    //w[1]=a[0][9];w[2]=a[1][9];w[3]=a[2][9];w[4]=a[3][9];z[1]=0;z[2]=0;z[3]=0,z[4]=0;
	goto end;
loop22:
	qs=a[i][9]/a[i][ys];
	j=i;
	for(i=0;i<4;i++)
		if(a[i][ys]>0&&a[i][9]/a[i][ys]<qs)
		{
		   qs=a[i][9]/a[i][ys];//按最小比值规则确定指标r
		   j=i;//记录主元所在的行
		}
	//j=j1;
	if(a[j][8]==1)
	{
		for(i=0;i<j;i++)
			if(a[i][8]=!0)
			   goto loop3;
		for(i=j+1;i<4;i++)
			if(a[i][8]=!0)
		       goto loop3;
		goto loop4;
	}
loop3://第三步
	//确定第j行的基变量
	for(k=0;k<9;k++)
	{
		if(a[j][k]==1)
		{
			for(i=0;i<j;i++)
				if(a[i][k]!=0)
					goto next;
			for(i=j+1;i<4;i++)
				if(a[i][k]!=0)
					goto next;
			goto loop31;
		}
next:;
	}
loop31:
	jj=k;//记录j行的基变量所在的列，即离基变量所在的列
	b=a[j][ys];//保存主元的值
	for(i=0;i<10;i++)
		a[j][i]=a[j][i]/b;
	for(i=0;i<j;i++)
		for(k=0;k<10;k++)
			a[i][k]=a[i][k]-a[i][k]*a[j][k];
	for(i=j+1;i<4;i++)
		for(k=0;k<10;k++)
			a[i][k]=a[i][k]-a[i][k]*a[j][k];
	if(jj<4)
		ys=4+jj;
	else
		ys=jj-4;
	goto loop2;
loop4://第四步
	for(i=0;i<10;i++)
		a[j][i]=a[j][i]/a[j][ys];

	for(i=0;i<j;i++)
		for(k=0;k<10;k++)
			a[i][k]=a[i][k]-a[i][k]*a[j][k];
	for(i=j+1;i<4;i++)
		for(k=0;k<10;k++)
			a[i][k]=a[i][k]-a[i][k]*a[j][k];
end:
	for(k=0;k<5;k++)
		z[k]=0;
	for(i=4;i<8;i++)//搜索z1->z4
	{
		if(a[0][i]==1&&a[1][i]==0&&a[2][i]==0&&a[3][i]==0)
		{
			z[i-3]=a[0][9];
			goto next1;
		}
		if(a[0][i]==0&&a[1][i]==1&&a[2][i]==0&&a[3][i]==0)
		{
			z[i-3]=a[1][9];
			goto next1;
		}
        if(a[0][i]==0&&a[1][i]==0&&a[2][i]==1&&a[3][i]==0)
		{
			z[i-3]=a[2][9];
			goto next1;
		}
		if(a[0][i]==0&&a[1][i]==0&&a[2][i]==0&&a[3][i]==1)
	       z[i-3]=a[3][9];
next1:;
	}

	fprintf(fp3,"参数a     ,b    ,c    分别为：\n");
	fprintf(fp3,"%f\t%f\t%f\n",z[1],z[2],z[3]);
	for(i=0;i<arraysize;i++)
	    fprintf(fp2,"%f\n",um[i]*um[i]*z[1]+um[i]*z[2]+z[3]);
	free(um);
	free(pm);
	free(qm);		

	return;   
    
  }