improvecholesky.cpp

用改进的乔里斯基方法解对称正定方程组

//=======================================================//
//函数名称：improvecholesky()                            //
//函数目的：用改进的cholesky算法解大型对称正定方程组     //
//参数说明：A：存放方程组系数矩阵                        //
//          B：存放方程组的右端项                        //
//          n：方程组的阶数                              //
//          w: 方程组的带宽                              //
//          m：右端项的组数                              //
//=======================================================// 
#include<math.h>
#include "iostream.h"
int improvecholesky(double *A,double *B,int n,int w,int m)
{
	int i,j,k,u,v;
	double temp;
	double *p=new double[w];
	///////////////////////////////////////////////////////
	///假如*A等于零，说明系数矩阵非正定，返回-1。
	if(fabs(A[0]+1.0)==1.0)return -1;
	///////////////////////////////////////////////////////
	//对系数矩阵A进行LDL'分解，求L和D。
	for(i=1;i<n;i++)
	{	
		///////////////////////////////////////////////////
		//if(i<w)元素的定位都从零开始
		if(i<w)
		{
			//求取U
			for(j=0;j<i;j++)
			{
				temp=0;
				for(k=0;k<j;k++)
				{
					u=(j-k)*n+k;
					temp=temp+A[u]*p[k];
				}
				u=(i-j)*n+j;
				p[j]=A[u]-temp;		
			}
			
			//求取L
			for(j=0;j<i;j++)
			{
				u=(i-j)*n+j;
				A[u]=p[j]/A[j];
				//cout<<A[(i-j)*n+j]<<endl;
			}
			//求取D
			temp=0.0;
			for(k=0;k<i;k++)
			{
				u=(i-k)*n+k;
				temp=temp+A[u]*p[k];
			}
			A[i]=A[i]-temp;
			//当D为零时返回-1
			if(fabs(A[i]+1.0)==1.0)return -1;
		}
		///////////////////////////////////////////////////
		//if(i>=w)元素的定位都从position开始
		else
		{
			int position;
			position=i-w+1;
			for(j=position;j<i;j++)
			{
				temp=0;
				for(k=position;k<j;k++)
				{
					u=(j-k)*n+k;
					v=k-position;
					temp=temp+A[u]*p[v];
				}
				u=(i-j)*n+j;
				v=j-position;
				p[v]=A[u]-temp;
			}
			//求取L
			for(j=position;j<i;j++)
			{
				u=(i-j)*n+j;
				v=j-position;
				A[u]=p[v]/A[j];
			}
			//求取D
			temp=0.0;
			for(k=position;k<i;k++)
			{
				u=(i-k)*n+k;
				v=k-position;
				temp=temp+A[u]*p[v];
			}
			A[i]=A[i]-temp;
			//当D为零时返回-1
			if(fabs(A[i]+1.0)==1.0)return -1;
		}
	} 
	delete p;
	///////////////////////////////////////////////////////
	//回代过程
	for(j=0;j<m;j++)
	{
		//求取Y
		for(i=0;i<n;i++)
		{
			temp=0;
			if(i>=w)k=i-w+1;
			else    k=0;
			for(;k<i;k++)
			{
				u=(i-k)*n+k;
				v=k*m+j;
				temp=temp+A[u]*B[v];
			}
			v=i*m+j;
			B[v]=B[v]-temp;
		}
		//求取x
		for(i=n-1;i>=0;i--)
		{
			int s;
			temp=0;
			if(i<=n-w-1)s=n-(n-w-i);
			else        s=n;
			for(k=i+1;k<s;k++)
			{
				u=(k-i)*n+i;
				v=k*m+j;
				temp=temp+A[u]*B[v];				
			}
			v=i*m+j;
			B[v]=B[v]/A[i]-temp;
		}
	}
	return 1;
}

#include "iostream.h"
#include "iomanip.h"
void main()

  { int i;
	double a[12]={ 5,6,6,5,-4,-4,-4,0,1,1,0,0};
     double c[8]={2,2,-1,-1,-1,-1,2,2};

    if (improvecholesky(a,c,4,3,2)>0)
      for (i=0; i<4; i++)
		 cout<<c[i*2]<<ends<<c[i*2+1]<<endl;
		cout<<endl;
	
  }