#include "stdlib.h"
#include "math.h"
#include "stdio.h"
int muav(a,m,n,u,v,eps,ka)
int m,n,ka;
double eps,a[],u[],v[];
{ 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;
void ppp();
void sss();
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;
评论0