factortablefornr.cpp

该程序是一个含直流输电系统的潮流计算程序 程序运行需要包含Boost库文件,需要在VS2005环境以上(VC6.0不完全符合C++标准,可能有些变量作用域问题,无法通过调试) 潮流程序计算主要参考

#include "FactorTableForNR.h"

JacNR::JacNR()
{
	bInitalized=false;
}
JacNR::JacNR(int r)
{
	Initialize(r);
}

JacNR::~JacNR()
{
  if(bInitalized)
   {
	delete [] D;
	delete [] Ljs;
	delete [] Ujs;
   }
 if(bReadyForFactorize)
	{
		delete [] Lj;
		delete [] Uj;
		delete [] LA;
		delete [] UA;
	}
}

void JacNR::Initialize(int r)
{
 Jac.Initialize(r,r);
 D=new double[r];
 Ljs=new int[r+1];
 Ujs=new int[r+1];
 rows=r;
 cols=r;
 bInitalized=true;
 bFirstTime=true;
 bFactorized=false;
 bReadyForFactorize=false;
 bJacBuilded=false;
 ZeroFlagValue=2*THRESHOLD;
}
//这是在需要多次进行符号分解时采用的过程
void JacNR::RestoreForChangeElem()
{
	bJacBuilded=false;
	bFactorized=false;
	bReadyForFactorize=false;
	bFirstTime=true;
}

//设置元素
 void JacNR::SetJacElement(int r,int c,double value)
  {
	 int i,j;

	 if(bFirstTime)
	 {
	 //为了确保只进行一次符号分解，在所有可能出现非零元素的位置设置不为零的值
	   if(fabs(value)<THRESHOLD) 
		   value=ZeroFlagValue;	
	   Jac.SetElement(r,c,value);
	 }
	 else
	 {
		if(r==c)
			D[r]=value;
		else if(r>c)  //下三角
		{
			for(i=Ljs[c];i<Ljs[c+1];i++)
				if(Lj[i]==r)	LA[i]=value;
		}
		else
		{
			 for(j=Ujs[r];j<Ujs[r+1];j++)
		      if(Uj[j]==c) UA[j]=value;
		}
		
	 }
  }
//提取元素
//分两种情况
//1，如果是第一次设置元素，即还没有形成稀疏排列则直接
//从稀疏矩阵里读元素，否则从几个稀疏数组中获得元素
 double JacNR::GetJacElement(int r,int c)
 {
	int i,j;
	double value=0.0;
    if(!bReadyForFactorize)
	  value=Jac.GetElement(r,c);
	else
	{
	if(r==c)
		  value=D[r];
	  else if(r>c) //下三角
	  {
		  for(i=Ljs[c];i<Ljs[c+1];i++)
		      if(Lj[i]==r)
				  value=LA[i];
	  }
	  else
	  {
		 for(j=Ujs[r];j<Ujs[r+1];j++)
		      if(Uj[j]==c)
				  value=UA[j];
	  }	
	}
	return value;
 }



//在第一次借助链表稀疏矩阵设置好所有元素后
//统计各非零元素并按顺序存储到相应的变量中
 //
  bool JacNR::CountElement()
  {
	 int i,j;
	 int Lnum, Unum;
	 double elemv;

	 if(!bJacBuilded) return false;
	 if(!bFirstTime) return false;

	 SumLnum=0;SumUnum=0;  //从0开始计数
//下三角按实际的列存
	 for(j=0;j<cols;j++)
	 {
		 Lnum=0;
		 for(i=j+1;i<rows;i++)
		 {
			elemv=Jac.GetElement(i,j);
			if(fabs(elemv)>THRESHOLD) Lnum++;
		 }
		 Ljs[j]=SumLnum;
	  	 SumLnum+=Lnum;
	 }
//上三角按实际的行存
	 for(i=0;i<rows;i++)
	 {
		Unum=0;
		 for(j=i+1;j<cols;j++)
		 {
			elemv=Jac.GetElement(i,j);
			if(fabs(elemv)>THRESHOLD) Unum++;
		 }
		Ujs[i]=SumUnum;
		SumUnum+=Unum;
	 }
//最后一个
	   Ljs[rows]=SumLnum;		
	   Ujs[rows]=SumUnum;


	LA=new double[SumLnum];	 
	Lj=new int[SumLnum];	 
	UA=new double[SumUnum];
	Uj=new int[SumUnum];
//现在开始填充元素
	 SumLnum=0;SumUnum=0;
	 
	for(j=0;j<cols;j++)
	{
		for(i=j+1;i<rows;i++)		//对于行
		{
			elemv=Jac.GetElement(i,j);  //注意
			if(fabs(elemv)<THRESHOLD) continue;
			if(elemv==ZeroFlagValue)
				LA[SumLnum]=0;
			else
				LA[SumLnum]=elemv;

			Lj[SumLnum]=i;  //行号！
			SumLnum++;
			if(SumLnum>Ljs[j+1]) break;  //该列所有的非零元素都找到了
		}
	}
	for(i=0;i<rows;i++)
	{
		for(j=i+1;j<cols;j++)
		{
			elemv=Jac.GetElement(i,j);
			if(fabs(elemv)<THRESHOLD) continue;
			if(elemv==ZeroFlagValue)
				UA[SumUnum]=0;	
			else
				UA[SumUnum]=elemv;
			Uj[SumUnum]=j;
			SumUnum++;
			if(SumUnum>Ujs[i+1]) break;  //该行所有的非零元素都找到了	
		}

		 D[i]=Jac.GetElement(i,i);
	}

	   
		bReadyForFactorize=true;
		bFirstTime=false;
		return true;
  }

  bool JacNR::BuildFactorTable()
  {
	if(!bReadyForFactorize) return false;
	if(bFactorized) return true;
	int maxnum=(rows*rows-rows)/2;
	EstimateElemNum=int((rows*rows)*0.01);
	if(EstimateElemNum>maxnum) EstimateElemNum=maxnum;
	if(EstimateElemNum<6*rows) EstimateElemNum=6*rows;
	SolveEqs.SpareLUDecompose(rows,EstimateElemNum,Ljs,Ujs,Lj,Uj);


	bFactorized=true;
	return true;
  }

bool JacNR::ResolveEquation(double *pB)
  {
	bool bSuc;
	bSuc=SolveEqs.SpareLUfilldata(rows,pB,LA,UA,D,Ljs,Ujs,Lj,Uj);
	return bSuc;
  }