fitting.cpp

c++的关于线性拟合的程序！ 采用了LU分解的方法来解方程组！

#include <iostream>
#include <cmath>
using namespace std;

void matric_multiple(double A[],double B[],int A1,int A2,int B1,int B2,double C[]);
void matric_transpose(double A[],int A1,int A2,double C[]);
void lubksb(double a[], int n, int indx[], double b[]);
void ludcmp(double a[], int n, int indx[], int d);

int main()
{
	double A[4][3]={1,2.5,3,4,5.3,6,7,8.8,9,10.4,11,12.5};
	double b[4]={2,4,5,9};
	int A1=4;
	int A2=3;
	double AA[12];
	int i,j,k=0,p=0;
	int indx[5];
	for(i=0;i<A1;i=i+1)
	{
		for(j=0;j<A2;j=j+1)
		{
			AA[k]=A[i][j];
			k=k+1;
		}
	}
	double AA_transpose[12];
	matric_transpose(AA,A1,A2,AA_transpose);
	double C[1000];
	double D[1000];
	matric_multiple(AA_transpose,AA,A2,A1,A1,A2,C);
	matric_multiple(AA_transpose,b,A2,A1,A1,1,D);
	ludcmp(C,A2,indx,p);
	lubksb(C,A2,indx,D);
	return 0;
}


void matric_multiple(double A[],double B[],int A1,int A2,int B1,int B2,double C[])
{
	// 本程序计算矩阵C=A*B，A1，A2为B1，B2分别为矩阵A，B的行数和列数，其中A2必须与B1相等。
	if(A2!=B1)
	{
		cout <<"number of columns of matric A must equal to the number of ranks of matric B";
		exit(0);
	}
	int i,j,k,q=0;
	double qq;
	for(i=0;i<A1;i=i+1)
	{
		for(k=0;k<B2;k=k+1)
		{
			C[q]=0;
			for(j=0;j<A2;j=j+1)
			{
				C[q]=C[q]+A[i*A2+j]*B[j*B2+k];
//				C[q]=C[q]+A[i*A2+j]*B[j*A1+k];
				qq=C[q];
			}
			q=q+1;
		}
	}
// 完毕
}

void matric_transpose(double A[],int A1,int A2,double C[])
{
// 本程序计算矩阵的转置
	int i,j,k=0;
	for(i=0;i<A2;i=i+1)
	{
		for(j=0;j<A1;j=j+1)
		{
			C[k]=A[i+j*A2];
			k=k+1;
		}
	}
//完毕
}


void ludcmp(double a[], int n, int indx[], int d)
{
	const int nmax=1000;
	const double tiny=1e-20;
	double vv[1000];
	double sum,dum,aamax;
	int i,j,k,imax;
	d=1;
	for (i=0; i<n; i++)
	{
		aamax=0.0;
		for (j=0; j<n; j++)
		{
			if (fabs(a[i*n+j])>aamax)
			{
				aamax=fabs(a[i*n+j]);
			}
		}
		if (aamax==0.0)
		{
			cout<<"singular matrix"<<endl;
		}
		vv[i]=1.0/aamax;
	}
	for (j=0; j<n; j++)
	{
		for (i=0; i<j; i++)
		{
            sum=a[i*n+j];
			for (int k=0; k<i; k++)
			{
				sum=sum-a[i*n+k]*a[k*n+j];
            }
            a[i*n+j]=sum;
		}
        aamax=0.0;
        for (i=j; i<n; i++)
		{
            sum=a[i*n+j];
			for (k=0; k<j; k++)
			{
				sum=sum-a[i*n+k]*a[k*n+j];
			}
            a[i*n+j]=sum;
            dum=vv[i]*fabs(sum);
            if (dum>=aamax)
			{
                imax=i;
                aamax=dum;
            }
        }
        if (j!=imax)
		{
            for (k=0; k<n; k++)
			{
                dum=a[imax*n+k];
                a[imax*n+k]=a[j*n+k];
                a[j*n+k]=dum;
            }
            d=-d;
            vv[imax]=vv[j];
        }
        indx[j]=imax;
		if(a[j*n+j]==0.0)
		{
			a[j*n+j]=tiny;
		}
        if (j!=n-1)
		{
            dum=1.0/a[j*n+j];
            for (i=j+1; i<n; i++)
			{
                a[i*n+j]=a[i*n+j]*dum;
            }
		}
    }
    if (a[n*n+n]==0.0)
	{
		a[n*n+n]=tiny;
	}
}

void lubksb(double a[], int n, int indx[], double b[])
{
	int i,j,ll;
	double sum;
    int ii = 0;
    for (i = 0; i<n; i++)
	{
        ll = indx[i];
        sum = b[ll];
        b[ll] = b[i];
        if (ii != 0 )
		{
            for (j = ii-1; j<i; j++)
			{
                sum = sum - a[i*n+j] * b[j];
            }
		}
        else if (sum != 0.0)
		{
			ii = i+1;
		}
        b[i] = sum;
    }
    for (i = n-1; i>=0; i--)
	{
        sum = b[i];
		for (j = i + 1; j<n; j++)
		{
            sum = sum - a[i*n+j] * b[j];
        }
        b[i] = sum / a[i*n+i];
    }
}