数值计算程序大放送线性代数方程组

//复系数方程组的全选主元高斯消去法
//ar-n*n存放复系数矩阵的实部&＃xff0c;要被破坏
//ai-n*n存放复系数矩阵的虚部&＃xff0c;要被破坏
//n-方程组的阶数
//br-存放方程组右端复常数向量的实部&＃xff0c;返回时存放解向量的实部
//bi-存放方程组右端复常数向量的虚部&＃xff0c;返回时存放解向量的虚部
//返回值0表示系数矩阵奇异
int cbcgaus(double ar[],double ai[],int n,double br[],double bi[])
{
int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
    js&＃61;malloc(n*sizeof(int));
    for (k&＃61;0;k<&＃61;n-2;k&＃43;&＃43;)
{
  d&＃61;0.0;
        for (i&＃61;k;i<&＃61;n-1;i&＃43;&＃43;)
        {
   for (j&＃61;k;j<&＃61;n-1;j&＃43;&＃43;)
   {
    u&＃61;i*n&＃43;j;
    p&＃61;ar[u]*ar[u]&＃43;ai[u]*ai[u];
    if (p>d)
    {
     d&＃61;p;
     js[k]&＃61;j;
     is&＃61;i;
    }
   }
        }
        if (d&＃43;1.0&＃61;&＃61;1.0)
  {
   free(js);
   printf("err**fail/n");
            return(0);
  }
        if (is!&＃61;k)
  {
   for (j&＃61;k;j<&＃61;n-1;j&＃43;&＃43;)
   {
    u&＃61;k*n&＃43;j;
    v&＃61;is*n&＃43;j;
                p&＃61;ar[u];
                ar[u]&＃61;ar[v];
                ar[v]&＃61;p;
                p&＃61;ai[u];
                ai[u]&＃61;ai[v];
                ai[v]&＃61;p;
   }
            p&＃61;br[k];
            br[k]&＃61;br[is];
            br[is]&＃61;p;
            p&＃61;bi[k];
            bi[k]&＃61;bi[is];
            bi[is]&＃61;p;
  }
        if (js[k]!&＃61;k)
        {
   for (i&＃61;0;i<&＃61;n-1;i&＃43;&＃43;)
            {
    u&＃61;i*n&＃43;k;
    v&＃61;i*n&＃43;js[k];
    p&＃61;ar[u];
    ar[u]&＃61;ar[v];
    ar[v]&＃61;p;
    p&＃61;ai[u];
    ai[u]&＃61;ai[v];
    ai[v]&＃61;p;
            }
        }
        v&＃61;k*n&＃43;k;
        for (j&＃61;k&＃43;1;j<&＃61;n-1;j&＃43;&＃43;)
  {
   u&＃61;k*n&＃43;j;
            p&＃61;ar[u]*ar[v];
   q&＃61;-ai[u]*ai[v];
            s&＃61;(ar[v]-ai[v])*(ar[u]&＃43;ai[u]);
            ar[u]&＃61;(p-q)/d;
   ai[u]&＃61;(s-p-q)/d;
  }
        p&＃61;br[k]*ar[v];
  q&＃61;-bi[k]*ai[v];
        s&＃61;(ar[v]-ai[v])*(br[k]&＃43;bi[k]);
        br[k]&＃61;(p-q)/d;
  bi[k]&＃61;(s-p-q)/d;
        for (i&＃61;k&＃43;1;i<&＃61;n-1;i&＃43;&＃43;)
  {
   u&＃61;i*n&＃43;k;
            for (j&＃61;k&＃43;1;j<&＃61;n-1;j&＃43;&＃43;)
   {
    v&＃61;k*n&＃43;j;
    l&＃61;i*n&＃43;j;
                p&＃61;ar[u]*ar[v];
    q&＃61;ai[u]*ai[v];
                s&＃61;(ar[u]&＃43;ai[u])*(ar[v]&＃43;ai[v]);
                ar[l]&＃61;ar[l]-p&＃43;q;
                ai[l]&＃61;ai[l]-s&＃43;p&＃43;q;
   }
            p&＃61;ar[u]*br[k];
   q&＃61;ai[u]*bi[k];
            s&＃61;(ar[u]&＃43;ai[u])*(br[k]&＃43;bi[k]);
            br[i]&＃61;br[i]-p&＃43;q;
   bi[i]&＃61;bi[i]-s&＃43;p&＃43;q;
  }
}
    u&＃61;(n-1)*n&＃43;n-1;
    d&＃61;ar[u]*ar[u]&＃43;ai[u]*ai[u];
    if (d&＃43;1.0&＃61;&＃61;1.0)
{
  free(js);
  printf("err**fail/n");
        return(0);
}
    p&＃61;ar[u]*br[n-1];
q&＃61;-ai[u]*bi[n-1];
    s&＃61;(ar[u]-ai[u])*(br[n-1]&＃43;bi[n-1]);
    br[n-1]&＃61;(p-q)/d;
bi[n-1]&＃61;(s-p-q)/d;
    for (i&＃61;n-2;i>&＃61;0;i--)
    {
  for (j&＃61;i&＃43;1;j<&＃61;n-1;j&＃43;&＃43;)
  {
   u&＃61;i*n&＃43;j;
   p&＃61;ar[u]*br[j];
   q&＃61;ai[u]*bi[j];
   s&＃61;(ar[u]&＃43;ai[u])*(br[j]&＃43;bi[j]);
   br[i]&＃61;br[i]-p&＃43;q;
   bi[i]&＃61;bi[i]-s&＃43;p&＃43;q;
  }
}
    js[n-1]&＃61;n-1;
    for (k&＃61;n-1;k>&＃61;0;k--)
{
  if (js[k]!&＃61;k)
        {
   p&＃61;br[k]; br[k]&＃61;br[js[k]]; br[js[k]]&＃61;p;
   p&＃61;bi[k]; bi[k]&＃61;bi[js[k]]; bi[js[k]]&＃61;p;
        }
}
free(js);
return(1);
}

//全选主元高斯-约当消去法
//a-n*n 存放方程组的系数矩阵&＃xff0c;返回时将被破坏
//b-n*m常数向量,每列为一组&＃xff0c;返回时存放m组解向量
//n-方程组的阶数
//m-方程组右端常数向量的个数
//返回0表示原方程组的系数矩阵奇异
int ccgj(double a[],double b[],int n,int m)
{
int *js,l,k,i,j,is,p,q;
    double d,t;
    js&＃61;malloc(n*sizeof(int));
    l&＃61;1;
    for (k&＃61;0;k<&＃61;n-1;k&＃43;&＃43;)
{
  d&＃61;0.0;
        for (i&＃61;k;i<&＃61;n-1;i&＃43;&＃43;)
  {
   for (j&＃61;k;j<&＃61;n-1;j&＃43;&＃43;)
            {
    t&＃61;fabs(a[i*n&＃43;j]);
    if (t>d)
    {
     d&＃61;t;js[k]&＃61;j;is&＃61;i;
    }
            }
  }
        if (d&＃43;1.0&＃61;&＃61;1.0)
  {
   l&＃61;0;
  }
        else
  {
   if (js[k]!&＃61;k)
   {
    for (i&＃61;0;i<&＃61;n-1;i&＃43;&＃43;)
                {
     p&＃61;i*n&＃43;k; q&＃61;i*n&＃43;js[k]; t&＃61;a[p]; a[p]&＃61;a[q]; a[q]&＃61;t;
                }
   }
   if (is!&＃61;k)
   {
    for (j&＃61;k;j<&＃61;n-1;j&＃43;&＃43;)
    {
     p&＃61;k*n&＃43;j; q&＃61;is*n&＃43;j; t&＃61;a[p]; a[p]&＃61;a[q]; a[q]&＃61;t;
    }
    for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
    {
     p&＃61;k*m&＃43;j; q&＃61;is*m&＃43;j; t&＃61;b[p]; b[p]&＃61;b[q]; b[q]&＃61;t;
    }
   }
  }
        if (l&＃61;&＃61;0)
  {
   free(js); printf("fail/n");
            return(0);
  }
        d&＃61;a[k*n&＃43;k];
        for (j&＃61;k&＃43;1;j<&＃61;n-1;j&＃43;&＃43;)
  {
   p&＃61;k*n&＃43;j; a[p]&＃61;a[p]/d;
  }
        for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
  {
   p&＃61;k*m&＃43;j; b[p]&＃61;b[p]/d;
  }
        for (j&＃61;k&＃43;1;j<&＃61;n-1;j&＃43;&＃43;)
  {
   for (i&＃61;0;i<&＃61;n-1;i&＃43;&＃43;)
            {
    p&＃61;i*n&＃43;j;
    if (i!&＃61;k)
    {
     a[p]&＃61;a[p]-a[i*n&＃43;k]*a[k*n&＃43;j];
    }
            }
  }
        for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
  {
   for (i&＃61;0;i<&＃61;n-1;i&＃43;&＃43;)
   {
    p&＃61;i*m&＃43;j;
    if (i!&＃61;k)
    {
     b[p]&＃61;b[p]-a[i*n&＃43;k]*b[k*m&＃43;j];
    }
   }
  }
}
    for (k&＃61;n-1;k>&＃61;0;k--)
{
  if (js[k]!&＃61;k)
  {
   for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
   {
    p&＃61;k*m&＃43;j; q&＃61;js[k]*m&＃43;j; t&＃61;b[p]; b[p]&＃61;b[q]; b[q]&＃61;t;
   }
  }
}
    free(js);
    return(1);
}

//复系数全选主元高斯-约当消去法
//ar-n*n 存放方程组的系数矩阵的实部&＃xff0c;返回时将被破坏
//ai-n*n 存放方程组的系数矩阵的虚部&＃xff0c;返回时将被破坏
//br-n*m常数向量的实部,每列为一组&＃xff0c;返回时存放m组解向量的实部
//bi-n*m常数向量的虚部,每列为一组&＃xff0c;返回时存放m组解向量的虚部
//n-方程组的阶数
//m-方程组右端常数向量的个数
//返回0表示原方程组的系数矩阵奇异
int cdcgj(double ar[],double ai[],double br[],double bi[],int n,int m)
{
int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
    js&＃61;malloc(n*sizeof(int));
    for (k&＃61;0;k<&＃61;n-1;k&＃43;&＃43;)
{
  d&＃61;0.0;
        for (i&＃61;k;i<&＃61;n-1;i&＃43;&＃43;)
  {
   for (j&＃61;k;j<&＃61;n-1;j&＃43;&＃43;)
   {
    u&＃61;i*n&＃43;j;
    p&＃61;ar[u]*ar[u]&＃43;ai[u]*ai[u];
    if (p>d)
    {
     d&＃61;p;js[k]&＃61;j;is&＃61;i;
    }
   }
  }
        if (d&＃43;1.0&＃61;&＃61;1.0)
  {
   free(js);
   printf("err**fail/n");
            return(0);
  }
        if (is!&＃61;k)
  {
   for (j&＃61;k;j<&＃61;n-1;j&＃43;&＃43;)
   {
    u&＃61;k*n&＃43;j; v&＃61;is*n&＃43;j;
                p&＃61;ar[u]; ar[u]&＃61;ar[v]; ar[v]&＃61;p;
                p&＃61;ai[u]; ai[u]&＃61;ai[v]; ai[v]&＃61;p;
   }
            for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
   {
    u&＃61;k*m&＃43;j; v&＃61;is*m&＃43;j;
                p&＃61;br[u]; br[u]&＃61;br[v]; br[v]&＃61;p;
                p&＃61;bi[u]; bi[u]&＃61;bi[v]; bi[v]&＃61;p;
   }
  }
        if (js[k]!&＃61;k)
  {
   for (i&＃61;0;i<&＃61;n-1;i&＃43;&＃43;)
            {
    u&＃61;i*n&＃43;k; v&＃61;i*n&＃43;js[k];
    p&＃61;ar[u]; ar[u]&＃61;ar[v]; ar[v]&＃61;p;
    p&＃61;ai[u]; ai[u]&＃61;ai[v]; ai[v]&＃61;p;
            }
  }
        v&＃61;k*n&＃43;k;
        for (j&＃61;k&＃43;1;j<&＃61;n-1;j&＃43;&＃43;)
  {
   u&＃61;k*n&＃43;j;
            p&＃61;ar[u]*ar[v]; q&＃61;-ai[u]*ai[v];
            s&＃61;(ar[v]-ai[v])*(ar[u]&＃43;ai[u]);
            ar[u]&＃61;(p-q)/d; ai[u]&＃61;(s-p-q)/d;
  }
        for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
  {
   u&＃61;k*m&＃43;j;
            p&＃61;br[u]*ar[v]; q&＃61;-bi[u]*ai[v];
            s&＃61;(ar[v]-ai[v])*(br[u]&＃43;bi[u]);
            br[u]&＃61;(p-q)/d; bi[u]&＃61;(s-p-q)/d;
  }
        for (i&＃61;0;i<&＃61;n-1;i&＃43;&＃43;)
  {
   if (i!&＃61;k)
            {
    u&＃61;i*n&＃43;k;
    for (j&＃61;k&＃43;1;j<&＃61;n-1;j&＃43;&＃43;)
                {
     v&＃61;k*n&＃43;j; l&＃61;i*n&＃43;j;
     p&＃61;ar[u]*ar[v]; q&＃61;ai[u]*ai[v];
     s&＃61;(ar[u]&＃43;ai[u])*(ar[v]&＃43;ai[v]);
     ar[l]&＃61;ar[l]-p&＃43;q;
     ai[l]&＃61;ai[l]-s&＃43;p&＃43;q;
                }
    for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
                {
     l&＃61;i*m&＃43;j; v&＃61;k*m&＃43;j;
     p&＃61;ar[u]*br[v]; q&＃61;ai[u]*bi[v];
     s&＃61;(ar[u]&＃43;ai[u])*(br[v]&＃43;bi[v]);
     br[l]&＃61;br[l]-p&＃43;q; bi[l]&＃61;bi[l]-s&＃43;p&＃43;q;
                }
            }
  }
}
    for (k&＃61;n-1;k>&＃61;0;k--)
{
  if (js[k]!&＃61;k)
  {
   for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
   {
    u&＃61;k*m&＃43;j;v&＃61;js[k]*m&＃43;j;
    p&＃61;br[u]; br[u]&＃61;br[v]; br[v]&＃61;p;
    p&＃61;bi[u]; bi[u]&＃61;bi[v]; bi[v]&＃61;p;
   }
  }
}
    free(js);
    return(1);
}

//求解三对角线方程组的追赶法
//存放三对角线矩阵中三对角线上的元素
//存放顺序为a00,a01,a10,a11,a12,a21,a22,a23,...,an-1,n-2,an-1,n-1
//就是按行依次存储
//n-方程组的阶数
//m-b的长度&＃xff0c;应为m&＃61;3n-2,本函数中用来做检验
//d-存放方程组右端常数向量&＃xff0c;返回时存放解向量
int cetrd(double b[],int n,int m,double d[])
{
int k,j;
    double s;
    if (m!&＃61;(3*n-2))
{
  printf("err/n"); return(-2);
}
    for (k&＃61;0;k<&＃61;n-2;k&＃43;&＃43;)
{
  j&＃61;3*k; s&＃61;b[j];
        if (fabs(s)&＃43;1.0&＃61;&＃61;1.0)
  {
   printf("fail/n"); return(0);
  }
        b[j&＃43;1]&＃61;b[j&＃43;1]/s;
        d[k]&＃61;d[k]/s;
        b[j&＃43;3]&＃61;b[j&＃43;3]-b[j&＃43;2]*b[j&＃43;1];
        d[k&＃43;1]&＃61;d[k&＃43;1]-b[j&＃43;2]*d[k];
}
    s&＃61;b[3*n-3];
    if (fabs(s)&＃43;1.0&＃61;&＃61;1.0)
{
  printf("fail/n"); return(0);
}
    d[n-1]&＃61;d[n-1]/s;
    for (k&＃61;n-2;k>&＃61;0;k--)
{
  d[k]&＃61;d[k]-b[3*k&＃43;1]*d[k&＃43;1];
}
    return(2);
}

//一般带型方程组的求解
//n-方程组的阶数
//l-系数矩阵的半带宽
//il-系数矩阵的带宽&＃xff0c;应为il&＃61;2l&＃43;1
//m-方程组右端的常数
int cfband(double b[],double d[],int n,int l,int il,int m)
{
int ls,k,i,j,is,u,v;
    double p,t;
    if (il!&＃61;(2*l&＃43;1))
{
  printf("fail/n"); return(-2);
}
    ls&＃61;l;
    for (k&＃61;0;k<&＃61;n-2;k&＃43;&＃43;)
{
  p&＃61;0.0;
        for (i&＃61;k;i<&＃61;ls;i&＃43;&＃43;)
  {
   t&＃61;fabs(b[i*il]);
            if (t>p)
   {
    p&＃61;t; is&＃61;i;
   }
  }
        if (p&＃43;1.0&＃61;&＃61;1.0)
  {
   printf("fail/n"); return(0);
  }
        for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
  {
   u&＃61;k*m&＃43;j; v&＃61;is*m&＃43;j;
            t&＃61;d[u]; d[u]&＃61;d[v]; d[v]&＃61;t;
  }
        for (j&＃61;0;j<&＃61;il-1;j&＃43;&＃43;)
  {
   u&＃61;k*il&＃43;j; v&＃61;is*il&＃43;j;
            t&＃61;b[u]; b[u]&＃61;b[v]; b[v]&＃61;t;
  }
        for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
  {
   u&＃61;k*m&＃43;j; d[u]&＃61;d[u]/b[k*il];
  }
        for (j&＃61;1;j<&＃61;il-1;j&＃43;&＃43;)
  {
   u&＃61;k*il&＃43;j; b[u]&＃61;b[u]/b[k*il];
  }
        for (i&＃61;k&＃43;1;i<&＃61;ls;i&＃43;&＃43;)
  {
   t&＃61;b[i*il];
            for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
   {
    u&＃61;i*m&＃43;j; v&＃61;k*m&＃43;j;
                d[u]&＃61;d[u]-t*d[v];
   }
            for (j&＃61;1;j<&＃61;il-1;j&＃43;&＃43;)
   {
    u&＃61;i*il&＃43;j; v&＃61;k*il&＃43;j;
                b[u-1]&＃61;b[u]-t*b[v];
   }
            u&＃61;i*il&＃43;il-1; b[u]&＃61;0.0;
  }
        if (ls!&＃61;(n-1))
  {
   ls&＃61;ls&＃43;1;
  }
}
    p&＃61;b[(n-1)*il];
    if (fabs(p)&＃43;1.0&＃61;&＃61;1.0)
{
  printf("fail/n"); return(0);
}
    for (j&＃61;0;j<&＃61;m-1;j&＃43;&＃43;)
{
  u&＃61;(n-1)*m&＃43;j; d[u]&＃61;d[u]/p;
}
    ls&＃61;1;
    for (i&＃61;n-2;i>&＃61;0;i--)
{
  for (k&＃61;0;k<&＃61;m-1;k&＃43;&＃43;)
  {
   u&＃61;i*m&＃43;k;
            for (j&＃61;1;j<&＃61;ls;j&＃43;&＃43;)
   {
    v&＃61;i*il&＃43;j; is&＃61;(i&＃43;j)*m&＃43;k;
                d[u]&＃61;d[u]-b[v]*d[is];
   }
  }
        if (ls!&＃61;(il-1)) ls&＃61;ls&＃43;1;
}
    return(2);
}

//求解对称方程组的分解法
//a n*n数组&＃xff0c;存放方程组的系数矩阵&＃xff08;应为对称矩阵&＃xff09;
//n 方程组的阶数
//m 方程组右端常数向量的个数
//c n*m 存放方程组右端的m组常数向量&＃xff1b;返回时存放方程组的m组解向量
//返回值小于0表示程序工作失败&＃xff0c;大于0表示正常返回
int cgldl(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])&＃43;1.0&＃61;&＃61;1.0)
{
  printf("fail/n");
  return(-2);
}
for (i&＃61;1; i<&＃61;n-1; i&＃43;&＃43;)
{
  u&＃61;i*n;
  a[u]&＃61;a[u]/a[0];
}
for (i&＃61;1; i<&＃61;n-2; i&＃43;&＃43;)
{
  u&＃61;i*n&＃43;i;
  for (j&＃61;1; j<&＃61;i; j&＃43;&＃43;)
  {
   v&＃61;i*n&＃43;j-1;
   l&＃61;(j-1)*n&＃43;j-1;
   a[u]&＃61;a[u]-a[v]*a[v]*a[l];
  }
  p&＃61;a[u];
  if (fabs(p)&＃43;1.0&＃61;&＃61;1.0)
  {
   printf("fail/n");
   return(-2);
  }
  for (k&＃61;i&＃43;1; k<&＃61;n-1; k&＃43;&＃43;)
  {
   u&＃61;k*n&＃43;i;
   for (j&＃61;1; j<&＃61;i; j&＃43;&＃43;)
   {
    v&＃61;k*n&＃43;j-1;
    l&＃61;i*n&＃43;j-1;
    w&＃61;(j-1)*n&＃43;j-1;
    a[u]&＃61;a[u]-a[v]*a[l]*a[w];
   }
   a[u]&＃61;a[u]/p;
  }
}
u&＃61;n*n-1;
for (j&＃61;1; j<&＃61;n-1; j&＃43;&＃43;)
{
  v&＃61;(n-1)*n&＃43;j-1;
  w&＃61;(j-1)*n&＃43;j-1;
  a[u]&＃61;a[u]-a[v]*a[v]*a[w];
}
p&＃61;a[u];
if (fabs(p)&＃43;1.0&＃61;&＃61;1.0)
{
  printf("fail/n");
  return(-2);
}
for (j&＃61;0; j<&＃61;m-1; j&＃43;&＃43;)
{
  for (i&＃61;1; i<&＃61;n-1; i&＃43;&＃43;)
  {
   u&＃61;i*m&＃43;j;
   for (k&＃61;1; k<&＃61;i; k&＃43;&＃43;)
   {
    v&＃61;i*n&＃43;k-1;
    w&＃61;(k-1)*m&＃43;j;
    c[u]&＃61;c[u]-a[v]*c[w];
   }
  }
}
for (i&＃61;1; i<&＃61;n-1; i&＃43;&＃43;)
{
  u&＃61;(i-1)*n&＃43;i-1;
  for (j&＃61;i; j<&＃61;n-1; j&＃43;&＃43;)
  {
   v&＃61;(i-1)*n&＃43;j;
   w&＃61;j*n&＃43;i-1;
   a[v]&＃61;a[u]*a[w];
  }
}
for (j&＃61;0; j<&＃61;m-1; j&＃43;&＃43;)
{
  u&＃61;(n-1)*m&＃43;j;
  c[u]&＃61;c[u]/p;
  for (k&＃61;1; k<&＃61;n-1; k&＃43;&＃43;)
  {
   k1&＃61;n-k;
   k3&＃61;k1-1;
   u&＃61;k3*m&＃43;j;
   for (k2&＃61;k1; k2<&＃61;n-1; k2&＃43;&＃43;)
   {
    v&＃61;k3*n&＃43;k2;
    w&＃61;k2*m&＃43;j;
    c[u]&＃61;c[u]-a[v]*c[w];
   }
   c[u]&＃61;c[u]/a[k3*n&＃43;k3];
  }
}
return(2);
}

//求解对称正定方程组的平方根法
//用Cholesky分解法&＃xff08;即平方根法&＃xff09;
//a n*n 存放系数矩阵&＃xff08;应为对称正定矩阵&＃xff09;&＃xff0c;返回时其上三角部分存放分解后的U矩阵
//n 方程的阶数
//m 方程组有短的常数向量的个数
//d n*m 存放方程组右端m组常数向量&＃xff1b;返回时&＃xff0c;存放方程组的m组解向量
//返回值小于0&＃xff0c;表示程序工作失败
int chchol(double a[],int n,int m,double d[])
{
int i,j,k,u,v;
    if ((a[0]&＃43;1.0&＃61;&＃61;1.0)||(a[0]<0.0))
{
  printf("fail/n"); return(-2);
}
    a[0]&＃61;sqrt(a[0]);
    for (j&＃61;1; j<&＃61;n-1; j&＃43;&＃43;)
{
  a[j]&＃61;a[j]/a[0];
}
    for (i&＃61;1; i<&＃61;n-1; i&＃43;&＃43;)
{
  u&＃61;i*n&＃43;i;
        for (j&＃61;1; j<&＃61;i; j&＃43;&＃43;)
  {
   v&＃61;(j-1)*n&＃43;i;
            a[u]&＃61;a[u]-a[v]*a[v];
  }
        if ((a[u]&＃43;1.0&＃61;&＃61;1.0)||(a[u]<0.0))
  {
   printf("fail/n"); return(-2);
  }
        a[u]&＃61;sqrt(a[u]);
        if (i!&＃61;(n-1))
  {
   for (j&＃61;i&＃43;1; j<&＃61;n-1; j&＃43;&＃43;)
   {
    v&＃61;i*n&＃43;j;
                for (k&＃61;1; k<&＃61;i; k&＃43;&＃43;)
    {
     a[v]&＃61;a[v]-a[(k-1)*n&＃43;i]*a[(k-1)*n&＃43;j];
    }
                a[v]&＃61;a[v]/a[u];
   }
  }
}
    for (j&＃61;0; j<&＃61;m-1; j&＃43;&＃43;)
{
  d[j]&＃61;d[j]/a[0];
        for (i&＃61;1; i<&＃61;n-1; i&＃43;&＃43;)
  {
   u&＃61;i*n&＃43;i; v&＃61;i*m&＃43;j;
            for (k&＃61;1; k<&＃61;i; k&＃43;&＃43;)
   {
    d[v]&＃61;d[v]-a[(k-1)*n&＃43;i]*d[(k-1)*m&＃43;j];
   }
            d[v]&＃61;d[v]/a[u];
  }
}
    for (j&＃61;0; j<&＃61;m-1; j&＃43;&＃43;)
{
  u&＃61;(n-1)*m&＃43;j;
        d[u]&＃61;d[u]/a[n*n-1];
        for (k&＃61;n-1; k>&＃61;1; k--)
  {
   u&＃61;(k-1)*m&＃43;j;
            for (i&＃61;k; i<&＃61;n-1; i&＃43;&＃43;)
   {
    v&＃61;(k-1)*n&＃43;i;
                d[u]&＃61;d[u]-a[v]*d[i*m&＃43;j];
   }
            v&＃61;(k-1)*n&＃43;k-1;
            d[u]&＃61;d[u]/a[v];
  }
}
    return(2);
}

//求解大型稀疏方程组的全选主元高斯-约当消去法
//只是在消去过程中避免了对零元素的运算&＃xff0c;从而可以大大减少运算次数
//本函数不考虑系数矩阵的压缩存储问题
//a n*n数组&＃xff0c;存放方程组的系数矩阵&＃xff0c;返回时要被破坏
//n 方程组的阶数
//d-存放方程组右端常数向量&＃xff0c;返回时存放解向量
//返回值为0表示系数矩阵奇异
int ciggj(double a[],int n,double b[])
{
int *js,i,j,k,is,u,v;
    double d,t;
    js&＃61;malloc(n*sizeof(int));
    for (k&＃61;0; k<&＃61;n-1; k&＃43;&＃43;)
{
  d&＃61;0.0;
        for (i&＃61;k; i<&＃61;n-1; i&＃43;&＃43;)
        {
   for (j&＃61;k; j<&＃61;n-1; j&＃43;&＃43;)
   {
    t&＃61;fabs(a[i*n&＃43;j]);
    if (t>d)
    {
     d&＃61;t; js[k]&＃61;j; is&＃61;i;
    }
   }
        }
        if (d&＃43;1.0&＃61;&＃61;1.0)
  {
   free(js); printf("fail/n"); return(0);
  }
        if (is!&＃61;k)
  {
   for (j&＃61;k; j<&＃61;n-1; j&＃43;&＃43;)
   {
    u&＃61;k*n&＃43;j; v&＃61;is*n&＃43;j;
                t&＃61;a[u]; a[u]&＃61;a[v]; a[v]&＃61;t;
   }
            t&＃61;b[k]; b[k]&＃61;b[is]; b[is]&＃61;t;
  }
        if (js[k]!&＃61;k)
        {
   for (i&＃61;0; i<&＃61;n-1; i&＃43;&＃43;)
            {
    u&＃61;i*n&＃43;k; v&＃61;i*n&＃43;js[k];
    t&＃61;a[u]; a[u]&＃61;a[v]; a[v]&＃61;t;
            }
        }
        t&＃61;a[k*n&＃43;k];
        for (j&＃61;k&＃43;1; j<&＃61;n-1; j&＃43;&＃43;)
  {
   u&＃61;k*n&＃43;j;
            if (a[u]!&＃61;0.0)
            {
    a[u]&＃61;a[u]/t;
   }
  }
        b[k]&＃61;b[k]/t;
        for (j&＃61;k&＃43;1; j<&＃61;n-1; j&＃43;&＃43;)
  {
   u&＃61;k*n&＃43;j;
            if (a[u]!&＃61;0.0)
   {
    for (i&＃61;0; i<&＃61;n-1; i&＃43;&＃43;)
    {
     v&＃61;i*n&＃43;k;
                    if ((i!&＃61;k)&&(a[v]!&＃61;0.0))
     {
      is&＃61;i*n&＃43;j;
                        a[is]&＃61;a[is]-a[v]*a[u];
     }
    }
   }
  }
        for (i&＃61;0; i<&＃61;n-1; i&＃43;&＃43;)
  {
   u&＃61;i*n&＃43;k;
            if ((i!&＃61;k)&&(a[u]!&＃61;0.0))
            {
    b[i]&＃61;b[i]-a[u]*b[k];
   }
  }
}
    for (k&＃61;n-1; k>&＃61;0; k--)
    {
  if (k!&＃61;js[k])
        {
   t&＃61;b[k]; b[k]&＃61;b[js[k]]; b[js[k]]&＃61;t;
  }
}
    free(js);
    return(1);
}

//求解Toeplitz型方程组的Levinson递推算法
//t 长度为n 的一维数组&＃xff0c;存放对称T型矩阵中的元素t0,t1,...tn-1
//n 方程组的阶数
//b 长度为n 的一维数组,存放方程组右端的常数向量
//x 长度为n 的一维数组,返回方程组的解
//返回值小于0表示程序工作失败
int cjtlvs(double t[],int n,double b[],double x[])
{
int i,j,k;
    double a,beta,q,c,h,*y,*s;
    s&＃61;malloc(n*sizeof(double));
    y&＃61;malloc(n*sizeof(double));
    a&＃61;t[0];
    if (fabs(a)&＃43;1.0&＃61;&＃61;1.0)
{
  free(s); printf("fail/n"); return(-1);
}
    y[0]&＃61;1.0; x[0]&＃61;b[0]/a;
    for (k&＃61;1; k<&＃61;n-1; k&＃43;&＃43;)
{
  beta&＃61;0.0; q&＃61;0.0;
        for (j&＃61;0; j<&＃61;k-1; j&＃43;&＃43;)
  {
   beta&＃61;beta&＃43;y[j]*t[j&＃43;1];
            q&＃61;q&＃43;x[j]*t[k-j];
  }
        if (fabs(a)&＃43;1.0&＃61;&＃61;1.0)
  {
   free(s); printf("fail/n"); return(-1);
  }
        c&＃61;-beta/a; s[0]&＃61;c*y[k-1]; y[k]&＃61;y[k-1];
        if (k!&＃61;1)
  {
   for (i&＃61;1; i<&＃61;k-1; i&＃43;&＃43;)
   {
    s[i]&＃61;y[i-1]&＃43;c*y[k-i-1];
   }
  }
        a&＃61;a&＃43;c*beta;
        if (fabs(a)&＃43;1.0&＃61;&＃61;1.0)
  {
   free(s); printf("fail/n"); return(-1);
  }
        h&＃61;(b[k]-q)/a;
        for (i&＃61;0; i<&＃61;k-1; i&＃43;&＃43;)
  {
   x[i]&＃61;x[i]&＃43;h*s[i]; y[i]&＃61;s[i];
  }
        x[k]&＃61;h*y[k];
}
    free(s); free(y);
    return(1);
}