inverssolve.cpp

用vc++实现广义的最小二乘法程序

  #include "math.h"
  #include "stdio.h"
  #include "stdlib.h"
 #include "iostream.h"
#include "malloc.h"
 /* #include "agmiv.c"
  #include "bginv.c"
  #include "bmuav.c"
  #include "dngin.c"*/
//方程左端的值
void least_dnginf( int m, int n,double *x,double *d);
//雅可比矩阵
void least_dngins(int m,int n,double *x,double *p);
///奇异值分解
int least_dngin(int m,int n,double eps1,double eps2,double *x,int ka);
//最小二乘分解
int least_agmiv(double *a,int m,int n,double *b,double *x,double *aa,double eps,double *u,double *v,int ka);

//奇异值分解中是否满足精度要求
int least_bginv(double *a,int m,int n,double *aa,double eps,double *u,double *v,int ka);
//奇异值分解
int least_bmuav(double *a,int m,int n,double *u,double *v,double eps,int ka);
void least_sss(double *fg,double *cs);
void least_ppp(double *a,double *e,double *s,double *v,int m,int n);
  main()
  { int m,n,i,ka;
    double eps1,eps2;
     double x[12]={0.8,-0,0,0,0,1,0,0,0,0,1,0.6};
    m=24; n=12; ka=4; eps1=0.000001; eps2=0.000001;
    i=least_dngin(m,n,eps1,eps2,x,ka);
    printf("\n");
    printf("i=%d\n",i);
    printf("\n");
    for (i=0; i<=12; i++)
      printf("x(%d)=%13.7e\n",i,x[i]);
    printf("\n");
  }
//方程左端的值
/*  void least_dnginf( int m, int n,double *x,double *d)
 // int m,n;
  //double x[],d[];
  { //m=m; n=n;
    d[0]=x[0]*x[0]+7.0*x[0]*x[1]+3.0*x[1]*x[1]+0.5;
    d[1]=x[0]*x[0]-2.0*x[0]*x[1]+x[1]*x[1]-1.0;
    d[2]=x[0]+x[1]+1.0;
    //return;
  }*/
 //雅可比矩阵

 /*void least_dngins(int m,int n,double *x,double *p)
 // int m,n;
 // double x[2],p[3][2];
  { //m=m; n=n;
    p[0*n+0]=2.0*x[0]+7.0*x[1];
    p[0*n+1]=7.0*x[0]+6.0*x[1];
    p[1*n+0]=2.0*x[0]-2.0*x[1];
    p[1*n+1]=-2.0*x[0]+2.0*x[1];
    p[2*n+0]=1.0;
    p[2*n+1]=1.0;
    //return;
  }*/
  ///////////////////////
  ///奇异值分解
    int least_dngin(int m,int n,double eps1,double eps2,double *x,int ka)
  //int m,n,ka;
  //double eps1,eps2,x[];
  {// extern void dnginf();
   // extern void dngins();
   // extern int agmiv();
    int i,j,k,l,kk,jt;
    double y[10],b[10],alpha,z,h2,y1,y2,y3,y0,h1;
   // double *p,*d,*pp,*dx,*u,*v,*w;
    double  p[4000],d[4000],pp[1000],dx[1000],u[1000],v[1000],w[100];
    /*p=malloc(m*n*sizeof(double));
    d=malloc(m*sizeof(double));
    pp=malloc(n*m*sizeof(double));
    dx=malloc(n*sizeof(double));
    u=malloc(m*m*sizeof(double));
    v=malloc(n*n*sizeof(double));
    w=malloc(ka*sizeof(double));*/
  /* p=new double[m*n];
   d=new double[m];
   pp=new double[n*m];
   dx=new double(n);
   u=new double[m*m];
   v=new double[n*n];
   w=new double[ka];*/
    l=60; alpha=1.0;
    while (l>0)
      { least_dnginf(m,n,x,d);
        least_dngins(m,n,x,p);
        jt=least_agmiv(p,m,n,d,dx,pp,eps2,u,v,ka);
        if (jt<0)
          { //free(p); free(d); free(pp); free(w);
            //free(dx); free(u); free(v); return(jt);
		/*	delete[]p;delete[]d;delete[]pp;
			delete[]dx;delete[] u;delete[]v;
			delete []w;*/
			return(jt);
          }
        j=0; jt=1; h2=0.0;
        while (jt==1)
          { jt=0;
            if (j<=2) z=alpha+0.01*j;
            else z=h2;
            for (i=0; i<=n-1; i++) w[i]=x[i]-z*dx[i];
            least_dnginf(m,n,w,d);
            y1=0.0;
            for (i=0; i<=m-1; i++) y1=y1+d[i]*d[i];
            for (i=0; i<=n-1; i++)
              w[i]=x[i]-(z+0.00001)*dx[i];
           least_dnginf(m,n,w,d);
            y2=0.0;
            for (i=0; i<=m-1; i++) y2=y2+d[i]*d[i];
            y0=(y2-y1)/0.00001;
            if (fabs(y0)>1.0e-10)
              { h1=y0; h2=z;
                if (j==0) { y[0]=h1; b[0]=h2;}
                else
                  { y[j]=h1; kk=0; k=0;
                    while ((kk==0)&&(k<=j-1))
                        { y3=h2-b[k];
                          if (fabs(y3)+1.0==1.0) kk=1;
                          else h2=(h1-y[k])/y3;
                          k=k+1;
                        }
                    b[j]=h2;
                    if (kk!=0) b[j]=1.0e+35;
                    h2=0.0;
                    for (k=j-1; k>=0; k--)
                      h2=-y[k]/(b[k+1]+h2);
                    h2=h2+b[0];
                  }
                j=j+1;
                if (j<=7) jt=1;
                else z=h2;
              }
          }
        alpha=z; y1=0.0; y2=0.0;
        for (i=0; i<=n-1; i++)
          { dx[i]=-alpha*dx[i];
            x[i]=x[i]+dx[i];
            y1=y1+fabs(dx[i]);
            y2=y2+fabs(x[i]);
          }
        if (y1<eps1*y2)
          { //free(p); free(pp); free(d); free(w);
            //free(dx); free(u); free(v); return(1);
		/*	delete [] p;delete [] pp;delete [] d;//delete[] pp;
			delete [] w;delete [] dx;delete [] u;
			delete [] v;*/
			return(1);
          }
        l=l-1;
      }
    //free(p); free(pp); free(d); free(dx);
   // free(u); free(v); free(w); return(0);
	/*delete [] p; delete[] pp;delete [] d;
	delete [] dx;delete[] u;delete [] v;
	delete [] w;*/
	return(0);
  }
  //////广义逆求解最小线性二乘问题
int least_agmiv(double *a,int m,int n,double *b,double *x,double *aa,double eps,double *u,double *v,int ka)
  //int m,n,ka;
 // double a[],b[],x[],aa[],u[],v[],eps;
  { int i,j;
    //extern int bginv();
    i=least_bginv(a,m,n,aa,eps,u,v,ka);
    if (i<0) return(-1);
    for (i=0; i<=n-1; i++)
      { x[i]=0.0;
        for (j=0; j<=m-1; j++)
          x[i]=x[i]+aa[i*m+j]*b[j];
      }
    return(1);
  }
  ///////////广义逆的函数
    int least_bginv(double *a,int m,int n,double *aa,double eps,double *u,double *v,int ka)
  //int m,n,ka;
 // double a[],aa[],u[],v[],eps;
  { int i,j,k,l,t,p,q,f;
  //  extern int bmuav();
    i=least_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 least_bmuav(double *a,int m,int n,double *u,double *v,double eps,int 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,*cs;
    //double *s,*e,*w;
    //void ppp();
    //void sss();
    /*s=malloc(ka*sizeof(double));
    e=malloc(ka*sizeof(double));
    w=malloc(ka*sizeof(double));*/
	/*fg= new double[2];
	cs=new  double[2];
	s=new double[ka];
	e=new double[ka];
	w=new double[ka];*/
	double cs[2],fg[2],w[50],s[40],e[50];
    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)
          {
		least_ppp(a,e,s,v,m,n);
            //free(s); free(e); free(w); return(1);
	     //  delete [] s;delete [] e;delete[] w;
		   return(1);

          }
        if (it==0)
          { 
			least_ppp(a,e,s,v,m,n);
            //free(s); free(e); free(w); return(-1);
	//	delete []s;delete []e;delete[]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)