数值计算程序大放送-矩阵运算

dormouse_buaa · 发表于 2005-1-19 22:34:43

<

>选自<<徐世良数值计算程序集(C)>>
<

>每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。
<

>今天都给贴出来了
#include "stdlib.h"
#include "math.h"
#include "stdio.h"
//实矩阵相乘
//计算矩阵A（m*n）和B（n*k）的乘积，结果保存在C（m*k）中
//a-长度为m*n的数组
//b-长度为n*k的数组
//c-长度为m*k的数组，存放结果
void damul(double a[],double b[],int m,int n,int k,double c[])
{
int i,j,l,u;
for (i=0; i<m; i++)
{
 for (j=0; j<k; j++)
 {
u=i*k+j;
c=0.0;
for (l=0; l<n; l++)
{
c+=a[i*n+l]*b[l*k+j];
}
 }
}
return;
}
//计算矩阵A（m*n）的转置矩阵AT（n*m）和B（m*k）的乘积，结果保存在C（n*k）中
//添加的函数，非原书程序
//a-长度为m*n的数组
//b-长度为m*k的数组
//c-长度为n*k的数组，存放结果
void ATdotB(double a[],double b[],int m,int n,int k,double c[])
{
int i,j,l,u;
for (i=0; i<n; i++)
{
 for (j=0; j<k; j++)
 {
u=i*k+j;
c=0.0;
for (l=0; l<m; l++)
{
c+=a[l*n+i]*b[l*k+j];
}
 }
}
return;
}
//计算矩阵A（m*n）和B（k*n）的转置矩阵BT（n*k）的乘积，结果保存在C（m*k）中
//添加的函数，非原书程序
//a-长度为m*n的数组
//b-长度为k*n的数组
//c-长度为m*k的数组，存放结果
void AdotBT(double a[],double b[],int m,int n,int k,double c[])
{
int i,j,l,u;
for (i=0; i<m; i++)
{
 for (j=0; j<k; j++)
 {
u=i*k+j;
c=0.0;
for (l=0; l<n; l++)
{
c+=a[i*n+l]*b[j*n+l];
}
 }
}
return;
}
//实矩阵求逆
//a-长度为n*n的数组, n*n矩阵
//n 矩阵的维数
int dcinv(double a[],int n)
{
int *is,*js,i,j,k,l,u,v;
double d,p;
is=new int[n];
is=malloc(n*sizeof(int));
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++)
{
l=i*n+j; p=fabs(a[l]);
if (p>d)
{
 d=p;
 is[k]=i;
 js[k]=j;
}
}
 }
 if (d+1.0==1.0)
 {
free(is);
free(js);
printf("err**not inv\n");
return(0);
 }
 if (is[k]!=k)
 {
for (j=0; j<=n-1; j++)
{
u=k*n+j; v=is[k]*n+j;
p=a; a=a[v]; a[v]=p;
}
 }
 if (js[k]!=k)
 {
for (i=0; i<=n-1; i++)
{
u=i*n+k; v=i*n+js[k];
p=a; a=a[v]; a[v]=p;
}
 }
 l=k*n+k;
 a[l]=1.0/a[l];
 for (j=0; j<=n-1; j++)
 {
if (j!=k)
{
u=k*n+j; a=a*a[l];
}
 }
 for (i=0; i<=n-1; i++)
 {
if (i!=k)
{
for (j=0; j<=n-1; j++)
{
 if (j!=k)
 {
 u=i*n+j;
 a=a-a[i*n+k]*a[k*n+j];
 }
}
}
 }
 for (i=0; i<=n-1; i++)
 {
if (i!=k)
{
u=i*n+k; a=-a*a[l];
}

 }

}
for (k=n-1; k>=0; k--)
{
 if (js[k]!=k)
 {
for (j=0; j<=n-1; j++)
{
u=k*n+j; v=js[k]*n+j;
p=a; a=a[v]; a[v]=p;
}
 }
 if (is[k]!=k)
 {
for (i=0; i<=n-1; i++)
{
u=i*n+k; v=i*n+is[k];
p=a; a=a[v]; a[v]=p;
}
 }
}
free(is);
free(js);
return(1);
}

dormouse_buaa · 发表于 2005-1-19 22:35:19

//对称正定矩阵求逆
//a-长度为n*n的数组, n*n矩阵
//n 矩阵的维数
int desgj(double a[],int n)
{
int i,j,k,m;
double w,g,*b;
b=malloc(n*sizeof(double));
for (k=0; k<=n-1; k++)
{
 w=a[0];
 if (fabs(w)+1.0==1.0)
 {
free(b);
printf("fail\n");
return(-2);
 }
 m=n-k-1;
 for (i=1; i<=n-1; i++)
 {
g=a[i*n];
b=g/w;
 if (i<=m)
{
b=-b;
}
 for (j=1; j<=i; j++)
{
a[(i-1)*n+j-1]=a[i*n+j]+g*b[j];
}
 }
 a[n*n-1]=1.0/w;
 for (i=1; i<=n-1; i++)
 {
a[(n-1)*n+i-1]=b;
 }
}
for (i=0; i<=n-2; i++)
{
 for (j=i+1; j<=n-1; j++)
 {
a[i*n+j]=a[j*n+i];
 }
}
free(b);
return(2);
}
//托伯利兹(Toeplitz)矩阵求逆的特兰持(Trench)方法
//t-长度为n的数组,存放n阶T型矩阵中的上三角元素t0,t1,t2...tn-1
//tt-长度为n的数组,从tt[1]开始依次存放tt[1]...tt[n-1]
//n-矩阵的阶数
//b-长度为n*n的数组，返回时存放逆矩阵
int dftrn(double t[],double tt[],int n,double b[])
{
int i,j,k;
double a,s,*c,*r,*p;
c=malloc(n*sizeof(double));
r=malloc(n*sizeof(double));
p=malloc(n*sizeof(double));
if (fabs(t[0])+1.0==1.0)
{
 free(c);
 free(r);
 free(p);
 printf("fail\n");
 return(-1);
}
a=t[0];
c[0]=tt[1]/t[0];
r[0]=t[1]/t[0];
for (k=0; k<=n-3; k++)
{
 s=0.0;
 for (j=1; j<=k+1; j++)
 {
s=s+c[k+1-j]*tt[j];
 }
 s=(s-tt[k+2])/a;
 for (i=0; i<=k; i++)
 {
p=c+s*r[k-i];
 }
 c[k+1]=-s;
 s=0.0;
 for (j=1; j<=k+1; j++)
 {
s=s+r[k+1-j]*t[j];
 }
 s=(s-t[k+2])/a;
 for (i=0; i<=k; i++)
 {
r=r+s*c[k-i];
 c[k-i]=p[k-i];
 }
 r[k+1]=-s;
 a=0.0;
 for (j=1; j<=k+2; j++)
 {
a=a+t[j]*c[j-1];
 }
 a=t[0]-a;
 if (fabs(a)+1.0==1.0)
 {
free(c);
free(r);
free(p);
 printf("fail\n");
return(-1);
 }
}
b[0]=1.0/a;
for (i=0; i<=n-2; i++)
{
 k=i+1;
 j=(i+1)*n;
 b[k]=-r/a;
 b[j]=-c/a;
}
for (i=0; i<=n-1; i++)
{
 for (j=0; j<=n-2; j++)
 {
k=(i+1)*n+j+1;
b[k]=b[i*n+j]-c*b[j+1];
b[k]=b[k]+c[n-j-2]*b[n-i-1];
 }
}
free(c);
free(r);
free(p);
return(1);
}

dormouse_buaa · 发表于 2005-1-19 22:36:28

//求矩阵的行列式值
//全选主元高斯消去法
//a-长度为n*n的数组
//n-矩阵的阶数
double dhdet(double a[],int n)
{
int i,j,k,is,js,l,u,v;
double f,det,q,d;
f=1.0; det=1.0;
for (k=0; k<=n-2; k++)
{
 q=0.0;
 for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{
l=i*n+j; d=fabs(a[l]);
if (d>q)
{
 q=d;
 is=i;
 js=j;
}
}
if (q+1.0==1.0)
{
det=0.0;
return(det);
}
if (is!=k)
{
f=-f;
for (j=k; j<=n-1; j++)
{
 u=k*n+j; v=is*n+j;
 d=a; a=a[v]; a[v]=d;
}
}
if (js!=k)
{
f=-f;
for (i=k; i<=n-1; i++)
{
 u=i*n+js; v=i*n+k;
 d=a; a=a[v]; a[v]=d;
}
}
l=k*n+k;
det=det*a[l];
for (i=k+1; i<=n-1; i++)
{
d=a[i*n+k]/a[l];
for (j=k+1; j<=n-1; j++)
{
 u=i*n+j;
 a=a-d*a[k*n+j];
}
}
}
det=f*det*a[n*n-1];
return(det);
}
//对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值
//返回值小于0表示程序工作失败，还输出"fail";
//返回值大于0表示正常返回
//a-长度为n*n的数组，存放正定矩阵，
// 返回时下三角部分存放分解后的下三角矩阵L,其余元素为0
//n-正定矩阵的阶数
//det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值
int dicll(double a[],int n,double *det)
{
int i,j,k,u,v,l;
double d;
if ((a[0]+1.0==1.0)||(a[0]<0.0))
{
 printf("fail\n");
 return(-2);
}
a[0]=sqrt(a[0]);
d=a[0];
for (i=1; i<=n-1; i++)
{
 u=i*n;
 a=a/a[0];
}
for (j=1; j<=n-1; j++)
{
 l=j*n+j;
 for (k=0; k<=j-1; k++)
 {
u=j*n+k;
a[l]=a[l]-a*a;
 }
 if ((a[l]+1.0==1.0)||(a[l]<0.0))
 {
printf("fail\n");
return(-2);
 }
 a[l]=sqrt(a[l]);
 d=d*a[l];
 for (i=j+1; i<=n-1; i++)
 {
u=i*n+j;
 for (k=0; k<=j-1; k++)
{
a=a-a[i*n+k]*a[j*n+k];
}
 a=a/a[l];
 }
}
*det=d*d;
for (i=0; i<=n-2; i++)
{
 for (j=i+1; j<=n-1; j++)
 {
a[i*n+j]=0.0;
 }
}
return(2);
}

dormouse_buaa · 发表于 2005-1-19 22:36:29

//求矩阵的行列式值
//全选主元高斯消去法
//a-长度为n*n的数组
//n-矩阵的阶数
double dhdet(double a[],int n)
{
int i,j,k,is,js,l,u,v;
double f,det,q,d;
f=1.0; det=1.0;
for (k=0; k<=n-2; k++)
{
 q=0.0;
 for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{
l=i*n+j; d=fabs(a[l]);
if (d>q)
{
 q=d;
 is=i;
 js=j;
}
}
if (q+1.0==1.0)
{
det=0.0;
return(det);
}
if (is!=k)
{
f=-f;
for (j=k; j<=n-1; j++)
{
 u=k*n+j; v=is*n+j;
 d=a; a=a[v]; a[v]=d;
}
}
if (js!=k)
{
f=-f;
for (i=k; i<=n-1; i++)
{
 u=i*n+js; v=i*n+k;
 d=a; a=a[v]; a[v]=d;
}
}
l=k*n+k;
det=det*a[l];
for (i=k+1; i<=n-1; i++)
{
d=a[i*n+k]/a[l];
for (j=k+1; j<=n-1; j++)
{
 u=i*n+j;
 a=a-d*a[k*n+j];
}
}
}
det=f*det*a[n*n-1];
return(det);
}
//对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值
//返回值小于0表示程序工作失败，还输出"fail";
//返回值大于0表示正常返回
//a-长度为n*n的数组，存放正定矩阵，
// 返回时下三角部分存放分解后的下三角矩阵L,其余元素为0
//n-正定矩阵的阶数
//det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值
int dicll(double a[],int n,double *det)
{
int i,j,k,u,v,l;
double d;
if ((a[0]+1.0==1.0)||(a[0]<0.0))
{
 printf("fail\n");
 return(-2);
}
a[0]=sqrt(a[0]);
d=a[0];
for (i=1; i<=n-1; i++)
{
 u=i*n;
 a=a/a[0];
}
for (j=1; j<=n-1; j++)
{
 l=j*n+j;
 for (k=0; k<=j-1; k++)
 {
u=j*n+k;
a[l]=a[l]-a*a;
 }
 if ((a[l]+1.0==1.0)||(a[l]<0.0))
 {
printf("fail\n");
return(-2);
 }
 a[l]=sqrt(a[l]);
 d=d*a[l];
 for (i=j+1; i<=n-1; i++)
 {
u=i*n+j;
 for (k=0; k<=j-1; k++)
{
a=a-a[i*n+k]*a[j*n+k];
}
 a=a/a[l];
 }
}
*det=d*d;
for (i=0; i<=n-2; i++)
{
 for (j=i+1; j<=n-1; j++)
 {
a[i*n+j]=0.0;
 }
}
return(2);
}

dormouse_buaa · 发表于 2005-1-19 22:37:28

<

>
//矩阵的三角分解(LU)
//其中下三角阵L的主对角元素为1。
//a-长度为N*N的矩阵，返回时为L+U-I
//n-矩阵的阶数
//l-返回下三角矩阵
//u-返回上三角矩阵
int djlu(double a[],int n,double l[],double u[])
{
int i,j,k,w,v,ll;
for (k=0; k<=n-2; k++)
{
 ll=k*n+k;
 if (fabs(a[ll])+1.0==1.0)
 {
printf("fail\n");
return(0);
 }
 for (i=k+1; i<=n-1; i++)
 {
w=i*n+k;
a[w]=a[w]/a[ll];
 }
 for (i=k+1; i<=n-1; i++)
 {
w=i*n+k;
for (j=k+1; j<=n-1; j++)
{
v=i*n+j;
a[v]=a[v]-a[w]*a[k*n+j];
}
 }
}
for (i=0; i<=n-1; i++)
{
 for (j=0; j<i; j++)
 {
w=i*n+j;
l[w]=a[w];
u[w]=0.0;
 }
 w=i*n+i;
 l[w]=1.0;
 u[w]=a[w];
 for (j=i+1; j<=n-1; j++)
 {
w=i*n+j;
l[w]=0.0;
u[w]=a[w];
 }
}
return(1);
}
//实数矩阵的QR分解法
//用Householder变换对一般m*n阶的实数矩阵进行QR分解
//a-长度为m*n的一维数组,返回时其左上三角部分存放上三角矩阵R
//m-矩阵的行数
//n-矩阵的列数
//q-长度为m*m的矩阵，返回时存放正交矩阵Q
int dkqr(double a[],int m,int n,double q[])
{
int i,j,k,l,nn,p,jj;
double u,alpha,w,t;
if (m<n)
{
 printf("fail\n");
 return(0);
}
for (i=0; i<=m-1; i++)
{
 for (j=0; j<=m-1; j++)
 {
l=i*m+j;
q[l]=0.0;
if (i==j)
{
q[l]=1.0;
}
 }
}
nn=n;
if (m==n)
{
 nn=m-1;
}
for (k=0; k<=nn-1; k++)
{
 u=0.0;
 l=k*n+k;
 for (i=k; i<=m-1; i++)
 {
w=fabs(a[i*n+k]);
if (w>u)
{
u=w;
}
 }
 alpha=0.0;
 for (i=k; i<=m-1; i++)
 {
t=a[i*n+k]/u;
alpha=alpha+t*t;
 }
 if (a[l]>0.0)
 {
u=-u;
 }
 alpha=u*sqrt(alpha);
 if (fabs(alpha)+1.0==1.0)
 {
printf("fail\n");
return(0);
 }
 u=sqrt(2.0*alpha*(alpha-a[l]));
 if ((u+1.0)!=1.0)
 {
a[l]=(a[l]-alpha)/u;
for (i=k+1; i<=m-1; i++)
{
p=i*n+k;
a[p]=a[p]/u;
}
for (j=0; j<=m-1; j++)
{
t=0.0;
for (jj=k; jj<=m-1; jj++)
{
 t=t+a[jj*n+k]*q[jj*m+j];
}
for (i=k; i<=m-1; i++)
{
 p=i*m+j;
 q[p]=q[p]-2.0*t*a[i*n+k];
}
}
for (j=k+1; j<=n-1; j++)
{
t=0.0;
for (jj=k; jj<=m-1; jj++)
{
 t=t+a[jj*n+k]*a[jj*n+j];
}
for (i=k; i<=m-1; i++)
{
 p=i*n+j;
 a[p]=a[p]-2.0*t*a[i*n+k];
}
}
a[l]=alpha;
for (i=k+1; i<=m-1; i++)
{
a[i*n+k]=0.0;
}
 }
}
for (i=0; i<=m-2; i++)
{
 for (j=i+1; j<=m-1;j++)
 {
p=i*m+j; l=j*m+i;
t=q[p]; q[p]=q[l]; q[l]=t;
 }
}
return(1);
}

dormouse_buaa · 发表于 2005-1-19 22:38:13

<

>
//奇异值分解法求广义逆
//本函数返回值小于0表示在奇异值分解过程,
//中迭代值超过了60次还未满足精度要求.
//返回值大于0表示正常返回。
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//aa-长度为n*m的数组，返回式存放A的广义逆
//eps-精度要求
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav()
int dginv(double a[],int m,int n,double aa[],double eps,double u[],double v[],int ka)
{
int i,j,k,l,t,p,q,f;
i=dluav(a,m,n,u,v,eps,ka);
if (i<0)
{
 return(-1);
}
j=n;
if (m<n)
{
 j=m;
}
j=j-1;
k=0;
while ((k<=j)&&(a[k*n+k]!=0.0))
{
 k=k+1;
}
k=k-1;
for (i=0; i<=n-1; i++)
{
 for (j=0; j<=m-1; j++)
 {
t=i*m+j;
aa[t]=0.0;
for (l=0; l<=k; l++)
{
f=l*n+i;
p=j*m+l;
q=l*n+l;
aa[t]=aa[t]+v[f]*u[p]/a[q];
}
 }
}
return(1);
}
//实数矩阵的奇异值分解
//利用Householder变换及变形QR算法
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//eps-精度要求
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav(),ppp(),sss()
static void ppp(double a[],double e[],double s[],double v[],int m,int n);
static void sss(double fg[2],double cs[2]);
int dluav(double a[],int m,int n,double u[],double v[],double eps,int ka)
{
int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,ks;
double d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];
double *s,*e,*w;
s=malloc(ka*sizeof(double));
e=malloc(ka*sizeof(double));
w=malloc(ka*sizeof(double));
it=60;
k=n;
if (m-1<n)
{
 k=m-1;
}
l=m;
if (n-2<m)
{
 l=n-2;
}
if (l<0)
{
 l=0;
}
ll=k;
if (l>k)
{
 ll=l;
}
if (ll>=1)
{
 for (kk=1; kk<=ll; kk++)
 {
if (kk<=k)
{
d=0.0;
for (i=kk; i<=m; i++)
{
 ix=(i-1)*n+kk-1;
 d=d+a[ix]*a[ix];
}
s[kk-1]=sqrt(d);
if (s[kk-1]!=0.0)
{
 ix=(kk-1)*n+kk-1;
 if (a[ix]!=0.0)
 {
 s[kk-1]=fabs(s[kk-1]);
 if (a[ix]<0.0)
 {
 s[kk-1]=-s[kk-1];
 }
 }
 for (i=kk; i<=m; i++)
 {
 iy=(i-1)*n+kk-1;
 a[iy]=a[iy]/s[kk-1];
 }
 a[ix]=1.0+a[ix];
}
s[kk-1]=-s[kk-1];
}
if (n>=kk+1)
{
for (j=kk+1; j<=n; j++)
{
 if ((kk<=k)&&(s[kk-1]!=0.0))
 {
 d=0.0;
 for (i=kk; i<=m; i++)
 {
 ix=(i-1)*n+kk-1;
 iy=(i-1)*n+j-1;
 d=d+a[ix]*a[iy];
 }
 d=-d/a[(kk-1)*n+kk-1];
 for (i=kk; i<=m; i++)
 {
 ix=(i-1)*n+j-1;
 iy=(i-1)*n+kk-1;
 a[ix]=a[ix]+d*a[iy];
 }
 }
 e[j-1]=a[(kk-1)*n+j-1];
}
}
if (kk<=k)
{
for (i=kk; i<=m; i++)
{
 ix=(i-1)*m+kk-1;
 iy=(i-1)*n+kk-1;
 u[ix]=a[iy];
}
}
if (kk<=l)
{
d=0.0;
for (i=kk+1; i<=n; i++)
{
 d=d+e[i-1]*e[i-1];
}
e[kk-1]=sqrt(d);
if (e[kk-1]!=0.0)
{
 if (e[kk]!=0.0)
 {
 e[kk-1]=fabs(e[kk-1]);
 if (e[kk]<0.0)
 {
 e[kk-1]=-e[kk-1];
 }
 }
 for (i=kk+1; i<=n; i++)
 {
 e[i-1]=e[i-1]/e[kk-1];
 }
 e[kk]=1.0+e[kk];
}
e[kk-1]=-e[kk-1];
if ((kk+1<=m)&&(e[kk-1]!=0.0))
{
 for (i=kk+1; i<=m; i++)
 {
 w[i-1]=0.0;
 }
 for (j=kk+1; j<=n; j++)
 {
 for (i=kk+1; i<=m; i++)
 {
 w[i-1]=w[i-1]+e[j-1]*a[(i-1)*n+j-1];
 }
 }
 for (j=kk+1; j<=n; j++)
 {
 for (i=kk+1; i<=m; i++)
 {
 ix=(i-1)*n+j-1;
 a[ix]=a[ix]-w[i-1]*e[j-1]/e[kk];
 }
 }
}
for (i=kk+1; i<=n; i++)
{
 v[(i-1)*n+kk-1]=e[i-1];
}
}
 }
}
mm=n;
if (m+1<n)
{
 mm=m+1;
}
if (k<n)
{
 s[k]=a[k*n+k];
}
if (m<mm)
{
 s[mm-1]=0.0;
}
if (l+1<mm)
{
 e[l]=a[l*n+mm-1];
}
e[mm-1]=0.0;
nn=m;
if (m>n)
{
 nn=n;
}
if (nn>=k+1)
{
 for (j=k+1; j<=nn; j++)
 {
for (i=1; i<=m; i++)
{
u[(i-1)*m+j-1]=0.0;
}
u[(j-1)*m+j-1]=1.0;
 }
}
if (k>=1)
{
 for (ll=1; ll<=k; ll++)
 {
kk=k-ll+1; iz=(kk-1)*m+kk-1;
if (s[kk-1]!=0.0)
{
if (nn>=kk+1)
{
 for (j=kk+1; j<=nn; j++)
 {
 d=0.0;
 for (i=kk; i<=m; i++)
 {
 ix=(i-1)*m+kk-1;
 iy=(i-1)*m+j-1;
 d=d+u[ix]*u[iy]/u[iz];
 }
 d=-d;
 for (i=kk; i<=m; i++)
 {
 ix=(i-1)*m+j-1;
 iy=(i-1)*m+kk-1;
 u[ix]=u[ix]+d*u[iy];
 }
 }
}
for (i=kk; i<=m; i++)
{
 ix=(i-1)*m+kk-1;
 u[ix]=-u[ix];
}
u[iz]=1.0+u[iz];
if (kk-1>=1)
{
 for (i=1; i<=kk-1; i++)
 {
 u[(i-1)*m+kk-1]=0.0;
 }
}
}
else
{
for (i=1; i<=m; i++)
{
 u[(i-1)*m+kk-1]=0.0;
}
u[(kk-1)*m+kk-1]=1.0;
}
 }
}
for (ll=1; ll<=n; ll++)
{
 kk=n-ll+1;
 iz=kk*n+kk-1;
 if ((kk<=l)&&(e[kk-1]!=0.0))
 {
for (j=kk+1; j<=n; j++)
{
d=0.0;
for (i=kk+1; i<=n; i++)
{
 ix=(i-1)*n+kk-1;
 iy=(i-1)*n+j-1;
 d=d+v[ix]*v[iy]/v[iz];
}
d=-d;
for (i=kk+1; i<=n; i++)
{
 ix=(i-1)*n+j-1;
 iy=(i-1)*n+kk-1;
 v[ix]=v[ix]+d*v[iy];
}
}
 }
 for (i=1; i<=n; i++)
 {
v[(i-1)*n+kk-1]=0.0;
 }
 v[iz-n]=1.0;
}
for (i=1; i<=m; i++)
{
 for (j=1; j<=n; j++)
 {
a[(i-1)*n+j-1]=0.0;
 }
}
m1=mm;
it=60;
while (1==1)
{
 if (mm==0)
 {
ppp(a,e,s,v,m,n);
free(s);
free(e);
free(w);
return(1);
 }
 if (it==0)
 {
ppp(a,e,s,v,m,n);
free(s);
free(e);
free(w);
return(-1);
 }
 kk=mm-1;
 while ((kk!=0)&&(fabs(e[kk-1])!=0.0))
 {
d=fabs(s[kk-1])+fabs(s[kk]);
dd=fabs(e[kk-1]);
if (dd>eps*d)
{
kk=kk-1;
}
else
{
e[kk-1]=0.0;
}
 }
 if (kk==mm-1)
 {
kk=kk+1;
if (s[kk-1]<0.0)
{
s[kk-1]=-s[kk-1];
for (i=1; i<=n; i++)
{
 ix=(i-1)*n+kk-1;
 v[ix]=-v[ix];
}
}
while ((kk!=m1)&&(s[kk-1]<s[kk]))
{
d=s[kk-1];
s[kk-1]=s[kk];
s[kk]=d;
if (kk<n)
{
 for (i=1; i<=n; i++)
 {
 ix=(i-1)*n+kk-1;
 iy=(i-1)*n+kk;
 d=v[ix];
 v[ix]=v[iy];
 v[iy]=d;
 }
}
if (kk<m)
{
 for (i=1; i<=m; i++)
 {
 ix=(i-1)*m+kk-1; iy=(i-1)*m+kk;
 d=u[ix]; u[ix]=u[iy]; u[iy]=d;
 }
}
kk=kk+1;
}
it=60;
mm=mm-1;
 }
 else
 {
ks=mm;
while ((ks>kk)&&(fabs(s[ks-1])!=0.0))
{
d=0.0;
if (ks!=mm)
{
 d=d+fabs(e[ks-1]);
}
if (ks!=kk+1)
{
 d=d+fabs(e[ks-2]);
}
dd=fabs(s[ks-1]);
if (dd>eps*d)
{
 ks=ks-1;
}
else
{
 s[ks-1]=0.0;
}
}
if (ks==kk)
{
kk=kk+1;
d=fabs(s[mm-1]);
t=fabs(s[mm-2]);
if (t>d)
{
 d=t;
}
t=fabs(e[mm-2]);
if (t>d)
{
 d=t;
}
t=fabs(s[kk-1]);
if (t>d)
{
 d=t;
}
t=fabs(e[kk-1]);
if (t>d)
{
 d=t;
}
sm=s[mm-1]/d;
sm1=s[mm-2]/d;
em1=e[mm-2]/d;
sk=s[kk-1]/d;
ek=e[kk-1]/d;
b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0;
c=sm*em1;
c=c*c;
shh=0.0;
if ((b!=0.0)||(c!=0.0))
{
 shh=sqrt(b*b+c);
 if (b<0.0)
 {
 shh=-shh;
 }
 shh=c/(b+shh);
}
fg[0]=(sk+sm)*(sk-sm)-shh;
fg[1]=sk*ek;
for (i=kk; i<=mm-1; i++)
{
 sss(fg,cs);
 if (i!=kk)
 {
 e[i-2]=fg[0];
 }
 fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1];
 e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1];
 fg[1]=cs[1]*s;
 s=cs[0]*s;
 if ((cs[0]!=1.0)||(cs[1]!=0.0))
 {
 for (j=1; j<=n; j++)
 {
 ix=(j-1)*n+i-1;
 iy=(j-1)*n+i;
 d=cs[0]*v[ix]+cs[1]*v[iy];
 v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];
 v[ix]=d;
 }
 }
 sss(fg,cs);
 s[i-1]=fg[0];
 fg[0]=cs[0]*e[i-1]+cs[1]*s;
 s=-cs[1]*e[i-1]+cs[0]*s;
 fg[1]=cs[1]*e;
 e=cs[0]*e;
 if (i<m)
 {
 if ((cs[0]!=1.0)||(cs[1]!=0.0))
 {
 for (j=1; j<=m; j++)
 {
 ix=(j-1)*m+i-1;
 iy=(j-1)*m+i;
 d=cs[0]*u[ix]+cs[1]*u[iy];
 u[iy]=-cs[1]*u[ix]+cs[0]*u[iy];
 u[ix]=d;
 }
 }
 }
}
e[mm-2]=fg[0];
it=it-1;
}
else
{
if (ks==mm)
{
 kk=kk+1;
 fg[1]=e[mm-2];
 e[mm-2]=0.0;
 for (ll=kk; ll<=mm-1; ll++)
 {
 i=mm+kk-ll-1;
 fg[0]=s[i-1];
 sss(fg,cs);
 s[i-1]=fg[0];
 if (i!=kk)
 {
 fg[1]=-cs[1]*e[i-2];
 e[i-2]=cs[0]*e[i-2];
 }
 if ((cs[0]!=1.0)||(cs[1]!=0.0))
 {
 for (j=1; j<=n; j++)
 {
 ix=(j-1)*n+i-1;
 iy=(j-1)*n+mm-1;
 d=cs[0]*v[ix]+cs[1]*v[iy];
 v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];
 v[ix]=d;
 }
 }
 }
}
else
{
 kk=ks+1;
 fg[1]=e[kk-2];
 e[kk-2]=0.0;
 for (i=kk; i<=mm; i++)
 {
 fg[0]=s[i-1];
 sss(fg,cs);
 s[i-1]=fg[0];
 fg[1]=-cs[1]*e[i-1];
 e[i-1]=cs[0]*e[i-1];
 if ((cs[0]!=1.0)||(cs[1]!=0.0))
 {
 for (j=1; j<=m; j++)
 {
 ix=(j-1)*m+i-1;
 iy=(j-1)*m+kk-2;
 d=cs[0]*u[ix]+cs[1]*u[iy];
 u[iy]=-cs[1]*u[ix]+cs[0]*u[iy];
 u[ix]=d;
 }
 }
 }
}
}
 }
}
return(1);
}<>static void ppp(double a[],double e[],double s[],double v[],int m,int n)
{
int i,j,p,q;
double d;
if (m>=n)
{
 i=n;
}
else
{
 i=m;
}
for (j=1; j<=i-1; j++)
{
 a[(j-1)*n+j-1]=s[j-1];
 a[(j-1)*n+j]=e[j-1];
}
a[(i-1)*n+i-1]=s[i-1];
if (m<n)
{
 a[(i-1)*n+i]=e[i-1];
}
for (i=1; i<=n-1; i++)
{
 for (j=i+1; j<=n; j++)
 {
p=(i-1)*n+j-1;
q=(j-1)*n+i-1;
d=v[p];
v[p]=v[q];
v[q]=d;
 }
}
return;
}<>static void sss(double fg[2],double cs[2])
{
double r,d;
if ((fabs(fg[0])+fabs(fg[1]))==0.0)
{
 cs[0]=1.0;
 cs[1]=0.0;
 d=0.0;
}
else
{
 d=sqrt(fg[0]*fg[0]+fg[1]*fg[1]);
 if (fabs(fg[0])>fabs(fg[1]))
 {
d=fabs(d);
if (fg[0]<0.0)
{
d=-d;
}
 }
 if (fabs(fg[1])>=fabs(fg[0]))
 {
d=fabs(d);
if (fg[1]<0.0)
{
d=-d;
}
 }
 cs[0]=fg[0]/d;
 cs[1]=fg[1]/d;
}
r=1.0;
if (fabs(fg[0])>fabs(fg[1]))
{
 r=cs[1];
}
else
{
 if(cs[0]!=0.0)
 {
r=1.0/cs[0];
 }
 fg[0]=d;
 fg[1]=r;
 return;
}
}

dormouse_buaa · 发表于 2005-1-19 22:39:31

//复数矩阵相乘
//计算矩阵A（m*n）和B（n*k）的乘积，结果保存在C（m*k）中
//ar-长度为m*n的数组,存放A的实部
//ai-长度为m*n的数组,存放A的虚部
//br-长度为n*k的数组,存放B的实部
//bi-长度为n*k的数组,存放B的虚部
//cr-长度为m*k的数组，存放结果C的实部
//ci-长度为m*k的数组，存放结果C的虚部
void dbcmul(double ar[],double ai[],double br[],double bi[],int m,int n,int k,double cr[],double ci[])
{
int i,j,l,u,v,w;
double p,q,s;
for (i=0; i<=m-1; i++)
{
 for (j=0; j<=k-1; j++)
 {
u=i*k+j;
cr=0.0;
ci=0.0;
for (l=0; l<=n-1; l++)
{
v=i*n+l; w=l*k+j;
p=ar[v]*br[w];
q=ai[v]*bi[w];
s=(ar[v]+ai[v])*(br[w]+bi[w]);
cr=cr+p-q;
ci=ci+s-p-q;
}
 }
}
 return;
}
//复数矩阵求逆
//ar-长度为n*n的数组, n*n矩阵的实部
//ai-长度为n*n的数组, n*n矩阵的虚部
//n 矩阵的维数
int ddcinv(double ar[],double ai[],int n)
{
int *is,*js,i,j,k,l,u,v,w;
double p,q,s,t,d,b;
is=malloc(n*sizeof(int));
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*ar+ai*ai;
if (p>d)
{
 d=p;
 is[k]=i;
 js[k]=j;
}
}
 }
 if (d+1.0==1.0)
 {
free(is);
free(js);
printf("err**not inv\n");
return(0);
 }
 if (is[k]!=k)
 {
for (j=0; j<=n-1; j++)
{
u=k*n+j;
v=is[k]*n+j;
t=ar; ar=ar[v]; ar[v]=t;
t=ai; ai=ai[v]; ai[v]=t;
}
 }
 if (js[k]!=k)
 {
for (i=0; i<=n-1; i++)
{
u=i*n+k;
v=i*n+js[k];
t=ar; ar=ar[v]; ar[v]=t;
t=ai; ai=ai[v]; ai[v]=t;
}
 }
 l=k*n+k;
 ar[l]=ar[l]/d;
 ai[l]=-ai[l]/d;
 for (j=0; j<=n-1; j++)
 {
if (j!=k)
{
u=k*n+j;
p=ar*ar[l];
q=ai*ai[l];
s=(ar+ai)*(ar[l]+ai[l]);
ar=p-q;
ai=s-p-q;
}
 }
 for (i=0; i<=n-1; i++)
 {
if (i!=k)
{
v=i*n+k;
for (j=0; j<=n-1; j++)
{
 if (j!=k)
 {
 u=k*n+j;
 w=i*n+j;
 p=ar*ar[v];
 q=ai*ai[v];
 s=(ar+ai)*(ar[v]+ai[v]);
 t=p-q;
 b=s-p-q;
 ar[w]=ar[w]-t;
 ai[w]=ai[w]-b;
 }
}
}
 }
 for (i=0; i<=n-1; i++)
 {
if (i!=k)
{
u=i*n+k;
p=ar*ar[l];
q=ai*ai[l];
s=(ar+ai)*(ar[l]+ai[l]);
ar=q-p;
ai=p+q-s;
}
 }
}
for (k=n-1; k>=0; k--)
{
 if (js[k]!=k)
 {
for (j=0; j<=n-1; j++)
 {
u=k*n+j;
v=js[k]*n+j;
t=ar; ar=ar[v]; ar[v]=t;
t=ai; ai=ai[v]; ai[v]=t;
 }
 }
 if (is[k]!=k)
 {
for (i=0; i<=n-1; i++)
{
u=i*n+k;
v=i*n+is[k];
t=ar; ar=ar[v]; ar[v]=t;
t=ai; ai=ai[v]; ai[v]=t;
}
 }
}
free(is);
free(js);
return(1);
}

dormouse_buaa · 发表于 2005-1-19 22:39:53

<

>//////////////////////////////////////////////////////////////
//实矩阵相乘
//计算矩阵A（m*n）和B（n*k）的乘积，结果保存在C（m*k）中
//a-长度为m*n的数组
//b-长度为n*k的数组
//c-长度为m*k的数组，存放结果
void damul(double a[],double b[],int m,int n,int k,double c[]);
//////////////////////////////////////////////////////////////
//计算矩阵A（m*n）的转置矩阵AT（n*m）和B（m*k）的乘积，结果保存在C（n*k）中
//添加的函数，非原书程序
//a-长度为m*n的数组
//b-长度为m*k的数组
//c-长度为n*k的数组，存放结果
void ATdotB(double a[],double b[],int m,int n,int k,double c[]);
//////////////////////////////////////////////////////////////
//计算矩阵A（m*n）和B（k*n）的转置矩阵BT（n*k）的乘积，结果保存在C（m*k）中
//添加的函数，非原书程序
//a-长度为m*n的数组
//b-长度为k*n的数组
//c-长度为m*k的数组，存放结果
void AdotBT(double a[],double b[],int m,int n,int k,double c[]);
//////////////////////////////////////////////////////////////
//实矩阵求逆
//全选主元高斯-约当法
//a-长度为n*n的数组, n*n矩阵
//n 矩阵的维数
int dcinv(double a[],int n);
//////////////////////////////////////////////////////////////
//对称正定矩阵求逆
//a-长度为n*n的数组, n*n矩阵
//n 矩阵的维数
int desgj(double a[],int n);
//////////////////////////////////////////////////////////////
//托伯利兹(Toeplitz)矩阵求逆的特兰持(Trench)方法
//t-长度为n的数组,存放n阶T型矩阵中的上三角元素t0,t1,t2...tn-1
//tt-长度为n的数组,从tt[1]开始依次存放tt[1]...tt[n-1]
//n-矩阵的阶数
//b-长度为n*n的数组，返回时存放逆矩阵
int dftrn(double t[],double tt[],int n,double b[]);
//////////////////////////////////////////////////////////////
//求矩阵的行列式值
//全选主元高斯消去法
//a-长度为n*n的数组
//n-矩阵的阶数
double dhdet(double a[],int n);
//////////////////////////////////////////////////////////////
//对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值
//返回值小于0表示程序工作失败，还输出"fail";
//返回值大于0表示正常返回
//a-长度为n*n的数组，存放正定矩阵，
// 返回时下三角部分存放分解后的下三角矩阵L,其余元素为0
//n-正定矩阵的阶数
//det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值
int dicll(double a[],int n,double *det);
//////////////////////////////////////////////////////////////
//矩阵的三角分解(LU)
//其中下三角阵L的主对角元素为1。
//a-长度为N*N的矩阵，返回时为L+U-I
//n-矩阵的阶数
//l-返回下三角矩阵
//u-返回上三角矩阵
int djlu(double a[],int n,double l[],double u[]);
//////////////////////////////////////////////////////////////
//实数矩阵的QR分解法
//用Householder变换对一般m*n阶的实数矩阵进行QR分解
//a-长度为m*n的一维数组,返回时其左上三角部分存放上三角矩阵R
//m-矩阵的行数
//n-矩阵的列数
//q-长度为m*m的矩阵，返回时存放正交矩阵Q
int dkqr(double a[],int m,int n,double q[]);
//////////////////////////////////////////////////////////////
//奇异值分解法求广义逆
//本函数返回值小于0表示在奇异值分解过程,
//中迭代值超过了60次还未满足精度要求.
//返回值大于0表示正常返回。
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//aa-长度为n*m的数组，返回式存放A的广义逆
//eps-精度要求
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav()
int dginv(double a[],int m,int n,double aa[],double eps,double u[],double v[],int ka);
//////////////////////////////////////////////////////////////
//实数矩阵的奇异值分解
//利用Householder变换及变形QR算法
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//eps-精度要求
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav(),ppp(),sss()
int dluav(double a[],int m,int n,double u[],double v[],double eps,int ka);
//////////////////////////////////////////////////////////////
<

>//////////////////////////////////////////////////////////////
<

>//////////////////////////////////////////////////////////////

dormouse_buaa · 发表于 2005-1-19 22:40:22

//////////////////////////////////////////////////////////////
//复数矩阵相乘
//计算矩阵A（m*n）和B（n*k）的乘积，结果保存在C（m*k）中
//ai-长度为m*n的数组,存放A的虚部
//br-长度为n*k的数组,存放B的实部
//bi-长度为n*k的数组,存放B的虚部
//cr-长度为m*k的数组，存放结果C的实部
//ci-长度为m*k的数组，存放结果C的虚部
void dbcmul(double ar[],double ai[],double br[],double bi[],int m,int n,int k,double cr[],double ci[]);
//////////////////////////////////////////////////////////////
//复数矩阵求逆
//全选主元高斯-约当法
//ar-长度为n*n的数组, n*n矩阵的实部
//ai-长度为n*n的数组, n*n矩阵的虚部
//n 矩阵的维数
int ddcinv(double ar[],double ai[],int n);
//////////////////////////////////////////////////////////////

		自动登录	找回密码
密码			注-册-帐-号