[C语言]矩阵的LU分解

void main()
{
 float a[4][4]={{4,2,1,5},{8,7,2,10},{4,8,3,6},{12,6,11,20}};
 float L[4][4],U[4][4];
 int n=4;
 int k,i,j;
 float s,t;

 for(j=0;j<n;j++)
  a[0][j]=a[0][j];

 for(i=1;i<n;i++)
  a[i][0]=a[i][0]/a[0][0];

 for(k=1;k<n;k++)
 {
  for(j=k;j<n;j++)
  {
   s=0;
   for (i=0;i<k;i++)
    s=s+a[k][i]*a[i][j];
   a[k][j]=a[k][j]-s;
  }
  for(i=k+1;i<n;i++)
  {
   t=0;
   for(j=0;j<k;j++)
    t=t+a[i][j]*a[j][k];
   a[i][k]=(a[i][k]-t)/a[k][k];
  }
 }

 for(i=0;i<n;i++)
  for(j=0;j<n;j++)
  { if(i>j)
     { L[i][j]=a[i][j]; U[i][j]=0;}
    else
     { U[i][j]=a[i][j];
       if(i==j) L[i][j]=1;
       else L[i][j]=0;
     }
  }
 printf("/nL[4][4]=");
 for(i=0;i<n;i++)
  { printf("/n");
    for(j=0;j<n;j++)
     printf(" %0.0f",L[i][j]);
  }

 printf("/nU[4][4]=");
 for(i=0;i<n;i++)
 { printf("/n");
   for(j=0;j<n;j++)
    printf(" %0.0f",U[i][j]);
 }
}