主题:[讨论]很奇怪这个程序能编译通过??
。。。。
int bmuav(a,m,n,u,v,eps,ka);
main()
{ int m,n,ka,i,j;
double a[5][4]={ {1.0,2.0,3.0,4.0},{6.0,7.0,8.0,9.0},{1.0,2.0,13.0,0.0},{16.0,17.0,8.0,9.0},{2.0,4.0,3.0,4.0}};
static double aa[4][5],c[5][4],u[5][5],v[4][4];
double eps;
m=5; n=4; ka=6; eps=0.000001;
i=bginv(a,m,n,aa,eps,u,v,ka);//aa为伪逆矩阵
if (i>0)
{ for (i=0; i<=3; i++)
{ for (j=0; j<=4; j++)
printf("%1.4f ",aa[i][j]);
printf("\n");
}
printf("\n");
}
/*printf("MAT A++ IS:\n");
i=bginv(aa,n,m,c,eps,v,u,ka);
。。。。
。。。。
printf("\n");
}*/
}
int bginv(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=bmuav(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);
}
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=(double*)malloc(ka*sizeof(double));
double *e=(double*)malloc(ka*sizeof(double));
double *w=(double*)malloc(ka*sizeof(double));
。。。
。。。。。
。。。。。。。
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; 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;
}
。。。
。。。。。
。。。。。
}
今天用到了矩阵求伪逆的算法,从c代码中找到了这个程序,编译通过了,本来如果直接用到程序中也就算了,但是偶然的扫了一眼,发现有个地方不明白了,不知道编译器怎么通过的,还是我学习了几年c一直是在错误的理解,特地把整个函数代码发上来和大家讨论一下:
注意开始的时候我是定义了一个5行4列的二维矩阵,double a[5][4];用过这个函数的都应该熟悉,里面的bginv()函数是用来求伪逆的过程,在这个函数中又用到了奇异值分解的函数bmuav();不再一一介绍他们了.问题出在我在函数bginv()中声明的参数为double a[];而在bmuav()中同样用到的是double a[];也就是将一个矩阵看成是一个数组的形式,以前我也这样用过.并且这里的所有子函数中都是用到的a[]形式,但是主函数中将a声明为一个二维矩阵,在子函数中姑且将它们看成是一个一维矩阵(一维矩阵和数组是一样的),那么这里编译的话应该出现指针不一致的情况,奇怪是通过的呢?很不明白...
另外,我自己编译了一个函数声明了一个double[][]后,用a[]的形式进行输出,会出现重大错误.高手指点一下怎么回事呢??
int bmuav(a,m,n,u,v,eps,ka);
main()
{ int m,n,ka,i,j;
double a[5][4]={ {1.0,2.0,3.0,4.0},{6.0,7.0,8.0,9.0},{1.0,2.0,13.0,0.0},{16.0,17.0,8.0,9.0},{2.0,4.0,3.0,4.0}};
static double aa[4][5],c[5][4],u[5][5],v[4][4];
double eps;
m=5; n=4; ka=6; eps=0.000001;
i=bginv(a,m,n,aa,eps,u,v,ka);//aa为伪逆矩阵
if (i>0)
{ for (i=0; i<=3; i++)
{ for (j=0; j<=4; j++)
printf("%1.4f ",aa[i][j]);
printf("\n");
}
printf("\n");
}
/*printf("MAT A++ IS:\n");
i=bginv(aa,n,m,c,eps,v,u,ka);
。。。。
。。。。
printf("\n");
}*/
}
int bginv(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=bmuav(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);
}
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=(double*)malloc(ka*sizeof(double));
double *e=(double*)malloc(ka*sizeof(double));
double *w=(double*)malloc(ka*sizeof(double));
。。。
。。。。。
。。。。。。。
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; 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;
}
。。。
。。。。。
。。。。。
}
今天用到了矩阵求伪逆的算法,从c代码中找到了这个程序,编译通过了,本来如果直接用到程序中也就算了,但是偶然的扫了一眼,发现有个地方不明白了,不知道编译器怎么通过的,还是我学习了几年c一直是在错误的理解,特地把整个函数代码发上来和大家讨论一下:
注意开始的时候我是定义了一个5行4列的二维矩阵,double a[5][4];用过这个函数的都应该熟悉,里面的bginv()函数是用来求伪逆的过程,在这个函数中又用到了奇异值分解的函数bmuav();不再一一介绍他们了.问题出在我在函数bginv()中声明的参数为double a[];而在bmuav()中同样用到的是double a[];也就是将一个矩阵看成是一个数组的形式,以前我也这样用过.并且这里的所有子函数中都是用到的a[]形式,但是主函数中将a声明为一个二维矩阵,在子函数中姑且将它们看成是一个一维矩阵(一维矩阵和数组是一样的),那么这里编译的话应该出现指针不一致的情况,奇怪是通过的呢?很不明白...
另外,我自己编译了一个函数声明了一个double[][]后,用a[]的形式进行输出,会出现重大错误.高手指点一下怎么回事呢??