vector.cpp

矩阵运算

//////////////////////////////////////////////////////////////
//实矩阵相乘
//计算矩阵A（m*n）和B（n*k）的乘积，结果保存在C（m*k）中
//a-长度为m*n的数组
//b-长度为n*k的数组
//c-长度为m*k的数组，存放结果
void damul(double a[],double b[],int m,int n,int k,double c[]);
//////////////////////////////////////////////////////////////
//计算矩阵A（m*n）的转置矩阵AT（n*m）和B（m*k）的乘积，结果保存在C（n*k）中
//添加的函数，非原书程序
//a-长度为m*n的数组
//b-长度为m*k的数组
//c-长度为n*k的数组，存放结果
void ATdotB(double a[],double b[],int m,int n,int k,double c[]);
//////////////////////////////////////////////////////////////
//计算矩阵A（m*n）和B（k*n）的转置矩阵BT（n*k）的乘积，结果保存在C（m*k）中
//添加的函数，非原书程序
//a-长度为m*n的数组
//b-长度为k*n的数组
//c-长度为m*k的数组，存放结果
void AdotBT(double a[],double b[],int m,int n,int k,double c[]);
//////////////////////////////////////////////////////////////
//实矩阵求逆
//全选主元高斯-约当法
//a-长度为n*n的数组, n*n矩阵
//n 矩阵的维数
int dcinv(double a[],int n);
//////////////////////////////////////////////////////////////
//对称正定矩阵求逆
//a-长度为n*n的数组, n*n矩阵
//n 矩阵的维数
int desgj(double a[],int n);
//////////////////////////////////////////////////////////////
//托伯利兹(Toeplitz)矩阵求逆的特兰持(Trench)方法
//t-长度为n的数组,存放n阶T型矩阵中的上三角元素t0,t1,t2...tn-1
//tt-长度为n的数组,从tt[1]开始依次存放tt[1]...tt[n-1]
//n-矩阵的阶数
//b-长度为n*n的数组，返回时存放逆矩阵
int dftrn(double t[],double tt[],int n,double b[]);
//////////////////////////////////////////////////////////////
//求矩阵的行列式值
//全选主元高斯消去法
//a-长度为n*n的数组
//n-矩阵的阶数
double dhdet(double a[],int n);
//////////////////////////////////////////////////////////////
//对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值
//返回值小于0表示程序工作失败，还输出"fail";
//返回值大于0表示正常返回
//a-长度为n*n的数组，存放正定矩阵，
//  返回时下三角部分存放分解后的下三角矩阵L,其余元素为0
//n-正定矩阵的阶数
//det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值
int dicll(double a[],int n,double  *det);
//////////////////////////////////////////////////////////////
//矩阵的三角分解(LU)
//其中下三角阵L的主对角元素为1。
//a-长度为N*N的矩阵，返回时为L+U-I
//n-矩阵的阶数
//l-返回下三角矩阵
//u-返回上三角矩阵
int djlu(double a[],int n,double l[],double u[]);
//////////////////////////////////////////////////////////////
//实数矩阵的QR分解法
//用Householder变换对一般m*n阶的实数矩阵进行QR分解
//a-长度为m*n的一维数组,返回时其左上三角部分存放上三角矩阵R
//m-矩阵的行数
//n-矩阵的列数
//q-长度为m*m的矩阵，返回时存放正交矩阵Q
int dkqr(double a[],int m,int n,double q[]);
//////////////////////////////////////////////////////////////
//奇异值分解法求广义逆
//本函数返回值小于0表示在奇异值分解过程,
//中迭代值超过了60次还未满足精度要求.
//返回值大于0表示正常返回。
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//aa-长度为n*m的数组，返回式存放A的广义逆
//eps-精度要求
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav()
int dginv(double a[],int m,int n,double aa[],double eps,double u[],double v[],int ka);
//////////////////////////////////////////////////////////////
//实数矩阵的奇异值分解
//利用Householder变换及变形QR算法
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//eps-精度要求
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav(),ppp(),sss()
int dluav(double a[],int m,int n,double u[],double v[],double eps,int ka);
//////////////////////////////////////////////////////////////
#include "math.h"
#include "stdio.h"

//实矩阵相乘
//计算矩阵A（m*n）和B（n*k）的乘积，结果保存在C（m*k）中
//a-长度为m*n的数组
//b-长度为n*k的数组
//c-长度为m*k的数组，存放结果
void damul(double a[],double b[],int m,int n,int k,double c[])
{
int i,j,l,u;
    for (i=0; i<m; i++)
{
  for (j=0; j<k; j++)
  {
   u=i*k+j;
   c[u]=0.0;
   for (l=0; l<n; l++)
   {
    c[u]+=a[i*n+l]*b[l*k+j];
   }
  }
}
    return;
} //计算矩阵A（m*n）的转置矩阵AT（n*m）和B（m*k）的乘积，结果保存在C（n*k）中
//添加的函数，非原书程序
//a-长度为m*n的数组
//b-长度为m*k的数组
//c-长度为n*k的数组，存放结果
void ATdotB(double a[],double b[],int m,int n,int k,double c[])
{
int i,j,l,u;
    for (i=0; i<n; i++)
{
  for (j=0; j<k; j++)
  {
   u=i*k+j;
   c[u]=0.0;
   for (l=0; l<m; l++)
   {
    c[u]+=a[l*n+i]*b[l*k+j];
   }
  }
}
    return;
} //计算矩阵A（m*n）和B（k*n）的转置矩阵BT（n*k）的乘积，结果保存在C（m*k）中
//添加的函数，非原书程序
//a-长度为m*n的数组
//b-长度为k*n的数组
//c-长度为m*k的数组，存放结果
void AdotBT(double a[],double b[],int m,int n,int k,double c[])
{
int i,j,l,u;
    for (i=0; i<m; i++)
{
  for (j=0; j<k; j++)
  {
   u=i*k+j;
   c[u]=0.0;
   for (l=0; l<n; l++)
   {
    c[u]+=a[i*n+l]*b[j*n+l];
   }
  }
}
    return;
} //实矩阵求逆
//a-长度为n*n的数组, n*n矩阵
//n 矩阵的维数
int dcinv(double a[],int n)
{
int *is,*js,i,j,k,l,u,v;
double d,p;
is=new int[n];
is=malloc(n*sizeof(int));
    js=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);
   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;
    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);
}
            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);
} //托伯利兹(Toeplitz)矩阵求逆的特兰持(Trench)方法
//t-长度为n的数组,存放n阶T型矩阵中的上三角元素t0,t1,t2...tn-1
//tt-长度为n的数组,从tt[1]开始依次存放tt[1]...tt[n-1]
//n-矩阵的阶数
//b-长度为n*n的数组，返回时存放逆矩阵
int dftrn(double t[],double tt[],int n,double b[])
{
int i,j,k;
    double a,s,*c,*r,*p;
    c=malloc(n*sizeof(double));
    r=malloc(n*sizeof(double));
    p=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);
}


//求矩阵的行列式值
//全选主元高斯消去法
//a-长度为n*n的数组
//n-矩阵的阶数
double dhdet(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);
} //对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值
//返回值小于0表示程序工作失败，还输出"fail";
//返回值大于0表示正常返回
//a-长度为n*n的数组，存放正定矩阵，
//  返回时下三角部分存放分解后的下三角矩阵L,其余元素为0
//n-正定矩阵的阶数
//det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值
int dicll(double a[],int n,double  *det)
{ int i,j,k,u,v,l;
    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);
}


//求矩阵的行列式值
//全选主元高斯消去法
//a-长度为n*n的数组
//n-矩阵的阶数
double dhdet(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);
} //对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值
//返回值小于0表示程序工作失败，还输出"fail";
//返回值大于0表示正常返回
//a-长度为n*n的数组，存放正定矩阵，
//  返回时下三角部分存放分解后的下三角矩阵L,其余元素为0
//n-正定矩阵的阶数
//det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值
int dicll(double a[],int n,double  *det)
{ int i,j,k,u,v,l;
    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);
}


//矩阵的三角分解(LU)
//其中下三角阵L的主对角元素为1。
//a-长度为N*N的矩阵，返回时为L+U-I
//n-矩阵的阶数
//l-返回下三角矩阵
//u-返回上三角矩阵
int djlu(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);
} //实数矩阵的QR分解法
//用Householder变换对一般m*n阶的实数矩阵进行QR分解
//a-长度为m*n的一维数组,返回时其左上三角部分存放上三角矩阵R
//m-矩阵的行数
//n-矩阵的列数
//q-长度为m*m的矩阵，返回时存放正交矩阵Q
int dkqr(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;
   if (i==j)
   {
    q[l]=1.0;
   }
  }
}
nn=n;
if (m==n)
{
  nn=m-1;
}
for (k=0; k<=nn-1; k++)
{
  u=0.0;
  l=k*n+k;
  for (i=k; i<=m-1; i++)
  {
   w=fabs(a[i*n+k]);
   if (w>u)
   {
    u=w;
   }
  }
  alpha=0.0;
  for (i=k; i<=m-1; i++)
  {
   t=a[i*n+k]/u;
   alpha=alpha+t*t;
  }
  if (a[l]>0.0)
  {
   u=-u;
  }
  alpha=u*sqrt(alpha);
  if (fabs(alpha)+1.0==1.0)
  {
   printf("fail\n");
   return(0);
  }
  u=sqrt(2.0*alpha*(alpha-a[l]));
  if ((u+1.0)!=1.0)
  {
   a[l]=(a[l]-alpha)/u;
   for (i=k+1; i<=m-1; i++)
   {
    p=i*n+k;
    a[p]=a[p]/u;
   }
   for (j=0; j<=m-1; j++)
   {
    t=0.0;
    for (jj=k; jj<=m-1; jj++)
    {
     t=t+a[jj*n+k]*q[jj*m+j];
    }
    for (i=k; i<=m-1; i++)
    {
     p=i*m+j;
#   q[p]=q[p]-2.0*t*a[i*n+k];
    }
   }
   for (j=k+1; j<=n-1; j++)
   {
    t=0.0;
    for (jj=k; jj<=m-1; jj++)
    {
     t=t+a[jj*n+k]*a[jj*n+j];
    }
    for (i=k; i<=m-1; i++)
    {
     p=i*n+j;
     a[p]=a[p]-2.0*t*a[i*n+k];
    }
   }
   a[l]=alpha;
   for (i=k+1; i<=m-1; i++)
   {
    a[i*n+k]=0.0;
   }
  }
}
for (i=0; i<=m-2; i++)
{
  for (j=i+1; j<=m-1;j++)
  {
   p=i*m+j; l=j*m+i;
   t=q[p]; q[p]=q[l]; q[l]=t;
  }
}
return(1);
}


//奇异值分解法求广义逆
//本函数返回值小于0表示在奇异值分解过程,
//中迭代值超过了60次还未满足精度要求.
//返回值大于0表示正常返回。
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//aa-长度为n*m的数组，返回式存放A的广义逆
//eps-精度要求
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav()
int dginv(double a[],int m,int n,double aa[],double eps,double u[],double v[],int ka)
{ int i,j,k,l,t,p,q,f;
    i=dluav(a,m,n,u,v,eps,ka);
    if (i<0)
{
  return(-1);
}
    j=n;
    if (m<n)
{
  j=m;
}