factortableforpq.cpp

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

#include "FactorTableForPQ.h"


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

JacPQ::~JacPQ()
{
  if(bInitalized)
  {
	delete [] D;
	delete [] Ljs;
  }
 if(bReadyForFactorize)
 {
	delete [] Lj;
	delete [] LA;
 }
}

void JacPQ::Initialize(int r)
{
 B.Initialize(r,r);
 D=new double[r];
 Ljs=new int[r+1];
 rows=r;
 cols=r;
 bInitalized=true;
 bFirstTime=true;
 bFactorized=false;
 bReadyForFactorize=false;
 bBBuilded=false;
}

//设置元素
 void JacPQ::SetBElement(int r,int c,double value)
  {
	 int i,j;
	 if(bFirstTime)
	 {
	 //为了确保只进行一次符号分解，在所有可能出现非零元素的位置设置不为零的值
		  B.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=Ljs[c];j<Ljs[c+1];j++)
		      if(Lj[j]==r) LA[j]=value;
		}
		
	 }
	
  }
//提取元素
//分两种情况
//1，如果是第一次设置元素，即还没有形成稀疏排列则直接
//从稀疏矩阵里读元素，否则从几个稀疏数组中获得元素
 double JacPQ::GetBElement(int r,int c)
 {
	int i;
	double value=0.0;
    if(!bReadyForFactorize)
	  value=B.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(i=Ljs[r];i<Ljs[r+1];i++)
		      if(Lj[i]==c)
				  value=LA[i];
	  }	
	}
	return value;
 }


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

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

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

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

			Lj[SumLnum]=i;  //行号！
			SumLnum++;
			if(SumLnum>Ljs[j+1]) break;  //该列所有的非零元素都找到了
		}
		 D[j]=B.GetElement(j,j);	
	}

	   
		bReadyForFactorize=true;
		bFirstTime=false;

		return true;
  }

  bool JacPQ::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<5*rows) EstimateElemNum=5*rows;
	SolveEqs.SetDecompTyepFlag(false);
	SolveEqs.SpareLDLDecompose(rows,EstimateElemNum,Ljs,Lj);
	bFactorized=true;
	return true;
  }

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