levin.h

用C开发的莱文森算法

//t：对称LEVIN矩阵元素
//b:方程组右端常数
//X:解
int levin(double t[],double b[],int n,double x[])
{ int i,j,k;
  double a,beta,q,c,h,*y,*s;
  
  s=(double *)malloc(n*sizeof(double));
  y=(double *)malloc(n*sizeof(double));

  a=t[0];
  if((fabs(a)+1)==1.0)
  { free(s);
    free(y);
	printf("ill!\n");
	return -1;
  }

  y[0]=1.0;
  x[0]=b[0]/a;
  for(k=1;k<n;k++)
  { beta=0;
    q=0;
	 for(j=0;j<k;j++)
	 { 
	   beta+=y[j]*t[j+1];
	   q+=x[j]*t[k-j];

	 }

   if((fabs(a)+1)==1.0)
   { free(s);
    free(y);
	printf("ill!\n");
	return -1;
   }
	 
   c=-beta/a;
   s[0]=c*y[k-1];
   y[k]=y[k-1];

   if(k!=1)
   { for(i=1;i<k;i++)
      s[i]=y[i-1]+c*y[k-i-1];
   }

   s[k]=y[k-1];
   a+=c*beta;
   
   if((fabs(a)+1)==1.0)
   { free(s);
    free(y);
	printf("ill!\n");
	return -1;
   }
   
    h=(b[k]-q)/a;
 
    for(i=0;i<k;i++)
	{ x[i]+=h*s[i];
	  y[i]=s[i];
	}
    
	x[k]=h*y[k];
  }
   free(s);
   free(y);

   return(1);


}