matriang.cpp

电力系统稀疏矩阵计算类,用法范例代码,包括三个目录,MALIB(类),MATEST(使用范例),INCLUDE(*.H),是我编写的电力系统能量管理系统基础类之一.

//===========================================================================
//
//                            EMS高级应用软件
//
//===========================================================================
//  rev    0.0   8-9-2004      initial release   Qwbin
//  版本   0.0   8-9-2004      第一次发布        齐文斌
//---------------------------------------------------------------------------
//
//  稀疏矩阵计算-----对称矩阵分解为上下三角矩阵的积
//  数据来源: mamatrix,macmatrix 
//  结果数据: magmtrxpt
//  思路:   A=LU   U,L 分别为上/下三角阵把L+U-I 存入 magmtrxpt
//
//***************************************************************************
//matriang - triangulerize symmetric sparce matrix into upper/lower factors
//***************************************************************************

#include <malib.h>


//************************************************************************
int EmsMatrix::matriang (void)
{
   int rc = SUCCESS;        
   int i;                   
   struct workst work;         
 
   struct workst *workpt;      
  
   struct orderst   *orderpt;  
   
   workpt = &work;
   luinitialize (workpt );
   
   orderpt = workpt->orderpt;
   for ( i=0;  i<nmastate; i++,orderpt++ )
   {
      if ( orderpt->orderstat.active )
      {
         if ( orderpt->cmtrxcnt > AS_TRI_WK_ROW_SIZE )
         {
            rc = AS_TRI_WORK_ROW_OFL;  
            break;
         }

         initialize_row ( orderpt, workpt );
   
         evaluate_row ( orderpt, workpt );

         if ( workpt->nwork > AS_TRI_WK_ROW_SIZE )
         {
            rc = AS_TRI_WORK_ROW_OFL;     
            break;
         }
   
         if ( workpt->gmatrixitm + workpt->nwork > nmagmtrx)
         {
            rc = AS_TRI_GMATRIX_OFL;        
            break;
         }
         else
         {
            rc = normalize_row ( orderpt, workpt );
            if ( rc != SUCCESS )
               break;                   
   
            if ( workpt->headpt != NULL )
            {
               workpt->tailpt->linkpt = workpt->garbagept;
               workpt->garbagept      = workpt->headpt;
            }
         }
      }
   }
   
   luterminate ( workpt, rc );

   return ( rc );
}

//************************************************************************
// initialize - Initialize internal work areas
// 把计算数据放到workpt中.
//************************************************************************
void EmsMatrix::luinitialize(struct workst *workpt)
{
   int i;              
   int size;        
   int link_size;   
   int order_size;  
   int gmatrix_size;
   struct orderst *orderpt;     
   struct linkst *linkpt;      
   struct linkst *next_linkpt; 
   struct maorderst *maorderpt;      
   struct mastatest *mastatept;      
   
   memset ( workpt, 0, sizeof(struct workst) );
  
   cmtrxpt =(struct maamtrxst*)cmtrxft;
 
   workpt->singular_mastatept = divg_state_pt;

   link_size    = AS_TRI_WK_ROW_SIZE * sizeof ( struct linkst );
   order_size   = nmastate * sizeof ( struct orderst );
   gmatrix_size = nmagmtrx * sizeof ( struct gmatrixst );
   
   size         = link_size + order_size + gmatrix_size;
   workpt->size = size;
   workpt->basept =(char *)malloc ( size );
  
   memset ( workpt->basept, 0, size );
  
   workpt->linkpt = ( struct linkst * ) workpt->basept;
   workpt->garbagept = workpt->linkpt;
   linkpt            = workpt->linkpt;
   for ( i=0;  i<AS_TRI_WK_ROW_SIZE; i++ )
   {
      next_linkpt = linkpt + 1;
      linkpt->linkpt = next_linkpt;
      linkpt = next_linkpt;
   }
   linkpt--;
   linkpt->linkpt = NULL;

   workpt->orderpt = ( struct orderst * ) ( workpt->basept + link_size );

   orderpt = workpt->orderpt;
 
   for ( i=1;  i<=nmastate;i++,orderpt++)
   {
      maorderpt = maorder_vt[i-1];

	  orderpt->mastate_itm  = maorderpt->mastate_itm;

	  mastatept = mastate_vt[orderpt->mastate_itm - 1];
      mastatept->stat.partial = FALSE;

      orderpt->diag = mastatept->diag;

	  //修正PQ分解法用  mdn  04-11-11
	  orderpt->diag1 = mastatept->diag1;
	  //

      if ( mastatept->stat.active )
         orderpt->orderstat.active = TRUE;
   
      orderpt->item = i;
      orderpt->cmtrxpt = cmtrxpt + maorderpt->cmtrxitm - 1;
      orderpt->cmtrxcnt = maorderpt->cmtrxcnt;
   }

   workpt->gmatrixpt = ( struct gmatrixst * ) ( workpt->basept
   					    + link_size + order_size );
   workpt->gmatrix_pt = workpt->gmatrixpt;
   
   return;
}

//************************************************************************
// initialize_row  - Initialize working row from C-matrix coefficients
// 按行Link起来.
// 指针从大到小排列
//************************************************************************
void EmsMatrix::initialize_row (struct orderst *orderpt, struct workst *workpt)       // Shared work area ptr             
{
   int i, nwork; 
   struct   linkst    *linkpt;       
   struct   linkst    *headpt;       
   struct   linkst    *tailpt;       
   struct   maamtrxst *cmtrxpt;      
   struct   orderst   *col_orderpt;  
   struct   orderst *link_orderpt;
   struct   linkst *prev_linkpt;

   headpt = NULL;
   tailpt = NULL;
   nwork  = 0;
   
   cmtrxpt = orderpt->cmtrxpt;
   for ( i=0;i<orderpt->cmtrxcnt;i++,cmtrxpt++ )
   {
      col_orderpt = workpt->orderpt + cmtrxpt->colitm - 1;
      if ( col_orderpt->orderstat.active )
      {
         //register struct orderst *link_orderpt;
         linkpt = workpt->garbagept;
         workpt->garbagept = linkpt->linkpt;
         nwork++;
   
         linkpt->orderpt = col_orderpt;
         link_orderpt    = linkpt->orderpt;
         linkpt->element = cmtrxpt->element;

		 //支路电抗倒数  mdn  04-11-10
		 linkpt->xij = cmtrxpt->xij;
		 //
   
         if ( headpt == NULL )
         {
            headpt = linkpt;
            tailpt = linkpt;
            linkpt->linkpt = NULL;
         }
         else if ( link_orderpt < headpt->orderpt )
         {
            linkpt->linkpt = headpt;
            headpt = linkpt;
         }
         else if ( link_orderpt > tailpt->orderpt )
         {
            linkpt->linkpt = NULL;
            tailpt->linkpt = linkpt;
            tailpt         = linkpt;
         }
         else
         {
            //register struct linkst *prev_linkpt;
            prev_linkpt = headpt;
            while ( prev_linkpt->linkpt != NULL )
            {
               if ( link_orderpt < prev_linkpt->linkpt->orderpt )
               {
                  linkpt->linkpt      = prev_linkpt->linkpt;
                  prev_linkpt->linkpt = linkpt;
                  break;
               }
               prev_linkpt = prev_linkpt->linkpt;
            }
         }
      }
   }
   
   workpt->headpt = headpt;
   workpt->tailpt = tailpt;
   workpt->nwork  = nwork;
   
   return;
}

//************************************************************************
// evaluate_row  - evaluate working row elements
// 计算 -Aip*Apj和
//************************************************************************
void EmsMatrix::evaluate_row(struct orderst *orderpt, struct workst *workpt)
{
   struct   linkst    *linkpt;    
   struct   linkst    *new_linkpt;   
   struct   linkst    *search_linkpt;
   struct   linkst    *tailpt;  
   struct   linkst *prev_linkpt;
   struct orderst *gmatrix_order_rowpt;
   
   struct   gmatrixst *gmatrixpt;    
   
   int      nwork;                   
   double   product;                 
   double   zero = 0.0;    
   
   //修正PQ分解法用  mdn  04-11-10
   double   product1;
   //
  
   tailpt = workpt->tailpt;
   nwork  = workpt->nwork;
   
   linkpt = workpt->headpt;
   while ( linkpt != NULL )
   {
      search_linkpt = linkpt;
      gmatrixpt     = linkpt->orderpt->gmatrix_colpt;
      while ( gmatrixpt != NULL )
      {
         //register struct orderst *gmatrix_order_rowpt;
         gmatrix_order_rowpt = gmatrixpt->order_rowpt;
         product = linkpt->element * gmatrixpt->element;

		 //修正PQ分解法用  mdn  04-11-10
		 product1 = linkpt->xij * gmatrixpt->xij; 
		 //
          
         if ( gmatrix_order_rowpt > tailpt->orderpt )
         {
            nwork++;
            if ( nwork > AS_TRI_WK_ROW_SIZE )
               break;
            new_linkpt = workpt->garbagept;
            workpt->garbagept = new_linkpt->linkpt;
            tailpt->linkpt = new_linkpt;
            tailpt         = new_linkpt;
            new_linkpt->linkpt = NULL;
            new_linkpt->orderpt = gmatrix_order_rowpt;
            new_linkpt->element = zero - product;

			//修正PQ分解法用  mdn  04-11-10
			new_linkpt->xij = zero - product1;
			//

            search_linkpt = new_linkpt;
         }
         else
         {
            //register struct linkst *prev_linkpt;
            prev_linkpt = search_linkpt;
            while ( prev_linkpt != NULL )
            {
               if ( prev_linkpt->linkpt->orderpt == gmatrix_order_rowpt )
               {
                  prev_linkpt->linkpt->element -= product;

				  //修正PQ分解法用  mdn  04-11-10
				  prev_linkpt->linkpt->xij -= product1;
				  //

                  search_linkpt = prev_linkpt->linkpt;
                  break;
               }
          
               if ( gmatrix_order_rowpt < prev_linkpt->linkpt->orderpt )
               {
                  nwork++;
                  if ( nwork > AS_TRI_WK_ROW_SIZE )
                     break;
                  new_linkpt = workpt->garbagept;
                  workpt->garbagept = new_linkpt->linkpt;
                  new_linkpt->orderpt = gmatrix_order_rowpt;
                  new_linkpt->element = zero - product;

				  //修正PQ分解法用  mdn  04-11-10
				  new_linkpt->xij = zero - product1;
				  //

                  new_linkpt->linkpt  = prev_linkpt->linkpt;
                  prev_linkpt->linkpt = new_linkpt;
                  search_linkpt = new_linkpt;
          
                  break;
               }
               prev_linkpt = prev_linkpt->linkpt;
            }
            if ( nwork > AS_TRI_WK_ROW_SIZE )
               break;
         }
         gmatrixpt = gmatrixpt->gmatrixpt;
      }
   
      if ( nwork > AS_TRI_WK_ROW_SIZE )
         break;
       
      linkpt = linkpt->linkpt;
   }
   
   workpt->tailpt = tailpt;
   workpt->nwork  = nwork;
   
   return;
}
//************************************************************************
// normalize_row - normalize working row elements
// 
//  aij=aij-aip*apj
//
//************************************************************************
int EmsMatrix::normalize_row(struct orderst *orderpt, struct workst *workpt)
{
   struct   linkst    *linkpt;       
   struct   gmatrixst *gmatrix_pt;   
   int      rc          = SUCCESS;   
   int      gmatrixcnt  = 0;         
   double   dot_product = 0.0; 
   
   //修正PQ分解法用  mdn  04-11-10
   double   dot_product1 = 0.0;
   //
   
   orderpt->gmatrixpt = workpt->gmatrix_pt;
   gmatrix_pt = workpt->gmatrix_pt;
   linkpt     = workpt->headpt;
   while ( linkpt != NULL )
   {
      gmatrix_pt->element = linkpt->element * linkpt->orderpt->diag;

	  //修正PQ分解法用  mdn  04-11-10
	  gmatrix_pt->xij = linkpt->xij * linkpt->orderpt->diag1;
	  //

      workpt->gmatrixitm++;
      gmatrix_pt->item        = workpt->gmatrixitm;
      gmatrix_pt->order_rowpt = orderpt;
      gmatrix_pt->orderpt     = linkpt->orderpt;
      dot_product += gmatrix_pt->element * linkpt->element;

	  //修正PQ分解法用  mdn  04-11-10
	  dot_product1 += gmatrix_pt->xij * linkpt->xij;
	  //

      if ( linkpt->orderpt->gmatrix_colpt == NULL )
         linkpt->orderpt->gmatrix_colpt = gmatrix_pt;
      else
         linkpt->orderpt->gmatrix_last_colpt->gmatrixpt = gmatrix_pt;
   
      linkpt->orderpt->gmatrix_last_colpt = gmatrix_pt;
      gmatrix_pt++;
      gmatrixcnt++;
      linkpt = linkpt->linkpt;
   }
   
   orderpt->diag = orderpt->diag - dot_product;

   //修正PQ分解法用  mdn  04-11-11
   orderpt->diag1 = orderpt->diag1 - dot_product1;
   //
   orderpt->gmatrixcnt = gmatrixcnt;
   
   if ( fabs ( orderpt->diag ) < AS_TRI_TOL )
   {
      rc = AS_TRI_SINGULAR;
      workpt->singular_mastatept = orderpt->mastate_itm;
   }
   else
   {
	   orderpt->diag = 1.0 / orderpt->diag;

	   //修正PQ分解法用  mdn  04-11-11
	   orderpt->diag1 = 1.0 / orderpt->diag1;
	   //
   }
   workpt->gmatrix_pt  = gmatrix_pt;
   return ( rc );
}



//************************************************************************
// terminate - termination processing
//************************************************************************
void EmsMatrix::luterminate(struct workst *workpt, int rc)
{
   struct   maorderst  *maorderpt;       
   struct   magmtrxst *magmtrxpt;      
   struct   orderst   *orderpt;      
   struct   gmatrixst *gmatrixpt;    
   int i, j;                       
   
   if ( rc == SUCCESS )
   {
     
      orderpt = workpt->orderpt;
      for ( i=1;  i<=nmastate;  i++,orderpt++ )
      {
         maorderpt  =maorder_vt[i-1];

		 maorderpt->diag = orderpt->diag;

		 //修正PQ分解法用  mdn  04-11-11
		 maorderpt->diag1 = orderpt->diag1;
		 //

         if ( orderpt->orderstat.active )
   	      maorderpt->stat.active = TRUE;
         else
   	      maorderpt->stat.active = FALSE;
   
         if ( orderpt->gmatrixpt == NULL )
         {
            maorderpt->magmtrx_itm = 0;
            maorderpt->magmtrx_cnt = 0;
         }
         else
         {
            maorderpt->magmtrx_itm = orderpt->gmatrixpt->item;
            maorderpt->magmtrx_cnt = orderpt->gmatrixcnt;
         }
         if ( orderpt->gmatrix_colpt == NULL )
         {
            maorderpt->magmtrx_col = 0;
            maorderpt->link.magmtrx_lastcol = 0;
         }
         else
         {
            maorderpt->magmtrx_col = orderpt->gmatrix_colpt->item;
            maorderpt->link.magmtrx_lastcol = orderpt->gmatrix_last_colpt->item;
         }
   
         if ( orderpt->gmatrixpt != NULL )
         {
            gmatrixpt = orderpt->gmatrixpt;
            
            for ( j=0; j < orderpt->gmatrixcnt; j++,gmatrixpt++)
            {
               //qwbin 2006-2-18 
			   magmtrxpt = magmtrx_vt[j+orderpt->gmatrixpt->item-1];
			   //

			   magmtrxpt->rowitm = i;
               magmtrxpt->colitm = gmatrixpt->orderpt->item;
               if ( gmatrixpt->gmatrixpt == NULL )
                  magmtrxpt->magmtrx_itm = 0;
               else
                  magmtrxpt->magmtrx_itm = gmatrixpt->gmatrixpt->item;
               magmtrxpt->element = gmatrixpt->element;

			   //修正PQ分解法用  mdn  04-11-10
			   magmtrxpt->xij = gmatrixpt->xij;
			  //


            }   
         }
      }
   }
   
   free( workpt->basept );
   
   return;
}