ch1.h

数值处理算法程序

        { t=b[i*il];
          for (j=0;j<=m-1;j++)
            { u=i*m+j; v=k*m+j;
              d[u]=d[u]-t*d[v];
            }
          for (j=1;j<=il-1;j++)
            { u=i*il+j; v=k*il+j;
              b[u-1]=b[u]-t*b[v];
            }
          u=i*il+il-1; b[u]=0.0;
        }
      if (ls!=(n-1)) ls=ls+1;
    }
  p=b[(n-1)*il];
  if (fabs(p)+1.0==1.0)
    { printf("fail\n"); return(0);}
  for (j=0;j<=m-1;j++)
    { u=(n-1)*m+j; d[u]=d[u]/p;}
  ls=1;
  for (i=n-2;i>=0;i--)
    { for (k=0;k<=m-1;k++)
        { u=i*m+k;
          for (j=1;j<=ls;j++)
            { v=i*il+j; is=(i+j)*m+k;
              d[u]=d[u]-b[v]*d[is];
            }
        }
      if (ls!=(il-1)) ls=ls+1;
    }
  return(2);
}
/////////////////////////////////////////////////////////////
int aldle(double a[],int n,int m,double c[])
{ 
	int i,j,l,k,u,v,w,k1,k2,k3;
	double p;
	if (fabs(a[0])+1.0==1.0)
    {
		printf("fail\n"); return(-2);
	}
	for (i=1; i<=n-1; i++)
    { 
		u=i*n; a[u]=a[u]/a[0];
	}
	for (i=1; i<=n-2; i++)
    {
		u=i*n+i;
		for (j=1; j<=i; j++)
		{ 
			v=i*n+j-1; l=(j-1)*n+j-1;
			a[u]=a[u]-a[v]*a[v]*a[l];
		}
		p=a[u];
		if (fabs(p)+1.0==1.0)
        { 
			printf("fail\n"); return(-2);
		}
	    for (k=i+1; k<=n-1; k++)
        {
			u=k*n+i;
			for (j=1; j<=i; j++)
            {
				v=k*n+j-1; l=i*n+j-1; w=(j-1)*n+j-1;
				a[u]=a[u]-a[v]*a[l]*a[w];
            }
			a[u]=a[u]/p;
        }
    }
	u=n*n-1;
	for (j=1; j<=n-1; j++)
    {
		v=(n-1)*n+j-1; w=(j-1)*n+j-1;
		a[u]=a[u]-a[v]*a[v]*a[w];
    }
  p=a[u];
  if (fabs(p)+1.0==1.0)
  { 
	  printf("fail\n"); return(-2);
  }
  for (j=0; j<=m-1; j++)
	for (i=1; i<=n-1; i++)
    {
		u=i*m+j;
		for (k=1; k<=i; k++)
        {
			v=i*n+k-1; w=(k-1)*m+j;
			c[u]=c[u]-a[v]*c[w];
        }
    }
  for (i=1; i<=n-1; i++)
  {
	  u=(i-1)*n+i-1;
      for (j=i; j<=n-1; j++)
      {
		  v=(i-1)*n+j; w=j*n+i-1;
          a[v]=a[u]*a[w];
      }
  }
  for (j=0; j<=m-1; j++)
  {
	  u=(n-1)*m+j;
      c[u]=c[u]/p;
      for (k=1; k<=n-1; k++)
      {
		  k1=n-k; k3=k1-1; u=k3*m+j;
          for (k2=k1; k2<=n-1; k2++)
          {
			  v=k3*n+k2; w=k2*m+j;
              c[u]=c[u]-a[v]*c[w];
          }
          c[u]=c[u]/a[k3*n+k3];
      }
  }
  return(2);
}
/////////////////////////////////////////////////////////////
int achol(double a[],int n,int m,double d[])
{ int i,j,k,u,v;
  if ((a[0]+1.0==1.0)||(a[0]<0.0))
    { printf("fail\n"); return(-2);}
  a[0]=sqrt(a[0]);
  for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];
  for (i=1; i<=n-1; i++)
    { u=i*n+i;
      for (j=1; j<=i; j++)
        { v=(j-1)*n+i;
          a[u]=a[u]-a[v]*a[v];
        }
      if ((a[u]+1.0==1.0)||(a[u]<0.0))
        { printf("fail\n"); return(-2);}
      a[u]=sqrt(a[u]);
      if (i!=(n-1))
        { for (j=i+1; j<=n-1; j++)
            { v=i*n+j;
              for (k=1; k<=i; k++)
                a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
              a[v]=a[v]/a[u];
            }
        }
    }
  for (j=0; j<=m-1; j++)
    { d[j]=d[j]/a[0];
      for (i=1; i<=n-1; i++)
        { u=i*n+i; v=i*m+j;
          for (k=1; k<=i; k++)
            d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
          d[v]=d[v]/a[u];
        }
    }
  for (j=0; j<=m-1; j++)
    { u=(n-1)*m+j;
      d[u]=d[u]/a[n*n-1];
      for (k=n-1; k>=1; k--)
        { u=(k-1)*m+j;
          for (i=k; i<=n-1; i++)
            { v=(k-1)*n+i;
              d[u]=d[u]-a[v]*d[i*m+j];
            }
          v=(k-1)*n+k-1;
          d[u]=d[u]/a[v];
        }
    }
  return(2);
}
/////////////////////////////////////////////////////////////
int aggje(double a[],int n,double b[])
{ int *js,i,j,k,is,u,v;
  double d,t;
  js=(int *)malloc(n*sizeof(int));
  for (k=0; k<=n-1; k++)
    { d=0.0;
      for (i=k; i<=n-1; i++)
      for (j=k; j<=n-1; j++)
        { t=fabs(a[i*n+j]);
          if (t>d) {d=t; js[k]=j; is=i;}
        }
      if (d+1.0==1.0)
        { free(js); printf("fail\n"); return(0);}
      if (is!=k)
        { for (j=k; j<=n-1; j++)
            { u=k*n+j; v=is*n+j;
              t=a[u]; a[u]=a[v]; a[v]=t;
            }
          t=b[k]; b[k]=b[is]; b[is]=t;
        }
      if (js[k]!=k)
        for (i=0; i<=n-1; i++)
          { u=i*n+k; v=i*n+js[k];
            t=a[u]; a[u]=a[v]; a[v]=t;
          }
      t=a[k*n+k];
      for (j=k+1; j<=n-1; j++)
        { u=k*n+j;
          if (a[u]!=0.0) a[u]=a[u]/t;
        }
      b[k]=b[k]/t;
      for (j=k+1; j<=n-1; j++)
        { u=k*n+j;
          if (a[u]!=0.0)
            { for (i=0; i<=n-1; i++)
                { v=i*n+k;
                  if ((i!=k)&&(a[v]!=0.0))
                    { is=i*n+j;
                      a[is]=a[is]-a[v]*a[u];
                    }
                }
            }
        }
      for (i=0; i<=n-1; i++)
        { u=i*n+k;
          if ((i!=k)&&(a[u]!=0.0))
            b[i]=b[i]-a[u]*b[k];
        }
    }
  for (k=n-1; k>=0; k--)
    if (k!=js[k])
      { t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}
  free(js);
  return(1);
}
/////////////////////////////////////////////////////////////
int atlvs(double t[],int n,double b[],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.0==1.0)
    { free(s); free(y); printf("fail\n"); return(-1);}
  y[0]=1.0; x[0]=b[0]/a;
  for (k=1; k<=n-1; k++)
    { beta=0.0; q=0.0;
      for (j=0; j<=k-1; j++)
        { beta=beta+y[j]*t[j+1];
          q=q+x[j]*t[k-j];
        }
      if (fabs(a)+1.0==1.0)
        { free(s); free(y); printf("fail\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-1; i++)
          s[i]=y[i-1]+c*y[k-i-1];
      a=a+c*beta;
      if (fabs(a)+1.0==1.0)
        { free(s); free(y); printf("fail\n"); return(-1);}
      h=(b[k]-q)/a;
      for (i=0; i<=k-1; i++)
        { x[i]=x[i]+h*s[i]; y[i]=s[i];}
      x[k]=h*y[k];
    }
  free(s); free(y);
  return(1);
}
/////////////////////////////////////////////////////////////
int agsdl(double a[],double b[],int n,double x[],double eps)
{ int i,j,u,v;
  double p,t,s,q;
  for (i=0; i<=n-1; i++)
    { u=i*n+i; p=0.0; x[i]=0.0;
      for (j=0; j<=n-1; j++)
        if (i!=j)
          { v=i*n+j; p=p+fabs(a[v]);}
      if (p>=fabs(a[u]))
        { printf("fail\n"); return(-1);}
    }
  p=eps+1.0;
  while (p>=eps)
    { p=0.0;
      for (i=0; i<=n-1; i++)
        { t=x[i]; s=0.0;
          for (j=0; j<=n-1; j++)
            if (j!=i) s=s+a[i*n+j]*x[j];
          x[i]=(b[i]-s)/a[i*n+i];
          q=fabs(x[i]-t)/(1.0+fabs(x[i]));
          if (q>p) p=q;
        }
    }
  return(1);
}
/////////////////////////////////////////////////////////////
void agrad(double a[],int n,double b[],double eps,double x[])
{ int i,k;
  double *p,*r,*s,*q,alpha,beta,d,e;
  //extern void brmul(double* pa,double* pb,int m,int n,int k,double* pc);
  p=(double*)malloc(n*sizeof(double));
  r=(double*)malloc(n*sizeof(double));
  s=(double*)malloc(n*sizeof(double));
  q=(double*)malloc(n*sizeof(double));
  for (i=0; i<=n-1; i++)
    { x[i]=0.0; p[i]=b[i]; r[i]=b[i]; }
  i=0;
  while(i<=n-1)
  { 
	  brmul(a,p,n,n,1,s);
      d=0.0; e=0.0;
      for (k=0; k<=n-1; k++)
         { d=d+p[k]*b[k]; e=e+p[k]*s[k]; }
      alpha=d/e;
      for (k=0; k<=n-1; k++)
         x[k]=x[k]+alpha*p[k];
      brmul(a,x,n,n,1,q);
      d=0.0;
      for (k=0; k<=n-1; k++)
        { r[k]=b[k]-q[k]; d=d+r[k]*s[k]; }
      beta=d/e; d=0.0;
      for (k=0; k<=n-1; k++) d=d+r[k]*r[k];
      d=sqrt(d);
      if (d<eps) 
        { free(p); free(r); free(s); free(q);return;}
      for (k=0; k<=n-1; k++)
         p[k]=r[k]-beta*p[k];
      i=i+1;
  }
  free(p); free(r); free(s); free(q);
  return;
}
/////////////////////////////////////////////////////////////
int agmqr(double a[],int m,int n,double b[],double q[])
{ int i,j;
  double d,*c;
  //extern int bmaqr();
  c=(double*)malloc(n*sizeof(double));
  i=bmaqr(a,m,n,q);
  if (i==0) { free(c); return(0);}
  for (i=0; i<=n-1; i++)
    { d=0.0;
      for (j=0; j<=m-1; j++)
        d=d+q[j*m+i]*b[j];
      c[i]=d;
    }
  b[n-1]=c[n-1]/a[n*n-1];
  for (i=n-2; i>=0; i--)
    { d=0.0;
      for (j=i+1; j<=n-1; j++)
        d=d+a[i*n+j]*b[j];
      b[i]=(c[i]-d)/a[i*n+i];
    }
  free(c); return(1);
}
/////////////////////////////////////////////////////////////

int agmiv(double a[],int m,int n,double b[],double x[],double aa[],double eps,double u[],double v[],int ka)
{ int i,j;
  //extern int bginv();
  i=bginv(a,m,n,aa,eps,u,v,ka);
  if (i<0) return(-1);
  for (i=0; i<=n-1; i++)
    { x[i]=0.0;
      for (j=0; j<=m-1; j++)
        x[i]=x[i]+aa[i*m+j]*b[j];
    }
  return(1);
}
/////////////////////////////////////////////////////////////
int abint(double a[],int n,double b[],double eps,double x[])
{ int i,j,k,kk;
  double *p,*r,*e,q,qq;
  //extern int agaus();
  //extern void brmul();
  p=(double*)malloc(n*n*sizeof(double));
  r=(double*)malloc(n*sizeof(double));
  e=(double*)malloc(n*sizeof(double));
  i=60;
  for (k=0; k<=n-1; k++)
  for (j=0; j<=n-1; j++)
    p[k*n+j]=a[k*n+j];
  for (k=0; k<=n-1; k++) x[k]=b[k];
  kk=agaus(p,x,n);
  if (kk==0)
    { free(p); free(r); free(e); return(kk); }
  q=1.0+eps;
  while (q>=eps)
    { if (i==0)
        { free(p); free(r); free(e); return(i); }
      i=i-1;
      brmul(a,x,n,n,1,e);
      for ( k=0; k<=n-1; k++) r[k]=b[k]-e[k];
      for ( k=0; k<=n-1; k++)
      for ( j=0; j<=n-1; j++)
         p[k*n+j]=a[k*n+j];
      kk=agaus(p,r,n);
      if (kk==0)
        { free(p); free(r); free(e); return(kk); }
      q=0.0;
      for ( k=0; k<=n-1; k++)
        { qq=fabs(r[k])/(1.0+fabs(x[k]+r[k]));
          if (qq>q) q=qq;
        }
      for ( k=0; k<=n-1; k++) x[k]=x[k]+r[k];
    }
  free(p); free(r); free(e); return(1);
}
/////////////////////////////////////////////////////////////
#endif