分享
 
 
 

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

王朝other·作者佚名  2006-03-20
窄屏简体版  字體: |||超大  

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

//////////////////////////////////////////////////////////////

//实矩阵相乘

//计算矩阵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<m; i++)

{

for (j=0; j<k; j++)

{

u=i*k+j;

c[u]=0.0;

for (l=0; l<n; 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<n; i++)

{

for (j=0; j<k; j++)

{

u=i*k+j;

c[u]=0.0;

for (l=0; l<m; 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<m; i++)

{

for (j=0; j<k; j++)

{

u=i*k+j;

c[u]=0.0;

for (l=0; l<n; 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<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);

}

//奇异值分解法求广义逆

//本函数返回值小于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[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<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;

}

}

//复数矩阵相乘

//计算矩阵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);

}

 
 
 
免责声明:本文为网络用户发布,其观点仅代表作者个人观点,与本站无关,本站仅提供信息存储服务。文中陈述内容未经本站证实,其真实性、完整性、及时性本站不作任何保证或承诺,请读者仅作参考,并请自行核实相关内容。
2023年上半年GDP全球前十五强
 百态   2023-10-24
美众议院议长启动对拜登的弹劾调查
 百态   2023-09-13
上海、济南、武汉等多地出现不明坠落物
 探索   2023-09-06
印度或要将国名改为“巴拉特”
 百态   2023-09-06
男子为女友送行,买票不登机被捕
 百态   2023-08-20
手机地震预警功能怎么开?
 干货   2023-08-06
女子4年卖2套房花700多万做美容:不但没变美脸,面部还出现变形
 百态   2023-08-04
住户一楼被水淹 还冲来8头猪
 百态   2023-07-31
女子体内爬出大量瓜子状活虫
 百态   2023-07-25
地球连续35年收到神秘规律性信号,网友:不要回答!
 探索   2023-07-21
全球镓价格本周大涨27%
 探索   2023-07-09
钱都流向了那些不缺钱的人,苦都留给了能吃苦的人
 探索   2023-07-02
倩女手游刀客魅者强控制(强混乱强眩晕强睡眠)和对应控制抗性的关系
 百态   2020-08-20
美国5月9日最新疫情:美国确诊人数突破131万
 百态   2020-05-09
荷兰政府宣布将集体辞职
 干货   2020-04-30
倩女幽魂手游师徒任务情义春秋猜成语答案逍遥观:鹏程万里
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案神机营:射石饮羽
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案昆仑山:拔刀相助
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案天工阁:鬼斧神工
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案丝路古道:单枪匹马
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:与虎谋皮
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:李代桃僵
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:指鹿为马
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案金陵:小鸟依人
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案金陵:千金买邻
 干货   2019-11-12
 
推荐阅读
 
 
 
>>返回首頁<<
 
靜靜地坐在廢墟上,四周的荒凉一望無際,忽然覺得,淒涼也很美
© 2005- 王朝網路 版權所有