LU分解不分块的C程序

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<sys/time.h>
#include<math.h>
#define N                                       900
#define MAXRAND 3276700000.0
 
double sum;
long i=0,row=0,col=0,j=0,n=0;
double A[N][N];
double D[N][N];
double num;
double tstart,tstop,ttime;
double gflops = 0.0;


void InitA(double A[][N]);
void LU(double a[][N]);
void CheckResult (double a[][N],double b[][N] );
double dtime();




int main()
{


//初始化矩阵A
  InitA(A);


 tstart = dtime();


//对矩阵A进行LU分解col表示行坐标,row代表纵坐标


     LU(A);


     tstop = dtime();


     gflops = (double)(1.0e-9*(N-1)*N*((4*N+1)/6));


     if((ttime)>=0)
     {

       printf("Gflops=%lf,Secs = %lf, Gflops per sec = %lf\n",gflops,ttime,gflops/ttime);

     } 


//验证LU分解的结果是否正确
 CheckResult ( A,D);
   
  return 0;
}

double dtime()
{
   double tseconds = 0.0;
   struct timeval mytime;
   gettimeofday(&mytime,(struct timezone*)0);
   tseconds = (double)(mytime.tv_sec + mytime.tv_usec*1.0e-6);
   return tseconds;


}
void InitA(double A[][N])
{


  srand((unsigned int)time(NULL));
    for(i=0;i<N;i++)
     {


          for(j=0;j<N;j++)
          {
             num = ( double )rand()/RAND_MAX;
             A[i][j] = num;
             D[i][j] = num;
             if(i==j)
              {
                 A[i][j]=1.0;
                 D[i][j]=1.0;
              }
   
          }
    
     }
}

//对矩阵A进行LU分解col表示列坐标,row代表行坐标
void LU(double a[][N])
{


   for(i=0;i<N;i++)
   {
 
      for(row=i+1;row<N;row++)
       {
          a[row][i]=a[row][i]/a[i][i];
          for(col=i+1;col<N;col++)
          {
 
             a[row][col]=a[row][col]-a[row][i]*a[i][col];
          }
      
       }
    }

}


//验证LU分解的结果是否正确
void CheckResult (  double a[][N],double b[][N] )
{


  double max = 0.0000001;
  double B[N][N] = {0},C[N][N]={0};
   for(i=0;i<N;i++)
    {


      for(j=0;j<N;j++)
       {


           if(i<j)
            {
              C[i][j]=a[i][j];


            }else if(i>j)
            {


              B[i][j]=a[i][j];


            }else if(i==j)
           {


             C[i][j]=a[i][j];
             B[i][j]=1;
           }
       }
   
    }

  for(n=0;n<N;n++)
  {
       for(row=0; row<N; row++)
       {
            sum=0.0;
            for(col=0; col<N; col++)
            {
              sum += B[ n ][ col ]*C[ col ][ row ];
            }
              //printf("%lf %lf\n",b[n][row],sum);
              b[n][row] = (double)fabs(b[n][row]-sum);

       }
   }

//求误差的最大值

   for(n=0;n<N;n++)
    {


      for(i=0;i<N;i++)


         if(b[n][i]>max)


               max=b[n][i];


    }  


printf("最后差的最大值是%lf\n",max);


    
}