Matlab中QR分解的C++代码

void QR(double *a,int a_row,int a_colunm,double* q,double *r1)
{
 //QR分解

 int i,j,k,r,m;
 double temp,sum,dr,cr,hr;
 double *ur=new double[a_row*a_row];
 double *pr=new double[a_row*a_row];
 double *wr=new double[a_row*a_row];

 double *q1=new double[a_row*a_row];
 double *emp=new double[a_row*a_row];

 for(i=0;i<a_row;i++)//将a放入temp中

  for(j=0;j<a_colunm;j++)
  {
   emp[i*a_colunm+j]=a[i*a_colunm+j];
  };
 for(i=0;i<a_row;i++)//定义单位矩阵

  for(j=0;j<a_row;j++)
  {
   if(i==j)q[i*a_row+j]=1;

   else q[i*a_row+j]=0;
  };

 for(r=0;r<a_colunm;r++)
 {
  temp=0;

  for(k=r;k<a_row;k++)

   temp+=fabs(a[k*a_colunm+r]);

  if(temp>=0.0)
  {
   sum=0;

   for(k=r;k<a_row;k++)

    sum+=a[k*a_colunm+r]*a[k*a_colunm+r];
   dr=sqrt(sum);
   if(a[r*a_colunm+r]>0.0)m=-1;

   else m=1;
   cr=m*dr;
   hr=cr*(cr-a[r*a_colunm+r]);

   for(i=0;i<a_row;i++)//定义ur
   {
    if(i<r)ur[i]=0;

    if(i==r)ur[i]=a[r*a_colunm+r]-cr;

    if(i>r)ur[i]=a[i*a_colunm+r];
   };

   for(i=0;i<a_row;i++)//定义wr
   {
    sum=0;

    for(j=0;j<a_row;j++)

     sum+=q[i*a_row+j]*ur[j];
    wr[i]=sum;
   };

   for(i=0;i<a_row;i++)//定义qr
    for(j=0;j<a_row;j++)
    {
     q1[i*a_row+j]=q[i*a_row+j]-wr[i]*ur[j]/hr;
    };
   for(i=0;i<a_row;i++)//定义qr+1
    for(j=0;j<a_row;j++)
    {
     q[i*a_row+j]=q1[i*a_row+j];
    };
   for(i=0;i<a_row;i++)//定义pr
   {
    sum=0;
    for(j=0;j<a_row;j++)
     sum+=a[j*a_colunm+i]*ur[j];
    pr[i]=sum/hr;
   };

   for(i=0;i<a_row;i++)
    for(j=0;j<a_colunm;j++)
    {
     a[i*a_colunm+j]=a[i*a_colunm+j]-ur[i]*pr[j];
    };
  };
 };
 for(i=0;i<a_row;i++)
  for(j=0;j<a_colunm;j++)
  {
   if(fabs(a[i*a_colunm+j])<0.0)a[i*a_colunm+j]=0;
  };

 for(i=0;i<a_row;i++)
  for(j=0;j<a_colunm;j++)
  {
   r1[i*a_colunm+j]=a[i*a_colunm+j];
  };

 for(i=0;i<a_row;i++)//将a取出
  for(j=0;j<a_colunm;j++)
  {
   a[i*a_colunm+j]=emp[i*a_colunm+j];
  }
}