sparsematrix.h

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

#ifndef SPARSE_MATRIX_H
#define SPARSE_MATRIX_H

#include <cassert>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>

#include "Defs.h"
#include "Complex.h"

using namespace std;

//********************************************
//2008.5.29
//对稀疏矩阵，增加了备份矩阵元素的功能
//这样，因子表化后仍可以获得原来的矩阵元素
//********************************************

template < typename NumericType > class SparseMatrix;

// 矩阵行元素节点
template < typename Type > class CMNode {
  friend class SparseMatrix<Type>;

private:
	CMNode();
	CMNode(Type const& v);

	int col;
	Type value;
	CMNode<Type> * right;
};

// 矩阵行元素表头节点
template < typename Type > class CMtNode {
  friend class SparseMatrix<Type>;

private:
  CMtNode ();

  int nnz_row;                // 每行非零元素数
  CMNode<Type>  * right;      // 行元素链表指针
};

// 稀疏矩阵
template < typename Type > class SparseMatrix {
public:
  SparseMatrix();
  SparseMatrix(int rows, int cols);
  ~SparseMatrix();

	bool Initialize(int rows,int cols);

	Type GetElement(int r, int c);
  void SetElement(int r, int c, Type const& v);
	bool DeleteElement(int r, int c);

  int GetNum();                       //得到矩阵内非零元素个数
	int GetRNum(int r);                 //得到矩阵每行非零元素个数
	int GetLANum(int r);		            //获得下三角的元素个数
  void SetRNum(int r,int nums);       //设置每行非零元素个数

	int GetRows();                      //得到矩阵行数
	int GetCols();                      //得到矩阵列数

	inline	bool GetFactorFlag();
	bool BuildFactorTable();
	bool ResolveEquation(Type* pB);

	void Display();
	void DisplayInFile(char *pFileName);
	void DisplayMatlab(char *pFileName);

	bool IsBackuped() { return bBackuped; }
	bool BackupElements();                  //备份稀疏矩阵的元素
	bool RestoreBackupedElems();            //恢复成备份的状态
	bool DelBackupedElements();             //删除备份的元素
	void DeleteElements();                  //清空已有元素
	bool DeleteAllElements();               //清空存储的内容

private:
  template < typename T > 
  T Abs(T const& v) { return std::abs(v); }
  double Abs(Complex const& c) { return c.absCplx(); }

private:
	int rows,cols;        //矩阵的大小 

	CMtNode<Type> * H;    //指向行链表数组头 
	CMtNode<Type> * H0;   //因为因子表化后原来的矩阵会改变，用此指针来保存原来的元素
	
  bool bInitialized;    //是否初始化了
	bool bFactorBuilded;   //是否因子表化了
	bool bBackuped;       //是否对稀疏矩阵的元素做了备份 
  
  Type ** FullFactorTable;
};

// 矩阵行元素节点构造函数
template < class Type > CMNode<Type>::CMNode(Type const& v) {
   right = NULL;
   value = v;
   col = -1;
}

// 矩阵行元素表头节点构造函数
template < typename Type > CMtNode<Type>::CMtNode() {
	right = NULL;
	nnz_row = 0;		  
}

template < typename Type > CMNode<Type>::CMNode() {
	right = NULL;
	value = 0;
	col = -1;
}

// 稀疏矩阵构造函数
template < typename Type > SparseMatrix<Type>::SparseMatrix()
: bInitialized(false), bFactorBuilded(false), bBackuped(false)
{}

template < typename Type > SparseMatrix<Type>::SparseMatrix(int r, int c) {
	Initialize(r,c);
}

// 稀疏矩阵析构函数
template < typename Type > SparseMatrix<Type>::~SparseMatrix() {
  if(!bInitialized) return;

  for ( int I = 0; I < rows; ++I ) {
    CMNode<Type> * p = H[I].right;
    do {
      if ( p == NULL ) break;
      CMNode<Type> *pt = p;
      p = p->right;
      delete pt;
    } while ( true );
  }
  delete []H;

  if ( bBackuped ) {
    for ( int I = 0; I < rows; ++I ) {
      CMNode<Type> * p;
      p = H0[I].right;
      do {
        if ( p == NULL ) break;
        CMNode<Type> *pt = p;
        p = p->right;
        delete pt;
      } while ( true );
    }
    delete []H0;	
  }

  if ( bFactorBuilded ) {
	   for ( int i = 0; i < rows; ++i )
       delete [] FullFactorTable[i];
  }
}

// 稀疏矩阵初始化辅助函数
template < typename Type > bool SparseMatrix <Type>::Initialize(int r, int c) {
  rows = r;
	cols = c;
	H = new CMtNode<Type>[r]; //初始化行表头针数组指，指向各链表头,N行共N个表头；
  bFactorBuilded = false;
	bInitialized = true;
	bBackuped = false;
	return true;
}

template < typename Type > Type SparseMatrix <Type>::GetElement(int r, int c) { 
	CMNode <Type> *p = H[r].right;
	do {
	   if ( p == NULL ) return Type(0);   //没找到，返回零
	   if ( p->col > c ) return Type(0);  //没找到，返回零
	   if ( p->col == c ) return Type(p->value);

     p = p->right;
	} while ( true );
}


template < typename Type > bool SparseMatrix<Type>::DeleteElement(int r,int c) {
	CMNode <Type> *p1 = H[r].right;
	if ( p1 == NULL ) return false;
  if ( p1->col == c ) {
		H[r].right = p1->right;
		delete p1;
		return true;
	}
	CMNode <Type> *p2 = p1->right;
	do {
		if ( p2 == NULL ) return false;
		if( p2->col == c ) {
			p1->right = p2->right;
			H[r].nnz_row--;
			delete p2;
			break;
		}
		p1 = p2;
		p2 = p2->right;
	} while ( true );
	return true;
}

//-----------------------SetEllement
template < typename Type > void SparseMatrix<Type>::SetElement(int r, int c,Type const& v){
  assert( r < rows && c < cols && r >= 0 && c >= 0 );

  if(Abs(v)<THRESHOLD){
    DeleteElement(r,c);
    return;
  }

  CMNode <Type> *p2 ;
  CMNode <Type> *p1 ;
  p1 = H[r].right;
  if(p1 == NULL){
    CMNode<Type> * p3 = new CMNode<Type>(v);
    p3->col = c;
    H[r].right = p3;
    H[r].nnz_row++;
    return;
  }                             //如果该行无元素，加于末尾
  do{
    if(p1 == NULL){
      CMNode <Type> *p3 = new CMNode<Type>(v);
      p3->col = c;
      p2->right =p3;    //无此元素，且加入此元素至链表尾
      H[r].nnz_row++;
      return;
    }
    int k = p1->col - c;
    if(k<0){                       //继续找
      p2 = p1;
      p1 = p1->right;
    }

    else if(k == 0){
      p1->value = v;   //找到该元素，直接替换   
      return;
    }
    else if(k > 0){                //无则在p1和p2之间插入一元素
      CMNode <Type> *p3 = new CMNode<Type>(v);
      p3->col = c;
      if(H[r].right == p1){  //第一个元素
        H[r].right = p3;
        p3->right = p1;
      }
      else{
        p2->right = p3;
        p3->right = p1;
      }
      H[r].nnz_row++;
      return;
    }                   
  }while(true);

}

//-------------------------------得到矩阵非零元素个数
template < typename Type >  int SparseMatrix<Type>::GetNum(){
	int num = 0;
	for(int i = 0;i < rows; i++)num += H[i].nnz_row;
	return num;
};

//------------------------------得到矩阵每行非零元素个数
template < typename Type > int SparseMatrix<Type>::GetRNum(int r){
	return H[r].nnz_row;
};
//-----------------------------//获得行号大于列号，处于下三角的元素个数
template < typename Type > int SparseMatrix<Type>::GetLANum(int r){
	int counter=0;
	CMNode <Type> *p1;
	p1=H[r].right;
	while(p1)
	{
		if(p1->col<r)
		{
			counter++;
			p1=p1->right;
		}
		else
			break;
	}
	return counter;
};

//-----------------------------设置矩阵每行非零元素个数
template < typename Type > void SparseMatrix<Type>::SetRNum(int r,int nums){
	H[r].nnz_row = nums;
}

//-------------------------------得到矩阵行数
template < typename Type > int SparseMatrix<Type>::GetRows(){
	return rows;
};

//-------------------------------得到矩阵列数
template < typename Type > int SparseMatrix<Type>::GetCols(){
	return cols;
};
//-------------------------------得到factorbuild的值
template < typename Type > 
bool SparseMatrix<Type>::GetFactorFlag() {
	return bFactorBuilded;
}

/*!!!!!!!!!!!原方案的备份
//------------------------------建立因子表：BuildFactors()
template < typename Type > bool CsparMatrix<Type>::BuildFactorTable()
{
	if(rows!=cols) return false; 
 int begincol = 0;
 int i,j,k;                                           
 
 for(i=0;i<rows;i++)                                  
 { 
 	Type aii=GetElement(i,i);                         
 	for(j=i+1;j<cols;j++)	//归一化                    
 	{                                                   
 	Type temp;                                
 		temp=GetElement(i,j);                       
 		if(myabs(temp)<SMALLEST) continue;	    
 		temp/=aii;                                  
		SetElement(i,j,temp);
 	}                                                   
                                                      
 	for(j=i+1;j<rows;j++)	//消去                      
 	{                                                   
 		Type aji=GetElement(j,i);                 
 		if(myabs(aji)<SMALLEST) continue;            
 		for(k=i+1;k<rows;k++)                       
 		{       
			Type ajk;
 			Type aik=GetElement(i,k);         
 			if(myabs(aik)<SMALLEST) continue;         
 			ajk=GetElement(j,k);         
			ajk=ajk-aik*aji; 
 				SetElement(j,k,ajk);                

 		}                                           
 	}                                                   
 }                                                    
bFactorBuilded=true;                                   
 return true;                              
}
*/
//------------------------------建立因子表：BuildFactors()
template < typename Type > bool SparseMatrix<Type>::BuildFactorTable()
{
	if(rows!=cols) return false; 
 int begincol = 0;