ch1.h

数值处理算法程序

/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/*****  #include "CH1.h" 线性代数方程组的求解  *****/
#ifndef CH1_H_
#define CH1_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"

#include "CH2.h"
//*******************************************************************

int agaus(double a[],double b[],int n);//全选主元高斯消去法
int agjdn(double a[],double b[],int n,int m);//全选主元高斯约当消去法
int acgas(double ar[],double ai[],int n,double br[],double bi[]);//复系数方程组全选主元高斯消去法
int acjdn(double ar[],double ai[],double br[],double bi[],int n,int m);//复系数方程组全选主元高斯约当消去法
int atrde(double b[],int n,int m,double d[]);//三对角线方程的追赶法
int aband(double b[],double d[],int n,int l,int il,int m);//一般带型方程的求解--**例子程序有溢出,解的值就为零!**
int aldle(double a[],int n,int m,double c[]);//求解对称方程组的分解法
int achol(double a[],int n,int m,double d[]);//求解对称正定方程组的平方根法Cholesky
int aggje(double a[],int n,double b[]);//求解大型稀疏方程组的全选主元高斯约当消去法
int atlvs(double t[],int n,double b[],double x[]);//求解托伯利兹方程组的列文方法
int agsdl(double a[],double b[],int n,double x[],double eps);//高斯-赛德尔迭代法
void agrad(double a[],int n,double b[],double eps,double x[]);//求解对称正定方程组的共轭梯度法
int agmqr(double a[],int m,int n,double b[],double q[]);//求解最小二乘问题的豪斯荷尔德变换法
int agmiv(double a[],int m,int n,double b[],double x[],double aa[],double eps,double u[],double v[],int ka);//求解最小二乘问题的广义逆法
int abint(double a[],int n,double b[],double eps,double x[]);//病态方程组的求解

//*******************************************************************
//int agaus(double* a,double* b,int n)
int agaus(double a[],double b[],int n)
{ 
	int *js,l,k,i,j,is,p,q;
    double d,t;
    js=(int*)malloc(n*sizeof(int));
    l=1;
    for (k=0;k<=n-2;k++)
      { d=0.0;
        for (i=k;i<=n-1;i++)
          for (j=k;j<=n-1;j++)
            { t=fabs(a[i*n+j]);
              if (t>d) { d=t; js[k]=j; is=i;}
            }
        if (d+1.0==1.0) l=0;
        else
          { if (js[k]!=k)
              for (i=0;i<=n-1;i++)
                { p=i*n+k; q=i*n+js[k];
                  t=a[p]; a[p]=a[q]; a[q]=t;
                }
            if (is!=k)
              { for (j=k;j<=n-1;j++)
                  { p=k*n+j; q=is*n+j;
                    t=a[p]; a[p]=a[q]; a[q]=t;
                  }
                t=b[k]; b[k]=b[is]; b[is]=t;
              }
          }
        if (l==0)
          { free(js); //printf("fail\n");
            return(0);
          }
        d=a[k*n+k];
        for (j=k+1;j<=n-1;j++)
          { p=k*n+j; a[p]=a[p]/d;}
        b[k]=b[k]/d;
        for (i=k+1;i<=n-1;i++)
          { for (j=k+1;j<=n-1;j++)
              { p=i*n+j;
                a[p]=a[p]-a[i*n+k]*a[k*n+j];
              }
            b[i]=b[i]-a[i*n+k]*b[k];
          }
      }
    d=a[(n-1)*n+n-1];
    if (fabs(d)+1.0==1.0)
      { free(js); //printf("fail\n");
        return(0);
      }
    b[n-1]=b[n-1]/d;
    for (i=n-2;i>=0;i--)
      { t=0.0;
        for (j=i+1;j<=n-1;j++)
          t=t+a[i*n+j]*b[j];
        b[i]=b[i]-t;
      }
    js[n-1]=n-1;
    for (k=n-1;k>=0;k--)
      if (js[k]!=k)
        { t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}
    free(js);
    return(1);
 }

/////////////////////////////////////////////////////////////////

int agjdn(double a[],double b[],int n,int m)
{ int *js,l,k,i,j,is,p,q;
  double d,t;
  js=(int*)malloc(n*sizeof(int));
  l=1;
  for (k=0;k<=n-1;k++)
    { d=0.0;
      for (i=k;i<=n-1;i++)
        for (j=k;j<=n-1;j++)
          { t=fabs(a[i*n+j]);
            if (t>d) { d=t; js[k]=j; is=i;}
          }
      if (d+1.0==1.0) l=0;
      else
        { if (js[k]!=k)
            for (i=0;i<=n-1;i++)
              { p=i*n+k; q=i*n+js[k];
                t=a[p]; a[p]=a[q]; a[q]=t;
              }
          if (is!=k)
            { for (j=k;j<=n-1;j++)
                { p=k*n+j; q=is*n+j;
                  t=a[p]; a[p]=a[q]; a[q]=t;
                }
              for (j=0;j<=m-1;j++)
                { p=k*m+j; q=is*m+j;
                  t=b[p]; b[p]=b[q]; b[q]=t;
                }
            }
        }
      if (l==0)
        { free(js); //printf("fail\n");
          return(0);
        }
      d=a[k*n+k];
      for (j=k+1;j<=n-1;j++)
        { p=k*n+j; a[p]=a[p]/d;}
      for (j=0;j<=m-1;j++)
        { p=k*m+j; b[p]=b[p]/d;}
      for (j=k+1;j<=n-1;j++)
        for (i=0;i<=n-1;i++)
          { p=i*n+j;
            if (i!=k)
              a[p]=a[p]-a[i*n+k]*a[k*n+j];
          }
      for (j=0;j<=m-1;j++)
      for (i=0;i<=n-1;i++)
        { p=i*m+j;
          if (i!=k)
            b[p]=b[p]-a[i*n+k]*b[k*m+j];
        }
    }
  for (k=n-1;k>=0;k--)
    if (js[k]!=k)
      for (j=0;j<=m-1;j++)
        { p=k*m+j; q=js[k]*m+j;
          t=b[p]; b[p]=b[q]; b[q]=t;
        }
  free(js);
  return(1);
}
///////////////////////acgas(&ar[0][0],&ai[0][0],4,&br[0],&bi[0])
int acgas(double ar[],double ai[],int n,double br[],double bi[])
{
	int *js,l,k,i,j,is,u,v;
	double p,q,s,d;
	js=(int*)malloc(n*sizeof(int));
	for (k=0;k<=n-2;k++)
    { d=0.0;
      for (i=k;i<=n-1;i++)
      for (j=k;j<=n-1;j++)
      {
		  u=i*n+j;
          p=ar[u]*ar[u]+ai[u]*ai[u];
          if (p>d) {d=p;js[k]=j;is=i;}
      }
      if (d+1.0==1.0)
      { 
		  free(js); //printf("err**fail\n");
          return(0);
      }
      if (is!=k)
      { 
		  for (j=k;j<=n-1;j++)
          { 
			  u=k*n+j; v=is*n+j;
              p=ar[u]; ar[u]=ar[v]; ar[v]=p;
              p=ai[u]; ai[u]=ai[v]; ai[v]=p;
		  }
          p=br[k]; br[k]=br[is]; br[is]=p;
          p=bi[k]; bi[k]=bi[is]; bi[is]=p;
       }
      if (js[k]!=k)
        for (i=0;i<=n-1;i++)
        { 
			u=i*n+k; v=i*n+js[k];
            p=ar[u]; ar[u]=ar[v]; ar[v]=p;
            p=ai[u]; ai[u]=ai[v]; ai[v]=p;
        }
      v=k*n+k;
      for (j=k+1;j<=n-1;j++)
      { 
		  u=k*n+j;
          p=ar[u]*ar[v]; q=-ai[u]*ai[v];
          s=(ar[v]-ai[v])*(ar[u]+ai[u]);
          ar[u]=(p-q)/d; ai[u]=(s-p-q)/d;
      }
      p=br[k]*ar[v]; q=-bi[k]*ai[v];
      s=(ar[v]-ai[v])*(br[k]+bi[k]);
      br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
      for (i=k+1;i<=n-1;i++)
      { 
		  u=i*n+k;
          for (j=k+1;j<=n-1;j++)
          { 
			  v=k*n+j; l=i*n+j;
              p=ar[u]*ar[v]; q=ai[u]*ai[v];
              s=(ar[u]+ai[u])*(ar[v]+ai[v]);
              ar[l]=ar[l]-p+q;
              ai[l]=ai[l]-s+p+q;
          }
		  p=ar[u]*br[k]; q=ai[u]*bi[k];
          s=(ar[u]+ai[u])*(br[k]+bi[k]);
          br[i]=br[i]-p+q; bi[i]=bi[i]-s+p+q;
	  }
  }
  u=(n-1)*n+n-1;
  d=ar[u]*ar[u]+ai[u]*ai[u];
  if (d+1.0==1.0)
  { 
	  free(js); //printf("err**fail\n");
      return(0);
  }
  p=ar[u]*br[n-1]; q=-ai[u]*bi[n-1];
  s=(ar[u]-ai[u])*(br[n-1]+bi[n-1]);
  br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
  for (i=n-2;i>=0;i--)
  for (j=i+1;j<=n-1;j++)
  { 
	  u=i*n+j;
      p=ar[u]*br[j]; q=ai[u]*bi[j];
      s=(ar[u]+ai[u])*(br[j]+bi[j]);
      br[i]=br[i]-p+q;
      bi[i]=bi[i]-s+p+q;
  }
  js[n-1]=n-1;
  for (k=n-1;k>=0;k--)
  if (js[k]!=k)
  {
	  p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
	  p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
  }
  free(js);
  return(1);
}
/////////////if (acjdn(&ar[0][0],&ai[0][0],&br[0][0],&bi[0][0],4,2)!=0)
int acjdn(double ar[],double ai[],double br[],double bi[],int n,int m)
{
	int *js,l,k,i,j,is,u,v;
	double p,q,s,d;
	js=(int*)malloc(n*sizeof(int));
	for (k=0;k<=n-1;k++)
    { 
		d=0.0;
		for (i=k;i<=n-1;i++)
			for (j=k;j<=n-1;j++)
			{ 
				u=i*n+j;
				p=ar[u]*ar[u]+ai[u]*ai[u];
				if (p>d) 
				{
					d=p;js[k]=j;is=i;
				}
			}
		if (d+1.0==1.0)
        {
			free(js); printf("err**fail\n");
			return(0);
        }
		if (is!=k)
        {
			for (j=k;j<=n-1;j++)
            {
				u=k*n+j; v=is*n+j;p=ar[u];
				ar[u]=ar[v]; ar[v]=p; p=ai[u];
				ai[u]=ai[v]; ai[v]=p;
            }
          for (j=0;j<=m-1;j++)
            { u=k*m+j; v=is*m+j;
              p=br[u]; br[u]=br[v]; br[v]=p;
              p=bi[u]; bi[u]=bi[v]; bi[v]=p;
            }
        }
      if (js[k]!=k)
        for (i=0;i<=n-1;i++)
          { u=i*n+k; v=i*n+js[k];
            p=ar[u]; ar[u]=ar[v]; ar[v]=p;
            p=ai[u]; ai[u]=ai[v]; ai[v]=p;
          }
      v=k*n+k;
      for (j=k+1;j<=n-1;j++)
        { u=k*n+j;
          p=ar[u]*ar[v]; q=-ai[u]*ai[v];
          s=(ar[v]-ai[v])*(ar[u]+ai[u]);
          ar[u]=(p-q)/d; ai[u]=(s-p-q)/d;
        }
      for (j=0;j<=m-1;j++)
        { u=k*m+j;
          p=br[u]*ar[v]; q=-bi[u]*ai[v];
          s=(ar[v]-ai[v])*(br[u]+bi[u]);
          br[u]=(p-q)/d; bi[u]=(s-p-q)/d;
        }
      for (i=0;i<=n-1;i++)
        if (i!=k)
          { u=i*n+k;
            for (j=k+1;j<=n-1;j++)
              { v=k*n+j; l=i*n+j;
                p=ar[u]*ar[v]; q=ai[u]*ai[v];
                s=(ar[u]+ai[u])*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
              }
            for (j=0;j<=m-1;j++)
              { l=i*m+j; v=k*m+j;
                p=ar[u]*br[v]; q=ai[u]*bi[v];
                s=(ar[u]+ai[u])*(br[v]+bi[v]);
                br[l]=br[l]-p+q; bi[l]=bi[l]-s+p+q;
              }
          }
    }
  for (k=n-1;k>=0;k--)
    if (js[k]!=k)
      for (j=0;j<=m-1;j++)
        { u=k*m+j;v=js[k]*m+j;
          p=br[u]; br[u]=br[v]; br[v]=p;
          p=bi[u]; bi[u]=bi[v]; bi[v]=p;
        }
  free(js);
  return(1);
}
///////////////////// if (atrde(b,5,13,d)>0)
int atrde(double b[],int n,int m,double d[])
{ int k,j;
  double s;
  if (m!=(3*n-2))
    { printf("err\n"); return(-2);}
  for (k=0;k<=n-2;k++)
    { j=3*k; s=b[j];
      if (fabs(s)+1.0==1.0)
        { printf("fail\n"); return(0);}
      b[j+1]=b[j+1]/s;
      d[k]=d[k]/s;
      b[j+3]=b[j+3]-b[j+2]*b[j+1];
      d[k+1]=d[k+1]-b[j+2]*d[k];
    }
  s=b[3*n-3];
  if (fabs(s)+1.0==1.0)
    { printf("fail\n"); return(0);}
  d[n-1]=d[n-1]/s;
  for (k=n-2;k>=0;k--)
    d[k]=d[k]-b[3*k+1]*d[k+1];
  return(2);
}
/////////////////////////////////////////////////////////////
int aband(double b[],double d[],int n,int l,int il,int m)
{ int ls,k,i,j,is,u,v;
  double p,t;
  if (il!=(2*l+1))
    { printf("fail\n"); return(-2);}
  ls=l;
  for (k=0;k<=n-2;k++)
    { p=0.0;
      for (i=k;i<=ls;i++)
        { t=fabs(b[i*il]);
          if (t>p) {p=t; is=i;}
        }
      if (p+1.0==1.0)
        { printf("fail\n"); return(0);}
      for (j=0;j<=m-1;j++)
        { u=k*m+j; v=is*m+j;
          t=d[u]; d[u]=d[v]; d[v]=t;
        }
      for (j=0;j<=il-1;j++)
        { u=k*il+j; v=is*il+j;
          t=b[u]; b[u]=b[v]; b[v]=t;
        }
      for (j=0;j<=m-1;j++)
        { u=k*m+j; d[u]=d[u]/b[k*il];}
      for (j=1;j<=il-1;j++)
        { u=k*il+j; b[u]=b[u]/b[k*il];}
      for (i=k+1;i<=ls;i++)