ch2.h

该文件包含《线性代数方程组的求解》、《矩阵运算》、《矩阵特征值与特征向量的计算》、《数值积分》等算法与C程序

/************************************************
 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 "CH2.h" 矩阵运算  *****/
#ifndef CH2_H_
#define CH2_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"
//*******************************************************************
void brmul(double a[],double b[],int m,int n,int k,double* c);//实矩阵相乘
void bcmul(double ar[],double ai[],double br[],double bi[],int m,int n,int k,double cr[],double ci[]);//复矩阵相乘
int brinv(double a[],int n);//实矩阵求逆的全选主元高斯约当法
int bcinv(double ar[],double ai[],int n);//复矩阵求逆的全选主元高斯约当法
int bssgj(double a[],int n);//对称正定矩阵的求逆
int btrch(double t[],double tt[],int n,double b[]);//托伯利兹矩阵求逆的特兰持方法
double bsdet(double a[],int n);//求行列式值的全选主元高斯消去法
int brank(double a[],int m,int n);//求矩阵秩的全选主元高斯消去法
int bchol(double a[],int n,double* det);//对称正定矩阵的乔里斯基分解与行列式的求值
int blluu(double a[],int n,double l[],double u[]);//矩阵的三角分解
int bmaqr(double a[],int m,int n,double q[]);//一般实矩阵的QR分解
int bmuav(double a[],int m,int n,double u[],double v[],double eps,int ka);//一般实矩阵的奇异值分解
int bginv(double a[],int m,int n,double aa[],double eps,double u[],double v[],int ka);//求广义逆的奇异值分解法
//*******************************************************************
void brmul(double a[],double b[],int m,int n,int k,double c[])
{ 
	int i,j,l,u;
    for (i=0; i<=m-1; i++)
		for (j=0; j<=k-1; j++)
		{
			u=i*k+j; c[u]=0.0;
			for (l=0; l<=n-1; l++)
			  c[u]=c[u]+a[i*n+l]*b[l*k+j];
		 }
    return;
}
/////////////////////////////////////////////////////////////
void bcmul(double ar[],double ai[],double br[],double bi[],int m,int n,int k,double cr[],double ci[])
{ int i,j,l,u,v,w;
  double p,q,s;
  for (i=0; i<=m-1; i++)
  for (j=0; j<=k-1; j++)
    { u=i*k+j;
      cr[u]=0.0; ci[u]=0.0;
      for (l=0; l<=n-1; l++)
        { v=i*n+l; w=l*k+j;
          p=ar[v]*br[w];
          q=ai[v]*bi[w];
          s=(ar[v]+ai[v])*(br[w]+bi[w]);
          cr[u]=cr[u]+p-q;
          ci[u]=ci[u]+s-p-q;
        }
    }
  return;
}
/////////////////////////////////////////////////////////////
int brinv(double a[],int n)
{ 
	int *is,*js,i,j,k,l,u,v;
	double d,p;
	is=(int*)malloc(n*sizeof(int));
	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++)
			{ 
				l=i*n+j; p=fabs(a[l]);
				if (p>d)
				{
					d=p; is[k]=i; js[k]=j;
				}
			}      
			if (d+1.0==1.0)
			{ 
				free(is); free(js); /*failure*/
				return(0);
			}
			if(is[k]!=k)
				for (j=0; j<=n-1; j++)
				{ 
					u=k*n+j; v=is[k]*n+j;
					p=a[u]; a[u]=a[v]; a[v]=p;
				}
				if(js[k]!=k)
					for (i=0; i<=n-1; i++)
					{ 
						u=i*n+k; v=i*n+js[k];
						p=a[u]; a[u]=a[v]; a[v]=p;
					}
					l=k*n+k;
					a[l]=1.0/a[l];
					for (j=0; j<=n-1; j++)
						if (j!=k)
						{ u=k*n+j; a[u]=a[u]*a[l];}
						for(i=0; i<=n-1; i++)
							if(i!=k)
								for(j=0; j<=n-1; j++)
									if(j!=k)
									{ 
										u=i*n+j;
										a[u]=a[u]-a[i*n+k]*a[k*n+j];
									}
									for(i=0; i<=n-1; i++)
										if(i!=k)
										{ u=i*n+k; a[u]=-a[u]*a[l];}
    }
	for (k=n-1; k>=0; k--)
	{ 
		if(js[k]!=k)
			for(j=0; j<=n-1; j++)
			{ 
				u=k*n+j; v=js[k]*n+j;
				p=a[u]; a[u]=a[v]; a[v]=p;
			}
			if(is[k]!=k)
				for(i=0; i<=n-1; i++)
				{
					u=i*n+k; v=i*n+is[k];
					p=a[u]; a[u]=a[v]; a[v]=p;
				}
    }
	free(is); free(js);
	return(1);
}
/////////////////////////////////////////////////////////////
int bcinv(double ar[],double ai[],int n)
{ int *is,*js,i,j,k,l,u,v,w;
  double p,q,s,t,d,b;
  is=(int*)malloc(n*sizeof(int));
  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; is[k]=i; js[k]=j;}
        }
      if (d+1.0==1.0)
        { free(is); free(js); printf("err**not inv\n");
          return(0);
        }
      if (is[k]!=k)
        for (j=0; j<=n-1; j++)
          { u=k*n+j; v=is[k]*n+j;
            t=ar[u]; ar[u]=ar[v]; ar[v]=t;
            t=ai[u]; ai[u]=ai[v]; ai[v]=t;
          }
      if (js[k]!=k)
        for (i=0; i<=n-1; i++)
          { u=i*n+k; v=i*n+js[k];
            t=ar[u]; ar[u]=ar[v]; ar[v]=t;
            t=ai[u]; ai[u]=ai[v]; ai[v]=t;
          }
      l=k*n+k;
      ar[l]=ar[l]/d; ai[l]=-ai[l]/d;
      for (j=0; j<=n-1; j++)
        if (j!=k)
          { u=k*n+j;
            p=ar[u]*ar[l]; q=ai[u]*ai[l];
            s=(ar[u]+ai[u])*(ar[l]+ai[l]);
            ar[u]=p-q; ai[u]=s-p-q;
          }
      for (i=0; i<=n-1; i++)
        if (i!=k)
          { v=i*n+k;
            for (j=0; j<=n-1; j++)
              if (j!=k)
                { u=k*n+j;  w=i*n+j;
                  p=ar[u]*ar[v]; q=ai[u]*ai[v];
                  s=(ar[u]+ai[u])*(ar[v]+ai[v]);
                  t=p-q; b=s-p-q;
                  ar[w]=ar[w]-t;
                  ai[w]=ai[w]-b;
                }
          }
      for (i=0; i<=n-1; i++)
        if (i!=k)
          { u=i*n+k;
            p=ar[u]*ar[l]; q=ai[u]*ai[l];
            s=(ar[u]+ai[u])*(ar[l]+ai[l]);
            ar[u]=q-p; ai[u]=p+q-s;
          }
    }
  for (k=n-1; k>=0; k--)
    { if (js[k]!=k)
        for (j=0; j<=n-1; j++)
          { u=k*n+j; v=js[k]*n+j;
            t=ar[u]; ar[u]=ar[v]; ar[v]=t;
            t=ai[u]; ai[u]=ai[v]; ai[v]=t;
          }
      if (is[k]!=k)
        for (i=0; i<=n-1; i++)
          { u=i*n+k; v=i*n+is[k];
            t=ar[u]; ar[u]=ar[v]; ar[v]=t;
            t=ai[u]; ai[u]=ai[v]; ai[v]=t;
          }
    }
  free(is); free(js);
  return(1);
}
/////////////////////////////////////////////////////////////
int bssgj(double a[],int n)
{ int i,j,k,m;
  double w,g,*b;
  b=(double*)malloc(n*sizeof(double));
  for (k=0; k<=n-1; k++)
    { w=a[0];
      if (fabs(w)+1.0==1.0)
        { free(b); printf("fail\n"); return(-2);}
      m=n-k-1;
      for (i=1; i<=n-1; i++)
        { g=a[i*n]; b[i]=g/w;
          if (i<=m) b[i]=-b[i];
          for (j=1; j<=i; j++)
            a[(i-1)*n+j-1]=a[i*n+j]+g*b[j];
        }
      a[n*n-1]=1.0/w;
      for (i=1; i<=n-1; i++)
        a[(n-1)*n+i-1]=b[i];
    }
  for (i=0; i<=n-2; i++)
  for (j=i+1; j<=n-1; j++)
    a[i*n+j]=a[j*n+i];
  free(b);
  return(2);
}
/////////////////////////////////////////////////////////////
int btrch(double t[],double tt[],int n,double b[])
{ int i,j,k;
  double a,s,*c,*r,*p;
  c=(double*)malloc(n*sizeof(double));
  r=(double*)malloc(n*sizeof(double));
  p=(double*)malloc(n*sizeof(double));
  if (fabs(t[0])+1.0==1.0)
    { free(c); free(r); free(p);
      printf("fail\n"); return(-1);
    }
  a=t[0]; c[0]=tt[1]/t[0]; r[0]=t[1]/t[0];
  for (k=0; k<=n-3; k++)
    { s=0.0;
      for (j=1; j<=k+1; j++)
        s=s+c[k+1-j]*tt[j];
      s=(s-tt[k+2])/a;
	for (i=0; i<=k; i++)
        p[i]=c[i]+s*r[k-i];
      c[k+1]=-s;
      s=0.0;
      for (j=1; j<=k+1; j++)
        s=s+r[k+1-j]*t[j];
      s=(s-t[k+2])/a;
      for (i=0; i<=k; i++)
        { r[i]=r[i]+s*c[k-i];
          c[k-i]=p[k-i];
        }
      r[k+1]=-s;
	a=0.0;
      for (j=1; j<=k+2; j++)
        a=a+t[j]*c[j-1];
      a=t[0]-a;
      if (fabs(a)+1.0==1.0)
        { free(c); free(r); free(p);
          printf("fail\n"); return(-1);
        }
    }
  b[0]=1.0/a;
  for (i=0; i<=n-2; i++)
    { k=i+1; j=(i+1)*n;
      b[k]=-r[i]/a; b[j]=-c[i]/a;
    }
  for (i=0; i<=n-1; i++)
  for (j=0; j<=n-2; j++)
    { k=(i+1)*n+j+1;
      b[k]=b[i*n+j]-c[i]*b[j+1];
      b[k]=b[k]+c[n-j-2]*b[n-i-1];
    }
  free(c); free(r); free(p);
  return(1);
}
/////////////////////////////////////////////////////////////
double bsdet(double a[],int n)
{ int i,j,k,is,js,l,u,v;
  double f,det,q,d;
  f=1.0; det=1.0;
  for (k=0; k<=n-2; k++)
    { q=0.0;
      for (i=k; i<=n-1; i++)
      for (j=k; j<=n-1; j++)
        { l=i*n+j; d=fabs(a[l]);
	    if (d>q) { q=d; is=i; js=j;}
        }
      if (q+1.0==1.0)
        { det=0.0; return(det);}
      if (is!=k)
        { f=-f;
          for (j=k; j<=n-1; j++)
            { u=k*n+j; v=is*n+j;
              d=a[u]; a[u]=a[v]; a[v]=d;
            }
        }
      if (js!=k)
        { f=-f;
          for (i=k; i<=n-1; i++)
            { u=i*n+js; v=i*n+k;
              d=a[u]; a[u]=a[v]; a[v]=d;
            }
        }
      l=k*n+k;
      det=det*a[l];
      for (i=k+1; i<=n-1; i++)
        { d=a[i*n+k]/a[l];
          for (j=k+1; j<=n-1; j++)
            { u=i*n+j;
              a[u]=a[u]-d*a[k*n+j];
            }
        }
    }
  det=f*det*a[n*n-1];
  return(det);
}
/////////////////////////////////////////////////////////////
int brank(double a[],int m,int n)
{ int i,j,k,nn,is,js,l,ll,u,v;
  double q,d;
  nn=m;
  if (m>=n) nn=n;
  k=0;
  for (l=0; l<=nn-1; l++)
    { q=0.0;
      for (i=l; i<=m-1; i++)
      for (j=l; j<=n-1; j++)
        { ll=i*n+j; d=fabs(a[ll]);
	    if (d>q) { q=d; is=i; js=j;}
        }
      if (q+1.0==1.0) return(k);
      k=k+1;
      if (is!=l)
        { for (j=l; j<=n-1; j++)
            { u=l*n+j; v=is*n+j;
              d=a[u]; a[u]=a[v]; a[v]=d;
            }
        }
      if (js!=l)
        { for (i=l; i<=m-1; i++)
            { u=i*n+js; v=i*n+l;
              d=a[u]; a[u]=a[v]; a[v]=d;
            }
        }
      ll=l*n+l;
      for (i=l+1; i<=n-1; i++)
        { d=a[i*n+l]/a[ll];
          for (j=l+1; j<=n-1; j++)
            { u=i*n+j;
              a[u]=a[u]-d*a[l*n+j];
            }
        }
    }
  return(k);
}
/////////////////////////////////////////////////////////////
int bchol(double a[],int n,double* det)//det指向存放返回行列式值(双精度实型)的变量地址
{ int i,j,k,u,l;//,v定义而不使用？
  double d;
  if ((a[0]+1.0==1.0)||(a[0]<0.0))
    { printf("fail\n"); return(-2);}
  a[0]=sqrt(a[0]);
  d=a[0];
  for (i=1; i<=n-1; i++)
    { u=i*n; a[u]=a[u]/a[0];}
  for (j=1; j<=n-1; j++)
    { l=j*n+j;
      for (k=0; k<=j-1; k++)
        { u=j*n+k; a[l]=a[l]-a[u]*a[u];}
      if ((a[l]+1.0==1.0)||(a[l]<0.0))
        { printf("fail\n"); return(-2);}
      a[l]=sqrt(a[l]);
      d=d*a[l];
      for (i=j+1; i<=n-1; i++)
        { u=i*n+j;
          for (k=0; k<=j-1; k++)
            a[u]=a[u]-a[i*n+k]*a[j*n+k];
          a[u]=a[u]/a[l];
        }
    }
  *det=d*d;
  for (i=0; i<=n-2; i++)
  for (j=i+1; j<=n-1; j++)
    a[i*n+j]=0.0;
  return(2);
}
/////////////////////////////////////////////////////////////
int blluu(double a[],int n,double l[],double u[])
{ int i,j,k,w,v,ll;
  for (k=0; k<=n-2; k++)
    { ll=k*n+k;
	if (fabs(a[ll])+1.0==1.0)
        { printf("fail\n"); return(0);}
      for (i=k+1; i<=n-1; i++)
	  { w=i*n+k; a[w]=a[w]/a[ll];}
      for (i=k+1; i<=n-1; i++)
        { w=i*n+k;
          for (j=k+1; j<=n-1; j++)
            { v=i*n+j;
              a[v]=a[v]-a[w]*a[k*n+j];
            }
        }
    }
  for (i=0; i<=n-1; i++)
    { for (j=0; j<i; j++)
        { w=i*n+j; l[w]=a[w]; u[w]=0.0;}
      w=i*n+i;
      l[w]=1.0; u[w]=a[w];
      for (j=i+1; j<=n-1; j++)
        { w=i*n+j; l[w]=0.0; u[w]=a[w];}
    }
  return(1);
}
/////////////////////////////////////////////////////////////
int bmaqr(double a[],int m,int n,double q[])
{ int i,j,k,l,nn,p,jj;
  double u,alpha,w,t;
  if (m<n)
    { printf("fail\n"); return(0);}
  for (i=0; i<=m-1; i++)
  for (j=0; j<=m-1; j++)
    { l=i*m+j; q[l]=0.0;