《C常用算法程序集-徐士良》源代码

#include "stdlib.h" #include "math.h" int bmuav(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; for (i=kk;