lu.cpp

用于集成电路辅助设计的一写算法

#include "lu.h"
#include "stdio.h"
#include "stdlib.h"
#include "iostream.h"


void LU_decompose(float matrix_a[][500],float*matrix_L,float*matrix_U,float*vector_b,int order)
{
	int i,j,k;
	float *ptr=new float[order*order];
	for(i=0;i<order;i++) 
		for(j=0;j<order;j++)
			ptr[i*order+j]=matrix_a[i][j];
	for(i=0;i<order-1;i++)
	{
		gauss(ptr,vector_b,i,order);
		for(j=i+1;j<order;j++) 
		{
			ptr[i*order+j]=ptr[i*order+j]/ptr[i*order+i];

			for(k=i+1;k<order;k++)			
				ptr[k*order+j]=ptr[k*order+j]-ptr[k*order+i]*ptr[i*order+j];
			
		}
	}
	for(i=0;i<order;i++)
	{
		for(j=0;j<order;j++)
		{
			if(j>i) 
			{ 
				matrix_U[i*order+j]=ptr[i*order+j];
				matrix_L[i*order+j]=0;
			}
			else if(i==j) 
			{
				matrix_L[i*order+j]=ptr[i*order+j];
				matrix_U[i*order+j]=1;
			}
			else
			{
				matrix_L[i*order+j]=ptr[i*order+j]; 
				matrix_U[i*order+j]=0;
			}
		}
	}
	delete[] ptr;

}


void solve_L(float*matrix,float*vector,float*ans,int order)
{
	for(int i=0;i<order;i++)
	{
		for(int j=0;j<i;j++) vector[i]=vector[i]-matrix[i*order+j]*ans[j];
		ans[i]=vector[i]/matrix[i*order+i];
	}
}


void solve_U(float*matrix,float*vector,float*ans,int order)
{
	int i,j;
	for(i=order-1;i>=0;i--)
	{
		for(j=order-1;j>i;j--) vector[i]=vector[i]-matrix[i*order+j]*ans[j];
		ans[i]=vector[i];
	}
}


void gauss(float*matrix,float*vector,int n,int order)
{
	float max,tempf;
	int temp,i;
	temp=n;
	max=absolute(matrix[n*order+n]);
	for(i=n;i<order;i++)
		if(max<absolute(matrix[i*order+n]))
		{
			max=absolute(matrix[i*order+n]);
			temp=i;
		}

		if(temp!=n)
		{
			tempf=vector[n];
			vector[n]=vector[temp];
			vector[temp]=tempf;
			for(i=0;i<order;i++)
			{
				tempf=matrix[n*order+i];
				matrix[n*order+i]=matrix[temp*order+i];
				matrix[temp*order+i]=tempf;
			}
		}
}

float absolute(float num)
{
	if(num>=0) return num;
	else return(-num);
}

void displaymatrix(float* ptr,int order)
{
	float num;
	int i;
	for(i=0;i<order;i++)
	{
		for(int j=0;j<order;j++)
		{	
			num=ptr[i*order+j];
			num=absolute(num)-int(absolute(num));
			if(num==0) printf("%.0f   ",ptr[i*order+j]);
			else printf("%.4f   ",ptr[i*order+j]);						
		}
		printf("\n\n");
	}
}