maatrian.cpp

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

//************************************************************************
// maatrian - triangularize sparse asymmetric matrix into factors
// for opf 需要用时校对
//************************************************************************
#include <malib.h>

//************************************************************************
int EmsMatrix::maatrian (void)     
{
   int rc = SUCCESS;             
   int i;                        
   struct workst    work;          
   struct workst    *workpt;        
   struct orderst   *orderpt;       
   
   workpt = &work;
   opfinitialize (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; 
            printf("maatrian: working row overflow! \n");
            break;
         }
         opfinitialize_row ( orderpt, workpt );
   
         rc = opfevaluate_row ( orderpt, workpt );
         if ( rc != SUCCESS )
         {
            printf("maatrian: singular matrix-LF bus_itm=%d \n ",orderpt->mastate_itm);
            break;      
         }
         if ( workpt->nwork > AS_TRI_WK_ROW_SIZE )
         {
            rc = AS_TRI_WORK_ROW_OFL; 
            printf("maatrian: working row overflow \n");
            break;
         }
   
         if ( workpt->gmatrixitm + workpt->nwork > nmagmtrx )
         {
            rc = AS_TRI_GMATRIX_OFL;         
            printf("maatrian: gmatrixitm %d nwork %d > nmagmtrx %d \n",
                   workpt->gmatrixitm,workpt->nwork,nmagmtrx);
            break;
         }
         else
         {
            opfnormalize_row ( orderpt, workpt );
   
            if ( workpt->headpt != NULL )
            {
               workpt->tailpt->linkpt = workpt->new_itempt;
               workpt->new_itempt = workpt->headpt;
            }
         }
      }
   }
    
   opfterminate ( workpt, rc );
   return ( rc );
}

//************************************************************************
// initialize - initialize and build internal work areas
//************************************************************************
void EmsMatrix::opfinitialize(struct workst *workpt)
{
   int i, size;        
   int link_size, order_size;  
   int gmatrix_size;
   struct orderst *orderpt;     
   struct linkst  *linkpt,  *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->new_itempt = 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;
      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
//************************************************************************
void EmsMatrix::opfinitialize_row(struct orderst *orderpt,struct workst *workpt)
{
   int i, nwork = 0;                   
   struct linkst *linkpt, *headpt = NULL, *tailpt = NULL;       
   struct maamtrxst *cmtrxpt;      
   struct orderst *col_orderpt,*link_orderpt;  
   struct linkst *prev_linkpt;
   
   cmtrxpt = orderpt->cmtrxpt;

   for ( i=0; i<orderpt->cmtrxcnt; i++, cmtrxpt++ )
   {
      col_orderpt = workpt->orderpt + cmtrxpt->colitm - 1;
      if ( col_orderpt->orderstat.active )
      {
     
		 linkpt = workpt->new_itempt;
         workpt->new_itempt = linkpt->linkpt;
         nwork++;
          
         linkpt->orderpt = col_orderpt;
         link_orderpt    = linkpt->orderpt;
         linkpt->element = cmtrxpt->element;
   
         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
         {
    
			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
//************************************************************************
int EmsMatrix::opfevaluate_row(struct orderst *orderpt, struct workst *workpt)
{
   struct linkst *linkpt, *new_linkpt, *search_linkpt, *tailpt;       
   struct gmatrixst *gmatrixpt;   
   int nwork;                 
   double product, zero = 0.0;            
   int rc = SUCCESS; 
   int j;
   struct linkst *prev_linkpt;

   tailpt = workpt->tailpt;
   nwork  = workpt->nwork;
   linkpt = workpt->headpt;
   while ( linkpt != NULL )
   {
      search_linkpt = linkpt;
      gmatrixpt     = linkpt->orderpt->gmatrixpt;
      for ( j = 0; j < linkpt->orderpt->gmatrixcnt; j++, gmatrixpt++)
      {
         if ( gmatrixpt->orderpt != orderpt )
         {
            if ( gmatrixpt->orderpt > linkpt->orderpt )
            {
               product = linkpt->element * gmatrixpt->element;
               if ( gmatrixpt->orderpt > tailpt->orderpt )
               {
                  nwork++;
                  if ( nwork > AS_TRI_WK_ROW_SIZE )
                  {
                     printf("maatrian: Row overfull due to excess fill-in\n");
                     break;
                  }
                  new_linkpt = workpt->new_itempt;
                  workpt->new_itempt = new_linkpt->linkpt;
                  tailpt->linkpt = new_linkpt;
                  tailpt = new_linkpt;
                  new_linkpt->linkpt = NULL;
                  new_linkpt->orderpt = gmatrixpt->orderpt;
                  new_linkpt->element = zero - product;
                  search_linkpt = new_linkpt;
               }
               else
               {
                  //register struct linkst *prev_linkpt;
                  prev_linkpt = search_linkpt;
                  while ( prev_linkpt != NULL )
                  {
                     if ( prev_linkpt->linkpt->orderpt == gmatrixpt->orderpt)
                     {
                        prev_linkpt->linkpt->element -= product;
                        search_linkpt = prev_linkpt->linkpt;
                        break;
                     }
          
                     if ( gmatrixpt->orderpt < prev_linkpt->linkpt->orderpt)
                     {
                        nwork++;
                        if ( nwork > AS_TRI_WK_ROW_SIZE )
                        {
                           printf("maatrian: Row overfull due to excessive fill-in\n");
                           break;
                        }
                        new_linkpt = workpt->new_itempt;
                        workpt->new_itempt = new_linkpt->linkpt;
                        new_linkpt->orderpt = gmatrixpt->orderpt;
                        new_linkpt->element = zero - product;
                        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 )
                  {
                     printf("maatrian: Row overfull due to excessive fill-in \n");
                     break;
                  }
               }
            }
         }
         else
         {
            orderpt->diag -= linkpt->element * gmatrixpt->element;
         }
      }
   
      if ( nwork > AS_TRI_WK_ROW_SIZE )
      {
         printf("maatrian: Row overfull due to excessive fill-in\n");
         break;
      }
      linkpt = linkpt->linkpt;
   }
    
   if ( fabs ( orderpt->diag ) < AS_TRI_TOL )
   {
      rc = AS_TRI_SINGULAR;
      workpt->singular_mastatept = orderpt->mastate_itm;
   }
   else
      orderpt->diag = 1.0 / orderpt->diag;

   workpt->tailpt = tailpt;
   workpt->nwork  = nwork;
   
   return(rc);
}

//************************************************************************
// normalize_row - normalize working row elements for upper array
//************************************************************************
void EmsMatrix::opfnormalize_row(struct orderst *orderpt,struct workst *workpt)
{
   struct linkst *linkpt;       
   struct gmatrixst *gmatrix_pt; 
   int gmatrixcnt = 0;       
   double dot_product = 0.0;     
   
   orderpt->gmatrixpt = workpt->gmatrix_pt;
   gmatrix_pt = workpt->gmatrix_pt;
   linkpt = workpt->headpt;
   while ( linkpt != NULL )
   {
      if ( linkpt->orderpt > orderpt )
         gmatrix_pt->element = linkpt->element * orderpt->diag;
      else
         gmatrix_pt->element = linkpt->element * linkpt->orderpt->diag;
   
      workpt->gmatrixitm++;
      gmatrix_pt->item        = workpt->gmatrixitm;
      gmatrix_pt->order_rowpt = orderpt;
      gmatrix_pt->orderpt     = linkpt->orderpt;
   
      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->gmatrixcnt = gmatrixcnt;
   workpt->gmatrix_pt  = gmatrix_pt;
   
   return;
}

//************************************************************************
// terminate
//************************************************************************
void EmsMatrix::opfterminate(struct workst *workpt, int rc)
{
   struct maorderst  *maorderpt;       
   struct magmtrxst *magmtrxpt;      
   struct orderst   *orderpt;      
   struct gmatrixst *gmatrixpt;    
   int i, j;                       
   
   if ( rc == SUCCESS )
   {
     
	  //
	  if (magmtrxft==NULL)
	  {
	  magmtrxft = (char*) new struct magmtrxst[nmagmtrx];
	  }
	  else
	  {
		  magmtrx_vt.clear();
	  }
	  memset(magmtrxft,0,nmagmtrx* sizeof(struct magmtrxst));
      magmtrxpt = (struct magmtrxst *)magmtrxft;
	  // 
	
	  orderpt = workpt->orderpt;
      for ( i=1;  i<=nmastate; i++,orderpt++ )
      {
         maorderpt = maorder_vt[i-1];
		 maorderpt->diag = orderpt->diag;
         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;
            //magmtrxpt = magmtrxpt + orderpt->gmatrixpt->item - 1;
            for ( j=0; j < orderpt->gmatrixcnt; j++, gmatrixpt++, magmtrxpt++ )
            {
               //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;
               //              
			   magmtrx_vt.push_back(magmtrxpt);
               //
            }
         }
      }
   }
   
   free( workpt->basept );

   return;
}