lu文档.c

LU分解 用C的算法

#include<stdio.h>
#include<math.h>
#define N  3

double a[][N]={2,-2,3,1,1,1,1,3,-1};
void multmat(double a[N][N],double b[N][N],double c[N][N]) /*multiply two matrixex*/
{
  int i,j,k;
  for(i=0;i<N;i++)
	  for(j=0;j<N;j++)
	  {
		  c[i][j]=0;
		  for(k=0;k<N;k++)
			  c[i][j]+=a[i][k]*b[k][j];
	  }
}

int sign(double x)  /*seek the sign of a variable*/
{
  if(x<=0)
	  return(-1);
  else 
	  return(1);
}

double twomo(double x[N]) /*compute the "||x-y||*||x-y||"*/
{
  int i;
  double sum=0;
  for(i=0;i<N;i++)
	  sum+=pow(x[i],2);
  /*sum=sqrt(sum);*/
  return (sum);
}

void vectormat(double a[N],double c[N][N])   /*acquire (x-y)(x-y)T*/
{
int i,j;
  for(i=0;i<N;i++)
	  for(j=0;j<N;j++)
		  c[i][j]=a[i]*a[j];
	  
}


void main()
{
 int i,j,k;
 double H[N][N],temp[N][N],Q[N][N],h[N]={0,0,0};
 double sum;
 for(k=0;k<N-1;k++)
 {
   sum=0;
   for(i=k;i<N;i++)
	   sum+=a[i][k]*a[i][k];
   sum=sqrt(sum);
   for(i=k;i<N;i++)
   {	  
	   if(i==k)
		   h[i]=sum*sign(a[k][k])+a[k][k];
	   else
		   h[i]=a[i][k];
   }
	if(k>0)
		for(i=0;i<k;i++)
			h[i]=0;
 vectormat(h,H);
 for(i=0;i<N;i++)
	 for(j=0;j<N;j++)
		 if(i==j)
			 H[i][j]=1-2*H[i][j]/twomo(h);
		 else
             H[i][j]=-2*H[i][j]/twomo(h);
multmat(H,a,temp);
for(i=0;i<N;i++)
	 for(j=0;j<N;j++)
		 a[i][j]=temp[i][j];
if(k==0)
	 for(i=0;i<N;i++)
	    for(j=0;j<N;j++)
		    Q[i][j]=H[i][j];
else
{
    multmat(Q,H,temp);
	for(i=0;i<N;i++)
	 for(j=0;j<N;j++)
		 Q[i][j]=temp[i][j];
}

 }
 
printf("The answer Q  :\n");
for(i=0;i<N;i++)
{ 
	for(j=0;j<N;j++)
		 printf("%f\t",Q[i][j]);
	 printf("\n");
}

printf("\nThe answer R  :\n");
for(i=0;i<N;i++)
{ 
	for(j=0;j<N;j++)
		 printf("%f\t",a[i][j]);
	 printf("\n");
}
}