数值计算程序大放送-线性代数方程组
选自<<徐世良数值计算程序集(C)>>
每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。
今天都给贴出来了
#include "stdlib.h"
#include "math.h"
#include "stdio.h"
// 全选主元高斯消去法
//a-n*n 存放方程组的系数矩阵，返回时将被破坏
//b-常数向量
//x-返回方程组的解向量
//n-存放方程组的阶数
//返回0表示原方程组的系数矩阵奇异
int cagaus(double a[],double b[],int n,double x[])
{
int *js,l,k,i,j,is,p,q;
double d,t;
js=malloc(n*sizeof(int));
l=1;
for (k=0;k<=n-2;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)
{
l=0;
}
else
{
if (js[k]!=k)
{
for (i=0;i<=n-1;i++)
{
p=i*n+k;
q=i*n+js[k];
t=a[p];
a[p]=a[q];
a[q]=t;
}
}
if (is!=k)
{
for (j=k;j<=n-1;j++)
{
p=k*n+j;
q=is*n+j;
t=a[p];
a[p]=a[q];
a[q]=t;
}
t=b[k];
b[k]=b[is];
b[is]=t;
}
}
if (l==0)
{
free(js);
printf("fail\n");
return(0);
}
d=a[k*n+k];
for (j=k+1;j<=n-1;j++)
{
p=k*n+j;
a[p]=a[p]/d;
}
b[k]=b[k]/d;
for (i=k+1;i<=n-1;i++)
{
for (j=k+1;j<=n-1;j++)
{
p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[(n-1)*n+n-1];
if (fabs(d)+1.0==1.0)
{
free(js);
printf("fail\n");
return(0);
}
x[n-1]=b[n-1]/d;
for (i=n-2;i>=0;i--)
{
t=0.0;
for (j=i+1;j<=n-1;j++)
{
t=t+a[i*n+j]*x[j];
}
x[i]=b[i]-t;
}
js[n-1]=n-1;
for (k=n-1;k>=0;k--)
{
if (js[k]!=k)
{
t=x[k];
x[k]=x[js[k]];
x[js[k]]=t;
}
}
free(js);
//复系数方程组的全选主元高斯消去法
//ar-n*n存放复系数矩阵的实部，要被破坏
//ai-n*n存放复系数矩阵的虚部，要被破坏
//n-方程组的阶数
//br-存放方程组右端复常数向量的实部，返回时存放解向量的实部
//bi-存放方程组右端复常数向量的虚部，返回时存放解向量的虚部
//返回值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=malloc(n*sizeof(int));
for (k=0;k<=n-2;k++)
{
d=0.0;
for (i=k;i<=n-1;i++)
{
for (j=k;j<=n-1;j++)
{
u=i*n+j;
p=ar[u]*ar[u]+ai[u]*ai[u];
if (p>d)
{
d=p;
js[k]=j;
is=i;
}
}
}
if (d+1.0==1.0)
{
free(js);
printf("err**fail\n");
return(0);
}
if (is!=k)
{
for (j=k;j<=n-1;j++)
{
u=k*n+j;
v=is*n+j;
p=ar[u];
ar[u]=ar[v];
ar[v]=p;
p=ai[u];
ai[u]=ai[v];
ai[v]=p;
}
p=br[k];
br[k]=br[is];
br[is]=p;
p=bi[k];
bi[k]=bi[is];
bi[is]=p;
}
if (js[k]!=k)
{
for (i=0;i<=n-1;i++)
{
u=i*n+k;
v=i*n+js[k];
p=ar[u];
ar[u]=ar[v];
ar[v]=p;
p=ai[u];
ai[u]=ai[v];
ai[v]=p;
}
}
v=k*n+k;
for (j=k+1;j<=n-1;j++)
{
u=k*n+j;
p=ar[u]*ar[v];
q=-ai[u]*ai[v];
s=(ar[v]-ai[v])*(ar[u]+ai[u]);
ar[u]=(p-q)/d;
ai[u]=(s-p-q)/d;
}
p=br[k]*ar[v];
q=-bi[k]*ai[v];
s=(ar[v]-ai[v])*(br[k]+bi[k]);
br[k]=(p-q)/d;
bi[k]=(s-p-q)/d;
for (i=k+1;i<=n-1;i++)
{
u=i*n+k;
for (j=k+1;j<=n-1;j++)
{
v=k*n+j;
l=i*n+j;
p=ar[u]*ar[v];
q=ai[u]*ai[v];
s=(ar[u]+ai[u])*(ar[v]+ai[v]);
ar[l]=ar[l]-p+q;
ai[l]=ai[l]-s+p+q;
}
p=ar[u]*br[k];
q=ai[u]*bi[k];
s=(ar[u]+ai[u])*(br[k]+bi[k]);
br[i]=br[i]-p+q;
bi[i]=bi[i]-s+p+q;
}
}
u=(n-1)*n+n-1;
d=ar[u]*ar[u]+ai[u]*ai[u];
if (d+1.0==1.0)
{
free(js);
printf("err**fail\n");
return(0);
}
p=ar[u]*br[n-1];
q=-ai[u]*bi[n-1];
s=(ar[u]-ai[u])*(br[n-1]+bi[n-1]);
br[n-1]=(p-q)/d;
bi[n-1]=(s-p-q)/d;
for (i=n-2;i>=0;i--)
{
for (j=i+1;j<=n-1;j++)
{
u=i*n+j;
p=ar[u]*br[j];
q=ai[u]*bi[j];
s=(ar[u]+ai[u])*(br[j]+bi[j]);
br[i]=br[i]-p+q;
bi[i]=bi[i]-s+p+q;
}
}
js[n-1]=n-1;
for (k=n-1;k>=0;k--)
{
if (js[k]!=k)
{
p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
}
}
free(js);
return(1);
}
//全选主元高斯-约当消去法
//a-n*n 存放方程组的系数矩阵，返回时将被破坏
//b-n*m常数向量,每列为一组，返回时存放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=malloc(n*sizeof(int));
l=1;
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)
{
l=0;
}
else
{
if (js[k]!=k)
{
for (i=0;i<=n-1;i++)
{
p=i*n+k; q=i*n+js[k]; t=a[p]; a[p]=a[q]; a[q]=t;
}
}
if (is!=k)
{
for (j=k;j<=n-1;j++)
{
p=k*n+j; q=is*n+j; t=a[p]; a[p]=a[q]; a[q]=t;
}
for (j=0;j<=m-1;j++)
{
p=k*m+j; q=is*m+j; t=b[p]; b[p]=b[q]; b[q]=t;
}
}
}
if (l==0)
{
free(js); printf("fail\n");
return(0);
}
d=a[k*n+k];
for (j=k+1;j<=n-1;j++)
{
p=k*n+j; a[p]=a[p]/d;
}
for (j=0;j<=m-1;j++)
{
p=k*m+j; b[p]=b[p]/d;
}
for (j=k+1;j<=n-1;j++)
{
for (i=0;i<=n-1;i++)
{
p=i*n+j;
if (i!=k)
{
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
}
}
for (j=0;j<=m-1;j++)
{
for (i=0;i<=n-1;i++)
{
p=i*m+j;
if (i!=k)
{
b[p]=b[p]-a[i*n+k]*b[k*m+j];
}
}
}
}
for (k=n-1;k>=0;k--)
{
if (js[k]!=k)
{
for (j=0;j<=m-1;j++)
{
p=k*m+j; q=js[k]*m+j; t=b[p]; b[p]=b[q]; b[q]=t;
}
}
}
free(js);
return(1);
}
//复系数全选主元高斯-约当消去法
//ar-n*n 存放方程组的系数矩阵的实部，返回时将被破坏
//ai-n*n 存放方程组的系数矩阵的虚部，返回时将被破坏
//br-n*m常数向量的实部,每列为一组，返回时存放m组解向量的实部
//bi-n*m常数向量的虚部,每列为一组，返回时存放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=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++)
{
u=i*n+j;
p=ar[u]*ar[u]+ai[u]*ai[u];
if (p>d)
{
d=p;js[k]=j;is=i;
}
}
}
if (d+1.0==1.0)
{
free(js);
printf("err**fail\n");
return(0);
}
if (is!=k)
{
for (j=k;j<=n-1;j++)
{
u=k*n+j; v=is*n+j;
p=ar[u]; ar[u]=ar[v]; ar[v]=p;
p=ai[u]; ai[u]=ai[v]; ai[v]=p;
}
for (j=0;j<=m-1;j++)
{
u=k*m+j; v=is*m+j;
p=br[u]; br[u]=br[v]; br[v]=p;
p=bi[u]; bi[u]=bi[v]; bi[v]=p;
}
}
if (js[k]!=k)
{
for (i=0;i<=n-1;i++)
{
u=i*n+k; v=i*n+js[k];
p=ar[u]; ar[u]=ar[v]; ar[v]=p;
p=ai[u]; ai[u]=ai[v]; ai[v]=p;
}
}
v=k*n+k;
for (j=k+1;j<=n-1;j++)
{
u=k*n+j;
p=ar[u]*ar[v]; q=-ai[u]*ai[v];
s=(ar[v]-ai[v])*(ar[u]+ai[u]);
ar[u]=(p-q)/d; ai[u]=(s-p-q)/d;
}
for (j=0;j<=m-1;j++)
{
u=k*m+j;
p=br[u]*ar[v]; q=-bi[u]*ai[v];
s=(ar[v]-ai[v])*(br[u]+bi[u]);
br[u]=(p-q)/d; bi[u]=(s-p-q)/d;
}
for (i=0;i<=n-1;i++)
{
if (i!=k)
{
u=i*n+k;
for (j=k+1;j<=n-1;j++)
{
v=k*n+j; l=i*n+j;
p=ar[u]*ar[v]; q=ai[u]*ai[v];
s=(ar[u]+ai[u])*(ar[v]+ai[v]);
ar[l]=ar[l]-p+q;
ai[l]=ai[l]-s+p+q;
}
for (j=0;j<=m-1;j++)
{
l=i*m+j; v=k*m+j;
p=ar[u]*br[v]; q=ai[u]*bi[v];
s=(ar[u]+ai[u])*(br[v]+bi[v]);
br[l]=br[l]-p+q; bi[l]=bi[l]-s+p+q;
}
}
}
}
for (k=n-1;k>=0;k--)
{
if (js[k]!=k)
{
for (j=0;j<=m-1;j++)
{
u=k*m+j;v=js[k]*m+j;
p=br[u]; br[u]=br[v]; br[v]=p;
p=bi[u]; bi[u]=bi[v]; bi[v]=p;
}
}
}
free(js);
return(1);
}
//求解三对角线方程组的追赶法
//存放三对角线矩阵中三对角线上的元素
//存放顺序为a00,a01,a10,a11,a12,a21,a22,a23,...,an-1,n-2,an-1,n-1
//就是按行依次存储
//n-方程组的阶数
//m-b的长度，应为m=3n-2,本函数中用来做检验
//d-存放方程组右端常数向量，返回时存放解向量
int cetrd(double b[],int n,int m,double d[])
{
int k,j;
double s;
if (m!=(3*n-2))
{
printf("err\n"); return(-2);
}
for (k=0;k<=n-2;k++)
{
j=3*k; s=b[j];
if (fabs(s)+1.0==1.0)
{
printf("fail\n"); return(0);
}
b[j+1]=b[j+1]/s;
d[k]=d[k]/s;
b[j+3]=b[j+3]-b[j+2]*b[j+1];
d[k+1]=d[k+1]-b[j+2]*d[k];
}
s=b[3*n-3];
if (fabs(s)+1.0==1.0)
{
printf("fail\n"); return(0);
}
d[n-1]=d[n-1]/s;
for (k=n-2;k>=0;k--)
{
d[k]=d[k]-b[3*k+1]*d[k+1];
}
return(2);
}
//一般带型方程组的求解
//n-方程组的阶数
//l-系数矩阵的半带宽
//il-系数矩阵的带宽，应为il=2l+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!=(2*l+1))
{
printf("fail\n"); return(-2);
}
ls=l;
for (k=0;k<=n-2;k++)
{
p=0.0;
for (i=k;i<=ls;i++)
{
t=fabs(b[i*il]);
if (t>p)
{
p=t; is=i;
}
}
if (p+1.0==1.0)
{
printf("fail\n"); return(0);
}
for (j=0;j<=m-1;j++)
{
u=k*m+j; v=is*m+j;
t=d[u]; d[u]=d[v]; d[v]=t;
}
for (j=0;j<=il-1;j++)
{
u=k*il+j; v=is*il+j;
t=b[u]; b[u]=b[v]; b[v]=t;
}
for (j=0;j<=m-1;j++)
{
u=k*m+j; d[u]=d[u]/b[k*il];
}
for (j=1;j<=il-1;j++)
{
u=k*il+j; b[u]=b[u]/b[k*il];
}
for (i=k+1;i<=ls;i++)
{
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);
}
//求解对称方程组的分解法
//a n*n数组，存放方程组的系数矩阵（应为对称矩阵）
//n 方程组的阶数
//m 方程组右端常数向量的个数
//c n*m 存放方程组右端的m组常数向量；返回时存放方程组的m组解向量
//返回值小于0表示程序工作失败，大于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])+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);
}
//求解对称正定方程组的平方根法
//用Cholesky分解法（即平方根法）
//a n*n 存放系数矩阵（应为对称正定矩阵），返回时其上三角部分存放分解后的U矩阵
//n 方程的阶数
//m 方程组有短的常数向量的个数
//d n*m 存放方程组右端m组常数向量；返回时，存放方程组的m组解向量
//返回值小于0，表示程序工作失败
int chchol(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);
}
//求解大型稀疏方程组的全选主元高斯-约当消去法
//只是在消去过程中避免了对零元素的运算，从而可以大大减少运算次数
//本函数不考虑系数矩阵的压缩存储问题
//a n*n数组，存放方程组的系数矩阵，返回时要被破坏
//n 方程组的阶数
//d-存放方程组右端常数向量，返回时存放解向量
//返回值为0表示系数矩阵奇异
int ciggj(double a[],int n,double b[])
{
int *js,i,j,k,is,u,v;
double d,t;
js=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);
}
//求解Toeplitz型方程组的Levinson递推算法
//t 长度为n 的一维数组，存放对称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=malloc(n*sizeof(double));
y=malloc(n*sizeof(double));
a=t[0];
if (fabs(a)+1.0==1.0)
{
free(s); 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); 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); 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);
}
