数值计算程序大放送-矩阵运算
//////////////////////////////////////////////////////////////
//实矩阵相乘
//计算矩阵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);
//////////////////////////////////////////////////////////////
选自<<徐世良数值计算程序集(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 { for (j=0; j { u=i*k+j; c[u]=0.0; for (l=0; l { c[u]+=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 { for (j=0; j { u=i*k+j; c[u]=0.0; for (l=0; l { c[u]+=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 { for (j=0; j { u=i*k+j; c[u]=0.0; for (l=0; l { c[u]+=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[u]; a[u]=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[u]; a[u]=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[u]=a[u]*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[u]=a[u]-a[i*n+k]*a[k*n+j]; } } } } for (i=0; i<=n-1; i++) { if (i!=k) { u=i*n+k; a[u]=-a[u]*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[u]; a[u]=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[u]; a[u]=a[v]; a[v]=p; } } } free(is); free(js); return(1); //对称正定矩阵求逆 //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[i]=g/w; if (i<=m) { b[i]=-b[i]; } 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[i]; } } 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[i]=c[i]+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[i]=r[i]+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[i]/a; b[j]=-c[i]/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[i]*b[j+1]; b[k]=b[k]+c[n-j-2]*b[n-i-1]; } } free(c); free(r); free(p); return(1); } //求矩阵的行列式值 //全选主元高斯消去法 //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[u]; a[u]=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[u]; a[u]=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[u]=a[u]-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[u]=a[u]/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[u]*a[u]; } 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[u]=a[u]-a[i*n+k]*a[j*n+k]; } a[u]=a[u]/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); } //求矩阵的行列式值 //全选主元高斯消去法 //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[u]; a[u]=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[u]; a[u]=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[u]=a[u]-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[u]=a[u]/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[u]*a[u]; } 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[u]=a[u]-a[i*n+k]*a[j*n+k]; } a[u]=a[u]/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); } //矩阵的三角分解(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 { 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 { 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); } //奇异值分解法求广义逆 //本函数返回值小于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 { 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 { k=m-1; } l=m; if (n-2 { 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 { mm=m+1; } if (k { s[k]=a[k*n+k]; } if (m { s[mm-1]=0.0; } if (l+1 { 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] { d=s[kk-1]; s[kk-1]=s[kk]; s[kk]=d; if (kk { 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 { 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[i]; s[i]=cs[0]*s[i]; 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[i]; s[i]=-cs[1]*e[i-1]+cs[0]*s[i]; fg[1]=cs[1]*e[i]; e[i]=cs[0]*e[i]; if (i { 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 { 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; } } //复数矩阵相乘 //计算矩阵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[u]=0.0; ci[u]=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[u]=cr[u]+p-q; ci[u]=ci[u]+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[u]*ar[u]+ai[u]*ai[u]; 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[u]; ar[u]=ar[v]; ar[v]=t; t=ai[u]; ai[u]=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[u]; ar[u]=ar[v]; ar[v]=t; t=ai[u]; ai[u]=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[u]*ar[l]; q=ai[u]*ai[l]; s=(ar[u]+ai[u])*(ar[l]+ai[l]); ar[u]=p-q; ai[u]=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[u]*ar[v]; q=ai[u]*ai[v]; s=(ar[u]+ai[u])*(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[u]*ar[l]; q=ai[u]*ai[l]; s=(ar[u]+ai[u])*(ar[l]+ai[l]); ar[u]=q-p; ai[u]=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[u]; ar[u]=ar[v]; ar[v]=t; t=ai[u]; ai[u]=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[u]; ar[u]=ar[v]; ar[v]=t; t=ai[u]; ai[u]=ai[v]; ai[v]=t; } } } free(is); free(js); return(1); }