//奇异值分解
#include "stdlib.h"
#include"stdio.h"
#include"math.h"
#define pi 3.1415926
#define nx 225 //道数
#define nz 1501 //采样点数
void mbbub(double p[],int n);
int bmuav(double a[],int m,int n,double u[],double v[],double eps,int ka);
void damul(double a[],double b[],int m,int n,int k,double c[]);
void matr(double a[],double b[],int m,int n);
void main()
{
FILE *fp,*fp1,*fp2;
if(fopen_s(&fp,"02.dat","rb")!=NULL)// 600_2501
{
printf("\nfile can not be open\n");
exit(0);
}
fopen_s(&fp1,"px.txt","w");
fopen_s(&fp2,"11mm.dat","wb");
int i,j,k,b[nx],ij;
//由于奇异值分解是对矩阵进行的，所以读入时要以矩阵形式进行
//奇异值分解的餐参数设定
double *a,*u,*v,p[nx],*c,*d;
float *jilu=new float[nx*nz];
a=new double[nx*nz];
u=new double[nx*nx];
v=new double[nz*nz];
c=new double[nx*nz];
d=new double[nx*nz];
float temp=0.0;
for(i=0;i<NX;I++) if(a[i] { for(i="0;i<nx*nz;i++)" k="0;" printf(?='==bb====cond=%d\n",cond);' cond="bmuav(a,nx,nz,u,v,eps,nz+1);" cond; int printf(?---aa----\n?); eps="0.000001;" eps; double } jilu[i*nz+j]="temp;" a[i*nz+j]="temp;" fread(&temp,sizeof(float),1,fp); for(j="0;j<nz;j++)">0.0)
{
p[k]=a[i];
b[k]=i;
k++;
}
ij=k;
}
printf(" %d \n ",ij);
for(i=0;i<IJ;I++) { for(i="0;i<nx*nx;i++)" } if((fabs(u[i])<eps)||(fabs(u[i]) b[i],p[i]); %f\n?, %d fprintf(fp1,?>1/eps))
{
u[i]=0.0;
}
}
for(i=0;i<NZ*NZ;I++) { if((fabs(v[i])<eps)||(fabs(v[i])>1/eps))
{
v[i]=0.0;
}
}
//对输出值进行处理去除前半部分较大值//----结果是对原始数据顶端进行了除了
/* for(i=0;i<NX*NZ;I++) { if(fabs(a[i])>1000)
a[i]=0.0;
}*/
//对输出值进行处理去除后半部分较大值//----结果只对随机干扰有影响
for(i=0;i<NX*NZ;I++) { for(i="0;i<nx*nz;i++)" } if((fabs(a[i])<1000)&&(fabs(a[i]) * ---结果只含有面波成分 对输出值进行处理去除（100，1000）值 a[i]="0.0;" if(fabs(a[i])<10000)>100))
a[i]=0.0;
}*/
damul(u,a,nx,nx,nz,c);
damul(c,v,nx,nz,nz,d);
for(i=0;i<NX*NZ;I++) { for(i="0;i<=m-1;i++)" k="0;m=n-1;" int double } for(j="0;j<=k-1;j++)" * if(p[i] j="m-1;m=0;" while(k<m) d; m,k,j,i; p[]; n; n) p[],int mbbub(double void n—整型变量。待排序序列的长度。 p—双精度实型数组的指针。指示待排序序列（为顺序存储）的起始位置。 函数参数说明 return; c[u]="0.0;" for(l="0;l<=n-1;l++)" u="i*k+j;" i,j,l,u; c-为返回的乘积矩阵为m×k b-乘矩阵为n×k; a-被乘的矩阵为没m×n; c[]) k,double n,int m,int b[],int a[],double damul(double 矩阵相乘 delete c; v; u; jilu; a; fclose(fp2); fclose(fp1); fclose(fp); fwrite(&temp,sizeof(float),1,fp2); temp="d[i];">p[i+1])
{
d=p[i];
p[i]=p[i+1];
p[i+1]=d;
m=i;
}
j=k+1;
k=0;
for(i=m;i>=j;i--)
if(p[i-1]>p[i])
{
d=p[i];
p[i]=p[i-1];
p[i-1]=d;
k=i;
}
}
return;
}
//矩阵的奇异值分解
/*
函数参数说明
a—双精度实型二维数组，体积为m×n.存放m×n的实矩阵A；返回时其对角线给出奇异值
（已非递增次序排列），其余元素均为零。
m—整型变量。实矩阵A的行数。
n—整型变量。是矩阵A的列数。
u—双精度实型二维数组，体积为m×m.返回时存放左奇异向量。
v—双精度实型二维数组，体积为n×n.返回时存放右奇异向量。
eps—双精度实型变量。给定的精度要求。
ka—整型变量。其值为max(m,n)+1.
本函数的返回标志：若返回标志值小于0，则说明出现迭代60次还未完成求得某个奇异值
的情况，此时，矩阵A的分解式为UAV；若返回标志值大于0，则说明程序工作正常，二维
数组a中的对角线元素为奇异值。
*/
static void sss(double fg[2],double cs[2]);
static void ppp(double a[],double e[],double s[],double v[],int m,int n);
int bmuav(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=new double[ka];
e=new double[ka];
w=new double[ka];
it=60; k=n;
if (m-1<N) k="m-1;" (l if ll="k;" l="m;" (l<0) (n-2<m)>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++)
