单纯lu分解.cpp

下一部是计算求解方程组

#include <stdio.h>
#include <ctype.h>
#include <stdlib.h>
#include <iostream.h>
#include <malloc.h>
#define list_init_size 100
#define listincrement 10
#include <process.h>
#define overflow -2
#define ok 1
#define error 0
#define elemtype  float

typedef struct{
elemtype *elem;
int length;
int listsize;
}sqlist;

initlist_sq (sqlist &l)
 {
 l.elem=(elemtype *) malloc(list_init_size*sizeof(elemtype));
 if (!l.elem) exit(overflow);
 l.length=0;
 l.listsize=list_init_size;
 return 1;
 }

listinsert_sq(sqlist &l,int i,elemtype e)
{elemtype *newbase;
elemtype *p,*q;
 if (i<1||i>l.length+1) return error;
 if (l.length>=l.listsize)
    {
     newbase=(elemtype *) realloc(l.elem,(l.listsize+listincrement)*sizeof(elemtype));
     if (!newbase) exit (overflow);
     l.elem=newbase;
     l.listsize+=listincrement;
    }
 q=&(l.elem[i-1]);
 for (p=&(l.elem[l.length-1]);p>=q;--p)
     *(p+1)=*p;
 *q=e;
 ++l.length;
 return ok;
}




void main()
{ int n,t,i,j,k,r,temp=0;
elemtype *pa,*pl,*pu;
cout<<"Input n:"<endl; //A=n*n//
cin>>n;
sqlist a,l,u;
initlist_sq(a);
cout<<"Input A:";    //Input A//
for (i=0;i<n*n;i++)
    {
     cin>>t;
    listinsert_sq(a,i+1,t);
    }
cout<<endl;

pa=a.elem;
pl=l.elem=(elemtype *) malloc(a.listsize*sizeof(elemtype));
pu=u.elem=(elemtype *) malloc(a.listsize*sizeof(elemtype));
l.length=u.length=n*n;

for (i=0;i<n*n;i++)   pu[i]=pl[i]=0;

pu[0]=pa[0];  //Comput the value of L and U//
pl[0]=1;
 for (i=1;i<n;i++)
    {
     pl[i*n+i]=1;
     pu[i]=pa[i];
     pl[i*n]=pa[i*n]/pu[0];
    }

 for (k=1;k<n;k++)
     {
      for (j=k;j<n;j++)
	  {
	   for (r=0;r<=k-1;r++)
	       {
		 temp=temp+pl[k*n+r]*pu[r*n+j];
	       }
	   pu[k*n+j]=pa[k*n+j]-temp;
	   temp=0;
	  }
      for (i=k+1;i<n;i++)
	  {
	   for (r=0;r<=k-1;r++) temp=pl[i*n+r]*pu[r*n+k]+temp;
	   pl[i*n+k]=(pa[i*n+k]-temp)/pu[k*n+k];
	   temp=0;
	  }
     }

  cout<<"L"<<endl;  //Print L//
  for (i=0;i<n;i++)
    {
     for (j=0;j<n;j++)
     printf("%5.3f  ",pl[i*n+j]);
     cout <<"\n";
    }
 cout <<"U"<<"\n";     //Print U//
 for (i=0;i<n;i++)
    {
     for (j=0;j<n;j++)
     printf("%5.3f  ",pu[i*n+j]);
     cout <<"\n";
    }
 cout <<"\n";
}