7.10 求解非线性方程组最小二乘解的广义逆法 ngin.c

许士良常用算法程序集C语言,包括c++一些常用算法代码

#include "math.h"
#include "stdlib.h"
#include "6gmiv.c"
int ngin(m,n,eps1,eps2,x,ka,f,s)
void (*f)(),(*s)();
int m,n,ka;
double eps1,eps2,x[];
{ 
	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;
    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));
    l=60; alpha=1.0;
    while (l>0)
    { 
		(*f)(m,n,x,d);
        (*s)(m,n,x,p);
        jt=gmiv(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);
        }
        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];
            (*f)(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];
            (*f)(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);
        }
        l=l-1;
    }
    free(p); free(pp); free(d); free(dx);
    free(u); free(v); free(w); 
	return(0);
}