cresolveeqonludecomp.cpp

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

#include "CResolveEqOnLuDecomp.h"

void SolveEqOnLU::SpareLDLDecompose(int wide,int SpareElementNumber,int* Jjs,int* Jj)
{
    long n,i,j,k,kl;

    int *w=new int[wide];
	int *Z=new int[wide];
	int *jo=new int[wide];
    int *r=new int[SpareElementNumber];
	int *js=new int[wide];
	int *H=new int[SpareElementNumber];

	m_LSpare.Jj=new int[SpareElementNumber];
	m_LSpare.Jjs=new int[wide+1];
	m_LSpare.Ji=new int[SpareElementNumber];
	m_LSpare.Jis=new int[wide+1];
	m_LSpare.jtoi=new int[SpareElementNumber];

	for(i=0;i<wide+1;i++)
	{
        m_LSpare.Jis[i]=0;
		m_LSpare.Jjs[i]=0;
	}
    for(i=0;i<SpareElementNumber;i++)
	{
        m_LSpare.Ji[i]=0;
		m_LSpare.Jj[i]=0;
		m_LSpare.jtoi[i]=0;
	}
    for(i=0;i<wide;i++)
	{
        w[i]=0;
		Z[i]=0;
        jo[i]=0;
		js[i]=0;
	}
    for(i=0;i<SpareElementNumber;i++)
	{
        r[i]=0;
        H[i]=0;
	}
    //符号分解
    m_LSpare.Jis[0]=0;
    m_LSpare.Jis[1]=0;
    m_LSpare.Jjs[0]=0;
    m_LSpare.Jjs[1]=Jjs[1];
    n=m_LSpare.Jjs[1];
    for(j=0;j<wide;j++)
        m_LSpare.Jj[j]=Jj[j];

    for(k=0;k<=Jjs[1]-1;k++) //!!!
	{
        H[k]=0;
		r[k]=0;
        jo[Jj[k]]=k;
        js[Jj[k]]++;
    }

    for(j=1;j<wide;j++)
	{
        for(i=Jjs[j];i<=Jjs[j+1]-1;i++)
            w[Jj[i]]++;
        i=jo[j];
        m_LSpare.Jis[j+1]=m_LSpare.Jis[j]+js[j];
        for(kl=m_LSpare.Jis[j+1]-1;kl>=m_LSpare.Jis[j];kl--)
		{
            m_LSpare.Ji[kl]=H[i];
            i=r[i];
        }
        for(kl=m_LSpare.Jis[j];kl<=m_LSpare.Jis[j+1]-1;kl++)
		{
            k=m_LSpare.Ji[kl];
            for(i=m_LSpare.Jjs[k]+Z[k];i<=m_LSpare.Jjs[k+1]-1;i++)
                w[m_LSpare.Jj[i]]++;
            Z[k]++;
        }
        for(kl=j+1;kl<wide;kl++)
		{
            if(w[kl]!=0)
			{
                m_LSpare.Jj[n]=kl;
                H[n]=j; 
				r[n]=jo[kl];
				jo[kl]=n; 
                js[kl]++;
                n++;
                w[kl]=0;
            }
        }
        m_LSpare.Jjs[j+1]=n;
    }
//这里可以确定非零元素的个数了
	SpareElementNum=n+1;

    SolveEqOnLU::Formjtoi(wide,m_LSpare.Jjs,m_LSpare.Jj,m_LSpare.Jis,m_LSpare.jtoi);
  
    bDecomposed=true;  //标示已经分解

	delete []w;
	delete []Z;
	delete []jo;
    delete []r;
	delete []js;
	delete []H;

}

bool SolveEqOnLU::SpareLDLfilldata(int wide,double* inoutvector,int* Jjs,int* Jj,double* Ja,double* JD)
{
	int i,j,k,kl,km;
	double Sumx,ReadValue,ReadSum; 

	double  *midleVariable=new double[wide];
	double *unknow=new double[wide];
	double *SumL=new double[wide];
	double *diagonal=new double[wide];
	int *Z=new int[wide];
	double *Lje=new double[SpareElementNum];

	for(i=0;i<wide;i++)
	{
		midleVariable[i]=0;
		unknow[i]=0;
		Z[i]=0;
		diagonal[i]=0;
	}
	for(i=0;i<SpareElementNum;i++)
		Lje[i]=0;

    //求 L，D 的值
    
    for(j=0;j<wide;j++)
        SumL[j]=0;
    for(j=0;j<wide;j++)
	{
        diagonal[j]=JD[j];
        for(kl=Jjs[j];kl<=Jjs[j+1]-1;kl++)
		{
            k=Jj[kl];
            SumL[k]=Ja[kl];		//第j列原始的元素值
        }
        for(kl=m_LSpare.Jis[j];kl<=m_LSpare.Jis[j+1]-1;kl++)
		{
			k=m_LSpare.Ji[kl];
            ReadValue=Lje[m_LSpare.jtoi[kl]];
            ReadSum=diagonal[k]*ReadValue;   //Ljk*Dk
            diagonal[j]-=ReadValue*ReadSum;	 //Dj=Ljk*Dk*Ljk
            for(km=m_LSpare.Jjs[k]+Z[k];km<=m_LSpare.Jjs[k+1]-1;km++)
			{
				i=m_LSpare.Jj[km];		//获得非零元素行号i
                SumL[i]-=Lje[km]*ReadSum;	//计算k列对Lij的贡献
            }
            Z[k]++;		//动态识别第k列已经用过的元素
        }
        for(kl=m_LSpare.Jjs[j];kl<=m_LSpare.Jjs[j+1]-1;kl++)	//对j列所有非零元素
		{
            i=m_LSpare.Jj[kl];
            Lje[kl]=SumL[i]/diagonal[j];
            SumL[i]=0;
        }
    }
    
    //前代
    
    for(i=0;i<wide;i++)
        midleVariable[i]=inoutvector[i]; //赋右端向量的值
    for(i=0;i<wide;i++)//按行前代
	{
        for(k=m_LSpare.Jjs[i];k<=m_LSpare.Jjs[i+1]-1;k++)
		{
            j=m_LSpare.Jj[k];
            midleVariable[j]-=Lje[k]*midleVariable[i];	
        }
    }
    
    //回代
    
    for(i=wide-1;i>=0;i--)	//按列回代
    {
		Sumx=0;
        for(k=m_LSpare.Jjs[i];k<=m_LSpare.Jjs[i+1]-1;k++)
		{
            j=m_LSpare.Jj[k];
            Sumx-=Lje[k]*unknow[j];
        }
        unknow[i]=midleVariable[i]/diagonal[i]+Sumx;
		inoutvector[i]=unknow[i];
	}    

	delete []midleVariable;
	delete []SumL;
	delete []diagonal;
	delete []Z;
	delete []Lje;
/*
	delete []m_LSpare.Jj;
	delete []m_LSpare.Jjs;
	delete []m_LSpare.Ji;
	delete []m_LSpare.Jis;
	delete []m_LSpare.jtoi;
*/
	return true;

}


bool SolveEqOnLU::Formjtoi(int wide,int* jtoiJjs, int* jtoiJj,int* jtoiJis,int * jtoiArray)
{
	long k;
	int j,kl;

    int *Z;
	Z=new int[wide];


    for(k=0;k<wide;k++)
        Z[k]=0;

   
    jtoiArray[0]=1;

    for(k=0;k<wide;k++)
	{
        for(kl=jtoiJjs[k];kl<=jtoiJjs[k+1]-1;kl++)
		{
            j=jtoiJj[kl]; //表示按列存储的kl个元素的行号
		//jtoiJis[j]表示第j行的元素开始的列号，Z[j]表示第j行的元素个数
		//所以表示
            jtoiArray[jtoiJis[j]+Z[j]]=kl; 
            Z[j]++; 
        }
    }

	delete []Z;
	return true;
}


//LU分解
void SolveEqOnLU::SpareLUDecompose(int wide,int SpareElementNumber,int* Jlowjs, int* Jupjs,int* Jlowj,int* Jupj)
{
	int il,iu,j,l,k,i,nl,nu;	

    int* rl;
	rl=new int[SpareElementNumber];

	int *ru;
	ru=new int[SpareElementNumber];
	
	int* hl;
	hl=new int[SpareElementNumber];
	
	int*hu;
	hu=new int[SpareElementNumber];

    int *zl;
	zl=new int[wide+1];

	int *zu;
	zu=new int[wide+1];
	
	int *jlo;
	jlo=new int[wide];
	
	int *juo;
	juo=new int[wide];

    int *Jls;
	Jls=new int[wide];
	
	int *Jus;
	Jus=new int[wide];
	
	int *wl;
	wl=new int[wide];
	
	int *wu;
	wu=new int[wide];
	
	m_LSpare.Jj=new int[SpareElementNumber];
	m_LSpare.Jjs=new int[wide+1];
	m_LSpare.Ji=new int[SpareElementNumber];
	m_LSpare.Jis=new int[wide+1];
	m_LSpare.jtoi=new int[SpareElementNumber];

	m_USpare.Jj=new int[SpareElementNumber];
	m_USpare.Jjs=new int[wide+1];
	m_USpare.Ji=new int[SpareElementNumber];
	m_USpare.Jis=new int[wide+1];
	m_USpare.jtoi=new int[SpareElementNumber];


	for(i=0;i<SpareElementNumber;i++)
	{
		m_LSpare.Ji[i]=0;
		m_LSpare.Jj[i]=0;
		m_LSpare.jtoi[i]=0;
		m_USpare.Ji[i]=0;
		m_USpare.Jj[i]=0;
		m_USpare.jtoi[i]=0;
		rl[i]=0,ru[i]=0;
        hl[i]=0,hu[i]=0;
	}
	for(i=0;i<wide;i++)
	{
		jlo[i]=0,juo[i]=0;
        Jls[i]=0,Jus[i]=0;
        wl[i]=0,wu[i]=0;
	}

	for(i=0;i<=wide;i++)
	{
		m_LSpare.Jis[i]=0;
		m_LSpare.Jjs[i]=0;
		m_USpare.Jis[i]=0;
		m_USpare.Jjs[i]=0;
		zl[i]=0,zu[i]=0;
	}

    m_LSpare.Jis[0]=0;
	m_USpare.Jis[0]=0;
    m_LSpare.Jis[1]=0;
	m_USpare.Jis[1]=0;
    m_LSpare.Jjs[0]=0;
	m_USpare.Jjs[0]=0;
    m_LSpare.Jjs[1]=Jlowjs[1];
	m_USpare.Jjs[1]=Jupjs[1];
    nl=m_LSpare.Jjs[1];
	nu=m_USpare.Jjs[1];

    for(j=0;j<wide;j++)
	{
        wl[j]=0;
		wu[j]=0;
        zl[j]=0;
		zu[j]=0;
        jlo[j]=0;
		juo[j]=0;
        Jls[j]=0;
		Jus[j]=0;
        m_LSpare.Jj[j]=Jlowj[j];
		m_USpare.Jj[j]=Jupj[j];
    }
    for(k=0;k<=Jlowjs[1]-1;k++)  //第一列非零元素
	{
        hl[k]=0;		//h用于装列下标
        jlo[Jlowj[k]]=k;	//第一列非零元素行所对应的非零序号
        Jls[Jlowj[k]]++;	//第一列非零元素所对应的行作标记
    }
    for(k=0;k<=Jupjs[1]-1;k++)
	{
        hu[k]=0;
        juo[Jupj[k]]=k;
        Jus[Jupj[k]]++;
	}
    for(j=1;j<wide;j++)	//从第二列到wide列
	{
        for(i=Jlowjs[j];i<=Jlowjs[j+1]-1;i++)  //第j列
			wl[Jlowj[i]]++;				//标记第j列原有非零元素的行号
        for(i=Jupjs[j];i<=Jupjs[j+1]-1;i++)
			wu[Jupj[i]]++;
        il=jlo[j];		//第j行非零元素的起点
        iu=juo[j];
        m_LSpare.Jis[j+1]=m_LSpare.Jis[j]+Jls[j];  //统计非零元素的起点，定义第j+1行头元素的位置
        m_USpare.Jis[j+1]=m_USpare.Jis[j]+Jus[j];
        for(l=m_LSpare.Jis[j+1]-1;l>=m_LSpare.Jis[j];l--) //对于第j行！！！
		{ //逆推,这样第j行的非零元素对应列的列号是按从大到小排列的。
            m_LSpare.Ji[l]=hl[il];  //第j行非零元素的所对应的列     
            il=rl[il];		//与对il个非零元素同列的下一个元素的位置
		}
        for(l=m_USpare.Jis[j+1]-1;l>=m_USpare.Jis[j];l--)
		{
            m_USpare.Ji[l]=hu[iu];
            iu=ru[iu];
        }
        for(l=m_LSpare.Jis[j];l<=m_LSpare.Jis[j+1]-1;l++)//对于第j行
		{
            k=m_LSpare.Ji[l];   //取非零元素的列号
            for(i=m_USpare.Jjs[k]+zu[k];i<=m_USpare.Jjs[k+1]-1;i++) //对该列中行号大于j的非零元素
                wu[m_USpare.Jj[i]]++;
            zu[k]++;
        }
        for(l=m_USpare.Jis[j];l<=m_USpare.Jis[j+1]-1;l++)  //对于第j行
		{
            k=m_USpare.Ji[l];			//取非零元素的列号
            for(i=m_LSpare.Jjs[k]+zl[k];i<=m_LSpare.Jjs[k+1]-1;i++)
                wl[m_LSpare.Jj[i]]++;
	//动态地识别第j行的元素对应的列中被使用的元素
	//当前分解j列，用U的第j行的非零元素列，与L的k列非零元素相对应
	//zl用于识别k列中行号大于j的非零元素。
            zl[k]++;	
        }
        for(l=j+1;l<wide;l++)	//对j以下的行标示出注入元素
		{
            if(wl[l]!=0)
			{
                m_LSpare.Jj[nl]=l;
                hl[nl]=j;		//记录列号
                rl[nl]=jlo[l];	//下一个同列元素的位置
                jlo[l]=nl;	//l行的元素对应于hl的位置
                Jls[l]++;	//第l行非零元素的个数
                wl[l]=0;	//重新置0
                nl++;	//统计非零个数
            }
            if(wu[l]!=0)
			{
                m_USpare.Jj[nu]=l;
                hu[nu]=j;
                ru[nu]=juo[l];
                juo[l]=nu;
                Jus[l]++;
                wu[l]=0;
                nu++;
            }
       }
       m_LSpare.Jjs[j+1]=nl;
       m_USpare.Jjs[j+1]=nu;
    }
//这里可以确定非零元素的个数了	
//实际上nl应该等于nu的
	SpareElementNum=nl>nu?(nl+1):(nu+1);

	SolveEqOnLU::Formjtoi(wide,m_USpare.Jjs,m_USpare.Jj,m_USpare.Jis,m_USpare.jtoi);
	SolveEqOnLU::Formjtoi(wide,m_LSpare.Jjs,m_LSpare.Jj,m_LSpare.Jis,m_LSpare.jtoi);

	bDecomposed=true;  //标示已经分解

	delete []rl;
	delete []ru;;
	delete []hl;
	delete []hu;
	delete []zl;
	delete []zu;
	delete []jlo;
	delete []juo;
	delete []Jls;
	delete []Jus;
	delete []wl;
	delete []wu;
}

//LU分解填充元素
bool SolveEqOnLU::SpareLUfilldata(int wide,double* constant,double* Jlowa,double* Jupa,double* JD,int* Jlowjs,int* Jupjs,int* Jlowj,int* Jupj)
{
	int i,j,k,l,m;
	double ReadSum, Sumx;
    double *midleVariable;

	midleVariable=new double[wide];

	double *unknow;
	unknow=new double[wide];

    double *Lje;
	Lje=new double[SpareElementNum];

	double *Uje;
	Uje=new double[SpareElementNum];
	
	double *Ud;
	Ud=new double[wide];

    int *zl;
	zl=new int[wide];
	
	int *zu;
	zu=new int[wide];

    double *yl;
	yl=new double[wide];
	
	double *yu;
	yu=new double[wide];

    //double* pLU;
	//LU分解()
    for(i=0;i<SpareElementNum;i++)
	{
        Lje[i]=0;
		Uje[i]=0;
	}
    for(j=0;j<wide;j++)
	{
        zl[j]=0;
		zu[j]=0;
        yl[j]=0;
		yu[j]=0;
        midleVariable[j]=0;
		unknow[j]=0;
        Ud[j]=0;
    }
    for(j=0;j<wide;j++)
	{
		Ud[j]=JD[j];
        for(l=Jlowjs[j];l<=Jlowjs[j+1]-1;l++)
		{
            k=Jlowj[l];
            yl[k]=Jlowa[l];
        }
        for(l=Jupjs[j];l<=Jupjs[j+1]-1;l++)
		{
            k=Jupj[l];
            yu[k]=Jupa[l];
        }
        yu[j]=0;
        for(l=m_LSpare.Jis[j];l<=m_LSpare.Jis[j+1]-1;l++)
		{
            k=m_LSpare.Ji[l];
            ReadSum=Lje[m_LSpare.jtoi[l]];
            for(m=m_USpare.Jjs[k]+zu[k];m<=m_USpare.Jjs[k+1]-1;m++)
			{
                i=m_USpare.Jj[m];
                yu[i]-=Uje[m]*ReadSum;
            }
            zu[k]++;
        }
        Ud[j]+=yu[j];
        for(l=m_USpare.Jjs[j];l<=m_USpare.Jjs[j+1]-1;l++)
		{
            i=m_USpare.Jj[l];
            Uje[l]=yu[i];
            yu[i]=0;
        }
        for(l=m_USpare.Jis[j];l<=m_USpare.Jis[j+1]-1;l++)
		{
            k=m_USpare.Ji[l];
            ReadSum=Uje[m_USpare.jtoi[l]];
            for(m=m_LSpare.Jjs[k]+zl[k];m<=m_LSpare.Jjs[k+1]-1;m++)
			{
                i=m_LSpare.Jj[m];
                yl[i]-=Lje[m]*ReadSum;
            }
            zl[k]++;
        }
        for(l=m_LSpare.Jjs[j];l<=m_LSpare.Jjs[j+1]-1;l++)
		{
            i=m_LSpare.Jj[l];
            Lje[l]=yl[i]/Ud[j];
            yl[i]=0;
        }
    }

	//前代回代()
	//张伯明的书写的比较清楚
    for(i=0;i<wide;i++)
       midleVariable[i]=constant[i];
    for(i=0;i<wide;i++)		//前代
	{
        for(k=m_LSpare.Jjs[i];k<=m_LSpare.Jjs[i+1]-1;k++)
		{
            j=m_LSpare.Jj[k];
            midleVariable[j]-=Lje[k]*midleVariable[i];
        }
    }
    for(i=wide-1;i>=0;i--)		//回代
	{
		Sumx=0;
        for(k=m_USpare.Jjs[i];k<=m_USpare.Jjs[i+1]-1;k++)
		{
			j=m_USpare.Jj[k];
            Sumx-=Uje[k]*unknow[j];
        }
        unknow[i]=(midleVariable[i]+Sumx)/Ud[i];
		constant[i]=unknow[i];   //以输入值的空间返回输出值
    }

	delete []midleVariable;
	delete []Lje;
	delete []Uje;
	delete []Ud;
	delete []zl;
	delete []zu;
	delete []yl;
	delete []yu;
	delete []unknow;
	return true;

}

SolveEqOnLU::~SolveEqOnLU()
{
 if(bDecomposed)
 {
	delete []m_LSpare.Jj;
	delete []m_LSpare.Jjs;
	delete []m_LSpare.Ji;
	delete []m_LSpare.Jis;
	delete []m_LSpare.jtoi;

	if(bLUflag)
	{
	delete []m_USpare.Jj;
	delete []m_USpare.Jjs;
	delete []m_USpare.Ji;
	delete []m_USpare.Jis;
	delete []m_USpare.jtoi;
	}
 }

}