sirt.h

数值计算利用fft分解法来求解线性方程组,简明方便

//SIRT
void SIRT(double *a,double *T,double *s,double S0,int n,int m,double e)
{
    int i,j,k,*M=new int[m];
	double *ss=new double[m],temp0,temp1,temp2;
	bool judge=true;
    for(i=0;i<m;i++)    
	{
      s[i]=S0;
	}
        for(i=0;i<m;i++)
		{   M[i]=0;
	     	for(j=0;j<n;j++)
			{
             if(a[j*m+i]!=0.0)
				 M[i]++;
			}
		}
    while(judge==true)
	{
         judge=false;
	  for(j=0;j<m;j++)
	  {
		ss[j]=s[j];temp0=0.0;
		for(i=0;i<n;i++)
		{           
			temp1=0.0;temp2=0.0;
			for(k=0;k<m;k++)
			{
		    	temp1=temp1+a[i*m+k]*a[i*m+k];
				temp2=temp2+a[i*m+k]*s[k];
			}
		    temp0=temp0+(T[i]-temp2)*a[i*m+j]/temp1;
		}
		s[j]=s[j]+temp0/M[j];
		if(fabs(s[j]-ss[j])>e)
			judge=true;
		if(s[j]>1/1000.0)
			s[j]=1/1000.0;
		else if(s[j]<1/5000.0)
			s[j]=1/5000.0;
	  }

	}
}
void ART(double *a,double *T,double *s,double S0,int n,int m,double e,double R)
{
    int i,j,k=0;
	double temp1,temp2;
	bool judge=true;
    for(i=0;i<m;i++)    
	{
      s[i]=S0;
	}
while(judge==true)
{
    judge=false;
    i=k%n;
    for(j=0;j<m;j++)
	{
		temp1=temp2=0.0;
		for(int kk=0;kk<m;kk++)
		{
			temp1=temp1+a[i*m+kk]*s[kk];
			temp2=temp2+a[i*m+kk]*a[i*m+kk];
		}
		s[j]=s[j]+R*(T[i]-temp1)/(temp2)*a[i*m+j];
		if(fabs(R*(T[i]-temp1)/(temp2)*a[i*m+j])>e)
			judge=true;
        if(s[j]>1/1000.0)
			s[j]=1/1000.0;
		else if(s[j]<1/5000.0)
			s[j]=1/5000.0;
	}

	k++;
}

}