QR分解奇异值分解c语言编写

void QR(double a[N][N],double q[N][N],double r1[N][N],int n)      /*QR分解*/
{
   int i,j,k,r,m;
   double temp,sum,dr,cr,hr;
   double ur[N],pr[N],wr[N];
   double q1[N][N],emp[N][N];

   for(i=1;i<n+1;i++)//将a放入temp中
    for(j=1;j<n+1;j++)
    {
     emp[i][j]=a[i][j];
    };
   for(i=1;i<n+1;i++)//定义单位矩阵
    for(j=1;j<n+1;j++)
    {
     if(i==j)q[i][j]=1;
     else q[i][j]=0;
    };
   for(r=1;r<n;r++)
   {
    temp=0;
    for(k=r+1;k<n+1;k++)
    temp+=fabs(a[k][r]);
    if(temp>=ZERO)
     {
      sum=0;
      for(k=r;k<n+1;k++) 
      sum+=a[k][r]*a[k][r];
      dr=sqrt(sum);
      if(a[r][r]>ZERO)m=-1;